Merge pull request #11 from kenohori/cgal

Bugfixes in placements, half space solids and validation code

src/ifcgeom/kernels/cgal/CgalConversionFunctions.cpp

@@ -25,18 +25,11 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcRepresentation* l, Convers
 
 bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcCartesianPoint* l, cgal_point_t& point) {
   std::vector<double> xyz = l->Coordinates();
-  if (xyz.size() < 4) {
-    point = Kernel::Point_3(xyz.size()     ? (xyz[0]*getValue(GV_LENGTH_UNIT)) : 0.0f,
-                            xyz.size() > 1 ? (xyz[1]*getValue(GV_LENGTH_UNIT)) : 0.0f,
-                            xyz.size() > 2 ? (xyz[2]*getValue(GV_LENGTH_UNIT)) : 0.0f);
-//    std::cout << "Converted Point(" << point << ")" << std::endl;
-    return true;
-  } else {
-    std::cout << "Point(";
-    for (auto &coordinate: xyz) std::cout << coordinate << " ";
-    std::cout << ")";
-    throw std::runtime_error("Could not parse point");
-  }
+  point = Kernel::Point_3(xyz.size()     ? (xyz[0]*getValue(GV_LENGTH_UNIT)) : 0.0f,
+                          xyz.size() > 1 ? (xyz[1]*getValue(GV_LENGTH_UNIT)) : 0.0f,
+                          xyz.size() > 2 ? (xyz[2]*getValue(GV_LENGTH_UNIT)) : 0.0f);
+//  std::cout << "Converted Point(" << point << ")" << std::endl;
+  return true;
 }
 
 bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcDirection* l, cgal_direction_t& dir) {
@@ -63,15 +56,46 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcPlane* pln, cgal_plane_t&
   IfcSchema::IfcAxis2Placement3D* l = pln->Position();
   cgal_point_t o;
   cgal_direction_t axis = Kernel::Vector_3(0,0,1);
-  cgal_direction_t refDirection;
+  cgal_direction_t refDirection = Kernel::Vector_3(1,0,0);
   IfcGeom::CgalKernel::convert(l->Location(),o);
   bool hasRef = l->hasRefDirection();
   if ( l->hasAxis() ) IfcGeom::CgalKernel::convert(l->Axis(),axis);
   if ( hasRef ) IfcGeom::CgalKernel::convert(l->RefDirection(),refDirection);
+  Kernel::Vector_3 y = CGAL::cross_product(axis, refDirection);
+  Kernel::Vector_3 x = CGAL::cross_product(y, axis);
+  
   cgal_plane_t ax3;
-  if ( hasRef ) ax3 = Kernel::Plane_3(o,o+axis,o+refDirection);
+  if ( hasRef ) ax3 = Kernel::Plane_3(o,o+x,o+y);
   else ax3 = Kernel::Plane_3(o,axis);
   plane = ax3;
+  
+//  std::cout << "IfcPlane C = " << o << std::endl;
+//  std::cout << "IfcPlane z (axis, exact) = " << axis << std::endl;
+//  std::cout << "IfcPlane x (refDirection, approximate) = " << refDirection << std::endl;
+//  std::cout << "IfcPlane y (computed, exact) = " << y << std::endl;
+//  std::cout << "IfcPlane x (computed, exact) = " << x << std::endl;
+//  
+//  std::cout << "Plane_3 o = " << o << std::endl;
+//  std::cout << "Plane_3 o+x = " << o+x << std::endl;
+//  std::cout << "Plane_3 o+y = " << o+y << std::endl;
+  
+  // ax + by + cz + d = 0
+//  std::cout << "Plane: a = " << plane.a() << ", b = " << plane.b() << ", c = " << plane.c() << ", d = " << plane.d() << std::endl;
+  
+//  std::ofstream fresult;
+//  fresult.open("/Users/ken/Desktop/plane.obj");
+//  // x = -5, y = -5, z = (5a +5b -d)/c
+//  fresult << "v -5 -5 " << (5.0*CGAL::to_double(plane.a())+5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//  // x = -5, y = +5, z = (5a -5b -d)/c
+//  fresult << "v -5 5 " << (5.0*CGAL::to_double(plane.a())-5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//  // x = 5, y = -5, z = (-5a +5b -d)/c
+//  fresult << "v 5 -5 " << (-5.0*CGAL::to_double(plane.a())+5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//  // x = 5, y = +5, z = (-5a -5b -d)/c
+//  fresult << "v 5 5 " << (-5.0*CGAL::to_double(plane.a())-5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//  fresult << "f 1 2 3" << std::endl;
+//  fresult << "f 4 3 2" << std::endl;
+//  fresult.close();
+  
 //  CACHE(IfcPlane,pln,plane)
   return true;
 }
@@ -79,14 +103,13 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcPlane* pln, cgal_plane_t&
 bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcAxis2Placement2D* l, cgal_placement_t& trsf) {
   //  IN_CACHE(IfcAxis2Placement3D,l,gp_Trsf,trsf)
   cgal_point_t o;
-  cgal_direction_t axis = Kernel::Vector_3(0,0,1);
   cgal_direction_t refDirection = Kernel::Vector_3(1,0,0);
   IfcGeom::CgalKernel::convert(l->Location(),o);
   bool hasRef = l->hasRefDirection();
   if ( hasRef ) IfcGeom::CgalKernel::convert(l->RefDirection(),refDirection);
+  cgal_direction_t y = Kernel::Vector_3(-refDirection.y(), refDirection.x(), 0.0);
   
-  // TODO: From Thomas' email. Should be checked.
-  Kernel::Vector_3 y = CGAL::cross_product(Kernel::Vector_3(0.0, 0.0, 1.0), refDirection);
+  // TODO: Should be checked.
   trsf = Kernel::Aff_transformation_3(refDirection.cartesian(0), y.cartesian(0), 0.0, o.cartesian(0),
                                       refDirection.cartesian(1), y.cartesian(1), 0.0, o.cartesian(1),
                                       0.0, 0.0, 1.0, 0.0);
@@ -104,16 +127,17 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcAxis2Placement3D* l, cgal_
   bool hasRef = l->hasRefDirection();
   if ( l->hasAxis() ) IfcGeom::CgalKernel::convert(l->Axis(),axis);
   if ( hasRef ) IfcGeom::CgalKernel::convert(l->RefDirection(),refDirection);
+  Kernel::Vector_3 y = CGAL::cross_product(axis, refDirection);
+  Kernel::Vector_3 x = CGAL::cross_product(y, axis);
   
 //  std::cout << "Ref direction: " << refDirection << std::endl;
 //  std::cout << "Axis: " << axis << std::endl;
 //  std::cout << "Origin: " << o << std::endl;
   
-  // TODO: From Thomas' email. Should be checked.
-  Kernel::Vector_3 y = CGAL::cross_product(axis, refDirection);
-  trsf = Kernel::Aff_transformation_3(refDirection.cartesian(0), y.cartesian(0), axis.cartesian(0), o.cartesian(0),
-                                      refDirection.cartesian(1), y.cartesian(1), axis.cartesian(1), o.cartesian(1),
-                                      refDirection.cartesian(2), y.cartesian(2), axis.cartesian(2), o.cartesian(2));
+  // TODO: Should be checked.
+  trsf = Kernel::Aff_transformation_3(x.cartesian(0), y.cartesian(0), axis.cartesian(0), o.cartesian(0),
+                                      x.cartesian(1), y.cartesian(1), axis.cartesian(1), o.cartesian(1),
+                                      x.cartesian(2), y.cartesian(2), axis.cartesian(2), o.cartesian(2));
   
 //  for (int i = 0; i < 3; ++i) {
 //    for (int j = 0; j < 4; ++j) {
@@ -180,6 +204,54 @@ bool IfcGeom::CgalKernel::convert_wire_to_face(const cgal_wire_t& wire, cgal_fac
   return true;
 }
 
+bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcCartesianTransformationOperator2D* l, cgal_placement_t& trsf) {
+//  IN_CACHE(IfcCartesianTransformationOperator2D,l,cgal_placement_t,trsf)
+  
+  cgal_point_t origin;
+  cgal_direction_t axis1 (1.,0.,0.);
+  cgal_direction_t axis2 (0.,1.,0.);
+  
+  IfcGeom::CgalKernel::convert(l->LocalOrigin(),origin);
+  if ( l->hasAxis1() ) IfcGeom::CgalKernel::convert(l->Axis1(),axis1);
+  if ( l->hasAxis2() ) IfcGeom::CgalKernel::convert(l->Axis2(),axis2);
+  double scale = 1.0;
+  if (l->hasScale()) {
+    scale = l->Scale();
+  }
+  
+  // TODO: Untested
+  trsf = Kernel::Aff_transformation_3(scale*axis1.cartesian(0), axis2.cartesian(0), 0.0, origin.cartesian(0),
+                                      axis1.cartesian(1), scale*axis2.cartesian(1), 0.0, origin.cartesian(1),
+                                      0.0, 0.0, 1.0, 0.0);
+  
+//  CACHE(IfcCartesianTransformationOperator2D,l,trsf)
+  return true;
+}
+
+bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcCartesianTransformationOperator2DnonUniform* l, cgal_placement_t& gtrsf) {
+//  IN_CACHE(IfcCartesianTransformationOperator2DnonUniform,l,cgal_placement_t,gtrsf)
+  
+  cgal_placement_t trsf;
+  cgal_point_t origin;
+  cgal_direction_t axis1 (1.,0.,0.);
+  cgal_direction_t axis2 (0.,1.,0.);
+  
+  IfcGeom::CgalKernel::convert(l->LocalOrigin(),origin);
+  if ( l->hasAxis1() ) IfcGeom::CgalKernel::convert(l->Axis1(),axis1);
+  if ( l->hasAxis2() ) IfcGeom::CgalKernel::convert(l->Axis2(),axis2);
+  
+  const double scale1 = l->hasScale() ? l->Scale() : 1.0f;
+  const double scale2 = l->hasScale2() ? l->Scale2() : scale1;
+  
+  // TODO: Untested
+  trsf = Kernel::Aff_transformation_3(scale1*axis1.cartesian(0), axis2.cartesian(0), 0.0, origin.cartesian(0),
+                                      axis1.cartesian(1), scale2*axis2.cartesian(1), 0.0, origin.cartesian(1),
+                                      0.0, 0.0, 1.0, 0.0);
+  
+//  CACHE(IfcCartesianTransformationOperator2DnonUniform,l,gtrsf)
+  return true;
+}
+
 bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcCartesianTransformationOperator3D* l, cgal_placement_t& trsf) {
 //  IN_CACHE(IfcCartesianTransformationOperator3D,l,gp_Trsf,trsf)
   cgal_point_t origin;
@@ -257,3 +329,37 @@ void IfcGeom::CgalKernel::remove_duplicate_points_from_loop(cgal_wire_t& polygon
     }
   }
 }
+
+CGAL::Nef_polyhedron_3<Kernel> IfcGeom::CgalKernel::create_nef_polyhedron(std::list<cgal_face_t> &face_list) {
+  
+  // Naive creation
+  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
+  PolyhedronBuilder builder(&face_list);
+  polyhedron.delegate(builder);
+  
+  // Stitch edges
+  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
+  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
+  if (!polyhedron.is_valid()) {
+    std::cout << "create_nef_polyhedron: Polyhedron not valid!" << std::endl;
+//    std::ofstream fresult;
+//    fresult.open("/Users/ken/Desktop/invalid.off");
+//    fresult << polyhedron << std::endl;
+//    fresult.close();
+    return CGAL::Nef_polyhedron_3<Kernel>();
+  } if (!polyhedron.is_closed()) {
+    std::cout << "create_nef_polyhedron: Polyhedron not closed" << std::endl;
+//    std::ofstream fresult;
+//    fresult.open("/Users/ken/Desktop/open.off");
+//    fresult << polyhedron << std::endl;
+//    fresult.close();
+    // TODO: Nef constructor doesn't support open meshes
+    return CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  }
+  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
+    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
+  }
+  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
+  
+  return CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+}

src/ifcgeom/kernels/cgal/CgalEntityMapping.h

@@ -85,5 +85,7 @@ CLASS(IfcPlane,cgal_plane_t);
 CLASS(IfcAxis2Placement2D,cgal_placement_t);
 CLASS(IfcAxis2Placement3D,cgal_placement_t);
 CLASS(IfcObjectPlacement,cgal_placement_t);
-CLASS(IfcCartesianTransformationOperator3D,cgal_placement_t);
+CLASS(IfcCartesianTransformationOperator2DnonUniform,cgal_placement_t);
 CLASS(IfcCartesianTransformationOperator3DnonUniform,cgal_placement_t);
+CLASS(IfcCartesianTransformationOperator2D,cgal_placement_t);
+CLASS(IfcCartesianTransformationOperator3D,cgal_placement_t);

src/ifcgeom/kernels/cgal/CgalIfcGeomShapes.cpp

@@ -133,17 +133,6 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcExtrudedAreaSolid *l, cgal
   std::list<cgal_face_t> face_list;
   face_list.push_back(bottom_face);
   
-//  if (true) {
-//    CGAL::Polyhedron_3<Kernel> polyhedron;
-//    PolyhedronBuilder builder(&face_list);
-//    polyhedron.delegate(builder);
-//    
-//    std::ofstream fresult;
-//    fresult.open("/Users/ken/Desktop/profile.off");
-//    fresult << polyhedron << std::endl;
-//    fresult.close();
-//  }
-  
   for (std::vector<Kernel::Point_3>::const_iterator current_vertex = bottom_face.outer.begin();
        current_vertex != bottom_face.outer.end();
        ++current_vertex) {
@@ -166,33 +155,7 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcExtrudedAreaSolid *l, cgal
     top_face.outer.push_back(*vertex+height*dir);
   } face_list.push_back(top_face);
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polyhedron_3<Kernel> old_polyhedron(polyhedron);
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!polyhedron.is_valid()) {
-    std::cout << "Invalid polyhedron!" << std::endl;
-    std::ofstream fresult;
-    fresult.open("/Users/ken/Desktop/invalid.off");
-    fresult << old_polyhedron << std::endl;
-    fresult.close();
-  }
-  
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  } CGAL_postcondition(polyhedron.is_valid() && polyhedron.is_closed());
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-
-  for (auto &vertex : vertices(polyhedron)) {
-	  vertex->point() = vertex->point().transform(trsf);
-  }
-  
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
   
   // Inner
   // TODO: Would be faster to triangulate top/bottom face template rather than use Nef polyhedra for subtraction
@@ -226,48 +189,10 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcExtrudedAreaSolid *l, cgal
       hole_top_face.outer.push_back(*vertex+height*dir);
     } face_list.push_back(hole_top_face);
     
-    // Naive creation
-    CGAL::Polyhedron_3<Kernel> hole_polyhedron = CGAL::Polyhedron_3<Kernel>();
-    PolyhedronBuilder builder(&face_list);
-    hole_polyhedron.delegate(builder);
-    
-    // Stitch edges
-    //  std::cout << "Before: " << hole_polyhedron.size_of_vertices() << " vertices and " << hole_polyhedron.size_of_facets() << " facets" << std::endl;
-    CGAL::Polygon_mesh_processing::stitch_borders(hole_polyhedron);
-    if (!hole_polyhedron.is_valid()) {
-//      std::cout << "Invalid hole polyhedron!" << std::endl;
-//      std::ofstream fresult;
-//      fresult.open("/Users/ken/Desktop/invalid.off");
-//      fresult << hole_polyhedron << std::endl;
-//      fresult.close();
-      return false;
-    }
-
-	for (auto &vertex : vertices(hole_polyhedron)) {
-		vertex->point() = vertex->point().transform(trsf);
-	}
-    
-    if (!CGAL::Polygon_mesh_processing::is_outward_oriented(hole_polyhedron)) {
-      CGAL::Polygon_mesh_processing::reverse_face_orientations(hole_polyhedron);
-    } CGAL_postcondition(hole_polyhedron.is_valid() && hole_polyhedron.is_closed());
-    //  std::cout << "After: " << hole_polyhedron.size_of_vertices() << " vertices and " << hole_polyhedron.size_of_facets() << " facets" << std::endl;
-    
-    shape -= CGAL::Nef_polyhedron_3<Kernel>(hole_polyhedron);
+    shape -= create_nef_polyhedron(face_list);
   }
   
-//  std::cout << "trsf" << std::endl;
-//  for (int i = 0; i < 3; ++i) {
-//    for (int j = 0; j < 4; ++j) {
-//      std::cout << trsf.cartesian(i, j) << " ";
-//    } std::cout << std::endl;
-//  }
-  
-//  shape.convert_to_polyhedron(polyhedron);
-//  std::ofstream fresult;
-//  fresult.open("/Users/ken/Desktop/extrusion.off");
-//  fresult << polyhedron << std::endl;
-//  fresult.close();
-  
+  shape.transform(trsf);
   return true;
 }
 
@@ -296,20 +221,7 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcConnectedFaceSet* l, cgal_
     face_list.push_back(face);
   }
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  }
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
   return true;
 }
 
@@ -366,27 +278,12 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcBlock* l, cgal_shape_t& sh
   face_list.back().outer.push_back(Kernel::Point_3(dx, dy, dz));
   face_list.back().outer.push_back(Kernel::Point_3(dx, 0, dz));
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  } CGAL_postcondition(polyhedron.is_valid() && polyhedron.is_closed());
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  
   cgal_placement_t trsf;
   IfcGeom::CgalKernel::convert(l->Position(),trsf);
 
-  for (auto &vertex: vertices(polyhedron)) {
-    vertex->point() = vertex->point().transform(trsf);
-  }
-
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
+  shape.transform(trsf);
+  
   return true;
 }
 
@@ -405,11 +302,9 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcBooleanResult* l, cgal_sha
       return false;
 //    }
   } else if ( shape_type(operand1) == ST_SHAPE ) {
-    if ( ! convert_shape(operand1, s1) ) {
+    if (!convert_shape(operand1, s1) ) {
       return false;
     }
-//    TopoDS_Solid temp_solid;
-//    s1 = ensure_fit_for_subtraction(s1, temp_solid);
   } else {
     Logger::Message(Logger::LOG_ERROR, "s1: Invalid representation item for boolean operation", operand1->entity);
     return false;
@@ -425,18 +320,14 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcBooleanResult* l, cgal_sha
 //    shape2_processed = convert_shapes(operand2, items2) && flatten_shape_list(items2, s2, true);
   } else if ( shape_type(operand2) == ST_SHAPE ) {
     shape2_processed = convert_shape(operand2,s2);
-    if (shape2_processed && !is_halfspace) {
-//      TopoDS_Solid temp_solid;
-//      s2 = ensure_fit_for_subtraction(s2, temp_solid);
-    }
   } else {
     Logger::Message(Logger::LOG_ERROR, "s2: Invalid representation item for boolean operation", operand2->entity);
   }
 
   if (!shape2_processed) {
-//    shape = s1;
+    shape = s1;
     Logger::Message(Logger::LOG_ERROR,"Failed to convert SecondOperand of:",l->entity);
-//    return true;
+    return true;
   }
 
 //  if (!is_halfspace) {
@@ -450,37 +341,72 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcBooleanResult* l, cgal_sha
   if (!s1.is_simple()) {
     Logger::Message(Logger::LOG_ERROR, "s1: Not simple Nef?", operand1->entity);
     return false;
+  } else {
+//    std::ofstream f1;
+//    CGAL::Polyhedron_3<Kernel> p1;
+//    s1.convert_to_Polyhedron(p1);
+//    f1.open("/Users/ken/Desktop/s1.off");
+//    f1 << p1 << std::endl;
+//    f1.close();
   }
   
+  bool is_plane = false;
+  cgal_plane_t plane;
   if (!s2.is_simple()) {
     Logger::Message(Logger::LOG_ERROR, "s2: Not simple Nef?", operand2->entity);
     return false;
+  } else if (is_halfspace) {
+    //    std::cout << "s2: halfspace" << std::endl;
+    IfcSchema::IfcHalfSpaceSolid *hss = static_cast<IfcSchema::IfcHalfSpaceSolid *>(operand2);
+    IfcSchema::IfcSurface* surface = hss->BaseSurface();
+    if (surface->is(IfcSchema::Type::IfcPlane) ) {
+      is_plane = true;
+      IfcGeom::CgalKernel::convert((IfcSchema::IfcPlane *)surface, plane);
+      if (hss->AgreementFlag()) plane = plane.opposite();
+//      std::ofstream fresult;
+//      fresult.open("/Users/ken/Desktop/s2.off");
+//      fresult << "OFF" << std::endl << "4 2 4" << std::endl;
+//      // x = -5, y = -5, z = (5a +5b -d)/c
+//      fresult << "-5 -5 " << (5.0*CGAL::to_double(plane.a())+5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//      // x = -5, y = +5, z = (5a -5b -d)/c
+//      fresult << "-5 5 " << (5.0*CGAL::to_double(plane.a())-5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//      // x = 5, y = -5, z = (-5a +5b -d)/c
+//      fresult << "5 -5 " << (-5.0*CGAL::to_double(plane.a())+5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//      // x = 5, y = +5, z = (-5a -5b -d)/c
+//      fresult << "5 5 " << (-5.0*CGAL::to_double(plane.a())-5.0*CGAL::to_double(plane.b())-CGAL::to_double(plane.d()))/CGAL::to_double(plane.c()) << std::endl;
+//      fresult << "3 0 1 2" << std::endl;
+//      fresult << "3 3 2 1" << std::endl;
+//      fresult.close();
+    }
+  } else {
+//    std::ofstream f2;
+//    CGAL::Polyhedron_3<Kernel> p2;
+//    s2.convert_to_Polyhedron(p2);
+//    f2.open("/Users/ken/Desktop/s2.off");
+//    f2 << p2 << std::endl;
+//    f2.close();
   }
   
-//  std::ofstream f1;
-//  f1.open("/Users/ken/Desktop/s1.off");
-//  f1 << s1 << std::endl;
-//  f1.close();
-//  std::ofstream f2;
-//  f2.open("/Users/ken/Desktop/s2.off");
-//  f2 << s2 << std::endl;
-//  f2.close();
-  
   if (op == IfcSchema::IfcBooleanOperator::IfcBooleanOperator_DIFFERENCE) {
     
 //    std::cout << "Difference" << std::endl;
-    CGAL::Nef_polyhedron_3<Kernel> nef_result = s1-s2;
+    CGAL::Nef_polyhedron_3<Kernel> nef_result = s1;
+    if (is_halfspace) {
+      if (is_plane) nef_result = nef_result.intersection(plane, CGAL::Nef_polyhedron_3<Kernel>::Intersection_mode::CLOSED_HALFSPACE);
+    } else {
+      nef_result -= s2;
+    }
     if (!nef_result.is_simple()) {
       std::cout << "Not simple: " << nef_result.number_of_volumes() << " volumes" << std::endl;
       return false;
-    }
-//    cgal_shape_t result;
-//    nef_result.convert_to_polyhedron(result);
-//    std::ofstream fresult;
-//    fresult.open("/Users/ken/Desktop/result.off");
-//    fresult << result << std::endl;
-//    fresult.close();
-    shape = nef_result;
+    } else {
+//      CGAL::Polyhedron_3<Kernel> result;
+//      nef_result.convert_to_polyhedron(result);
+//      std::ofstream fresult;
+//      fresult.open("/Users/ken/Desktop/result.off");
+//      fresult << result << std::endl;
+//      fresult.close();
+    } shape = nef_result;
     return true;
     
   } else if (op == IfcSchema::IfcBooleanOperator::IfcBooleanOperator_UNION) {
@@ -490,14 +416,14 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcBooleanResult* l, cgal_sha
     if (!nef_result.is_simple()) {
       std::cout << "Not simple: " << nef_result.number_of_volumes() << " volumes" << std::endl;
       return false;
-    }
-//    cgal_shape_t result;
-//    nef_result.convert_to_polyhedron(result);
-//    std::ofstream fresult;
-//    fresult.open("/Users/ken/Desktop/result.off");
-//    fresult << result << std::endl;
-//    fresult.close();
-    shape = nef_result;
+    } else {
+//      CGAL::Polyhedron_3<Kernel> result;
+//      nef_result.convert_to_polyhedron(result);
+//      std::ofstream fresult;
+//      fresult.open("/Users/ken/Desktop/result.off");
+//      fresult << result << std::endl;
+//      fresult.close();
+    } shape = nef_result;
     return true;
     
   } else if (op == IfcSchema::IfcBooleanOperator::IfcBooleanOperator_INTERSECTION) {
@@ -507,19 +433,16 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcBooleanResult* l, cgal_sha
     if (!nef_result.is_simple()) {
       std::cout << "Not simple: " << nef_result.number_of_volumes() << " volumes" << std::endl;
       return false;
-    }
-//    cgal_shape_t result;
-//    nef_result.convert_to_polyhedron(result);
-//    std::ofstream fresult;
-//    fresult.open("/Users/ken/Desktop/result.off");
-//    fresult << result << std::endl;
-//    fresult.close();
-    shape = nef_result;
+    } else {
+//      CGAL::Polyhedron_3<Kernel> result;
+//      nef_result.convert_to_polyhedron(result);
+//      std::ofstream fresult;
+//      fresult.open("/Users/ken/Desktop/result.off");
+//      fresult << result << std::endl;
+//      fresult.close();
+    } shape = nef_result;
     return true;
-    
-  }
-  
-  return false;
+  } return false;
 }
 
 bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcSphere* l, cgal_shape_t& shape) {
@@ -691,39 +614,11 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcSphere* l, cgal_shape_t& s
     } face_list = refined_face_list;
   }
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-//  std::ofstream fresult;
-//  fresult.open("/Users/ken/Desktop/sphere.off");
-//  fresult << polyhedron << std::endl;
-//  fresult.close();
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  } CGAL_postcondition(polyhedron.is_valid() && polyhedron.is_closed());
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  
   cgal_placement_t trsf;
   IfcGeom::CgalKernel::convert(l->Position(),trsf);
   
-  for (auto &vertex: vertices(polyhedron)) {
-    vertex->point() = Kernel::Point_3(vertex->point().x()*r,
-                                      vertex->point().y()*r,
-                                      vertex->point().z()*r).transform(trsf);
-  }
-  
-//  std::ofstream fresult;
-//  fresult.open("/Users/ken/Desktop/sphere.off");
-//  fresult << polyhedron << std::endl;
-//  fresult.close();
-  
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
+  shape.transform(trsf);
   return true;
 }
 
@@ -762,27 +657,11 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcRectangularPyramid* l, cga
   face_list.back().outer.push_back(Kernel::Point_3(0, 0, 0));
   face_list.back().outer.push_back(Kernel::Point_3(0.5*dx, 0.5*dy, dz));
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  } CGAL_postcondition(polyhedron.is_valid() && polyhedron.is_closed());
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  
   cgal_placement_t trsf;
   IfcGeom::CgalKernel::convert(l->Position(),trsf);
   
-  for (auto &vertex: vertices(polyhedron)) {
-    vertex->point() = vertex->point().transform(trsf);
-  }
-  
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
+  shape.transform(trsf);
   return true;
 }
 
@@ -820,27 +699,11 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcRightCircularCylinder* l,
     face_list.back().outer.push_back(Kernel::Point_3(r*cos(current_angle), r*sin(current_angle), h));
   }
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  } CGAL_postcondition(polyhedron.is_valid() && polyhedron.is_closed());
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  
   cgal_placement_t trsf;
   IfcGeom::CgalKernel::convert(l->Position(),trsf);
   
-  for (auto &vertex: vertices(polyhedron)) {
-    vertex->point() = vertex->point().transform(trsf);
-  }
-  
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
+  shape.transform(trsf);
   return true;
 }
 
@@ -870,27 +733,11 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcRightCircularCone* l, cgal
     face_list.back().outer.push_back(Kernel::Point_3(0, 0, h));
   }
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  } CGAL_postcondition(polyhedron.is_valid() && polyhedron.is_closed());
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  
   cgal_placement_t trsf;
   IfcGeom::CgalKernel::convert(l->Position(),trsf);
   
-  for (auto &vertex: vertices(polyhedron)) {
-    vertex->point() = vertex->point().transform(trsf);
-  }
-  
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
+  shape.transform(trsf);
   return true;
 }
 
@@ -940,20 +787,7 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcTriangulatedFaceSet* l, cg
     face_list.back().outer.push_back(c);
   }
   
-  // Naive creation
-  CGAL::Polyhedron_3<Kernel> polyhedron = CGAL::Polyhedron_3<Kernel>();
-  PolyhedronBuilder builder(&face_list);
-  polyhedron.delegate(builder);
-  
-  // Stitch edges
-  //  std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  CGAL::Polygon_mesh_processing::stitch_borders(polyhedron);
-  if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
-    CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
-  } CGAL_postcondition(polyhedron.is_valid() && polyhedron.is_closed());
-  //  std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl;
-  
-  shape = CGAL::Nef_polyhedron_3<Kernel>(polyhedron);
+  shape = create_nef_polyhedron(face_list);
   return true;
 }
 #endif
@@ -967,11 +801,13 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcHalfSpaceSolid* l, cgal_sh
   cgal_plane_t pln;
   IfcGeom::CgalKernel::convert((IfcSchema::IfcPlane*)surface,pln);
   
-  // TODO: This might be the other way around?
-  if (!l->AgreementFlag()) pln = pln.opposite();
+  // TODO: Don't fully understand the logic here. Might be incorrect.
+  if (l->AgreementFlag()) pln = pln.opposite();
 //  const gp_Pnt pnt = pln.Location().Translated( l->AgreementFlag() ? -pln.Axis().Direction() : pln.Axis().Direction());
 //  shape = BRepPrimAPI_MakeHalfSpace(BRepBuilderAPI_MakeFace(pln),pnt).Solid();
   
-  shape = CGAL::Nef_polyhedron_3<Kernel>(pln);
+  shape = CGAL::Nef_polyhedron_3<Kernel>();
+  // TODO: We return an empty Nef polyhedron for now and handle halfspace differences in IfcBooleanResult. The other option would be to switch to an extended kernel.
+//  shape = CGAL::Nef_polyhedron_3<Kernel>(pln);
   return true;
 }

src/ifcgeom/kernels/cgal/CgalIfcGeomWires.cpp

@@ -176,7 +176,6 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcCompositeCurve* l, cgal_wi
   return true;
 }
 
-// TODO: Project points to closest point in curve?
 bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcTrimmedCurve* l, cgal_wire_t& wire) {
   IfcSchema::IfcCurve* basis_curve = l->BasisCurve();
   bool isConic = basis_curve->is(IfcSchema::Type::IfcConic);
@@ -217,6 +216,7 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcTrimmedCurve* l, cgal_wire
   trim_cartesian &= has_pnts[0] && has_pnts[1];
   bool trim_cartesian_failed = !trim_cartesian;
   if ( trim_cartesian ) {
+    // TODO: Project points to closest point in curve?
     if ( CGAL::squared_distance(pnts[0], pnts[1]) < getValue(GV_WIRE_CREATION_TOLERANCE)*getValue(GV_WIRE_CREATION_TOLERANCE) ) {
       Logger::Message(Logger::LOG_WARNING,"Skipping segment with length below tolerance level:",l->entity);
       return false;
@@ -272,58 +272,24 @@ bool IfcGeom::CgalKernel::convert(const IfcSchema::IfcTrimmedCurve* l, cgal_wire
       const double magnitude = line->Dir()->Magnitude();
       flts[0] *= magnitude; flts[1] *= magnitude;
     }
-    if ( basis_curve->is(IfcSchema::Type::IfcEllipse) ) {
-      IfcSchema::IfcEllipse* ellipse = static_cast<IfcSchema::IfcEllipse*>(basis_curve);
-      double x = ellipse->SemiAxis1() * getValue(GV_LENGTH_UNIT);
-      double y = ellipse->SemiAxis2() * getValue(GV_LENGTH_UNIT);
-      const bool rotated = y > x;
-      if (rotated) {
-        flts[0] -= M_PI / 2.;
-        flts[1] -= M_PI / 2.;
-      }
-    }
     if ( isConic && ALMOST_THE_SAME(fmod(flts[1]-flts[0],M_PI*2.),0.) ) {
       for (auto &point: curve) w.push_back(point);
     } else {
-      if (l->SenseAgreement()) {
-        bool found = false;
-        int loops_to_go = 2;
-        std::vector<Kernel::Point_3>::const_iterator point = curve.begin();
-        do {
-          if (!found) {
-            if (CGAL::squared_distance(*point, pnts[0]) < getValue(GV_WIRE_CREATION_TOLERANCE)*getValue(GV_WIRE_CREATION_TOLERANCE)) {
-              found = true;
-              w.push_back(*point);
-            }
-          } else {
-            w.push_back(*point);
-            if (CGAL::squared_distance(*point, pnts[1]) < getValue(GV_WIRE_CREATION_TOLERANCE)*getValue(GV_WIRE_CREATION_TOLERANCE)) {
-              break;
-            }
-          } ++point;
-          if (point == curve.end()) {
-            point = curve.begin();
-            --loops_to_go;
-          }
-        } while (point != curve.begin() && loops_to_go > 0);
-      } else {
-        bool found = false;
-        int loops_to_go = 2;
-        std::vector<Kernel::Point_3>::const_reverse_iterator point = curve.rbegin();
-        do {
-          if (!found) {
-            if (CGAL::squared_distance(*point, pnts[0]) < getValue(GV_WIRE_CREATION_TOLERANCE)*getValue(GV_WIRE_CREATION_TOLERANCE)) {
-              found = true;
-              w.push_back(*point);
-            }
-          } else {
-            w.push_back(*point);
-            if (CGAL::squared_distance(*point, pnts[1]) < getValue(GV_WIRE_CREATION_TOLERANCE)*getValue(GV_WIRE_CREATION_TOLERANCE)) {
-              break;
-            }
-          } ++point;
-          if (point == curve.rend() && loops_to_go > 0) point = curve.rbegin();
-        } while (point != curve.rbegin());
+      const int segments_of_full_curve = 12;
+      double segment_angle = 2.0*3.141592653589793/segments_of_full_curve;
+      if ( basis_curve->is(IfcSchema::Type::IfcEllipse) ) {
+        IfcSchema::IfcEllipse* ellipse = static_cast<IfcSchema::IfcEllipse*>(basis_curve);
+        double x = ellipse->SemiAxis1() * getValue(GV_LENGTH_UNIT);
+        double y = ellipse->SemiAxis2() * getValue(GV_LENGTH_UNIT);
+        for (double current_angle = flts[0]; current_angle < flts[1]; current_angle += segment_angle) {
+          w.push_back(Kernel::Point_3(x*cos(current_angle), y*sin(current_angle), 0));
+        } w.push_back(Kernel::Point_3(x*cos(flts[1]), y*sin(flts[1]), 0));
+      } if ( basis_curve->is(IfcSchema::Type::IfcCircle) ) {
+        IfcSchema::IfcCircle* circle = static_cast<IfcSchema::IfcCircle*>(basis_curve);
+        double r = circle->Radius() * getValue(GV_LENGTH_UNIT);
+        for (double current_angle = flts[0]; current_angle < flts[1]; current_angle += segment_angle) {
+          w.push_back(Kernel::Point_3(r*cos(current_angle), r*sin(current_angle), 0));
+        } w.push_back(Kernel::Point_3(r*cos(flts[1]), r*sin(flts[1]), 0));
       }
     }
   } else if ( trim_cartesian_failed && (has_pnts[0] && has_pnts[1]) ) {

src/ifcgeom/kernels/cgal/CgalKernel.h

@@ -138,7 +138,9 @@ namespace IfcGeom {
     
     void remove_duplicate_points_from_loop(cgal_wire_t& polygon, bool closed, double tol = -1.);
 
-		 bool convert_openings(const IfcSchema::IfcProduct* entity, const IfcSchema::IfcRelVoidsElement::list::ptr& openings, const ConversionResults& entity_shapes, const cgal_placement_t& entity_trsf, ConversionResults& cut_shapes);
+    bool convert_openings(const IfcSchema::IfcProduct* entity, const IfcSchema::IfcRelVoidsElement::list::ptr& openings, const ConversionResults& entity_shapes, const cgal_placement_t& entity_trsf, ConversionResults& cut_shapes);
+    
+    CGAL::Nef_polyhedron_3<Kernel> create_nef_polyhedron(std::list<cgal_face_t> &face_list);
 
 		void purge_cache() {
 			// Rather hack-ish, but a stopgap solution to keep memory under control