Add experimental functions to compute the kernel of a triangle mesh

Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h

@@ -28,49 +28,59 @@ namespace internal {
 //  triple plane intersection
 template <class K>
 std::optional<typename K::Point_3>
-intersection_point(const typename K::Plane_3& plane1,
-                   const typename K::Plane_3& plane2,
-                   const typename K::Plane_3& plane3,
-                   const K&)
+intersection_point_RT(const typename K::Plane_3& plane1,
+                      const typename K::Plane_3& plane2,
+                      const typename K::Plane_3& plane3,
+                      const K&)
 {
-  typedef typename K::FT FT;
-
-  const FT& m00 = plane1.a();
-  const FT& m01 = plane1.b();
-  const FT& m02 = plane1.c();
-  const FT& b0  = - plane1.d();
-  const FT& m10 = plane2.a();
-  const FT& m11 = plane2.b();
-  const FT& m12 = plane2.c();
-  const FT& b1  = - plane2.d();
-  const FT& m20 = plane3.a();
-  const FT& m21 = plane3.b();
-  const FT& m22 = plane3.c();
-  const FT& b2  = - plane3.d();
+  typedef typename K::RT RT;
+
+  const RT& m00 = plane1.a();
+  const RT& m01 = plane1.b();
+  const RT& m02 = plane1.c();
+  const RT& b0  = - plane1.d();
+  const RT& m10 = plane2.a();
+  const RT& m11 = plane2.b();
+  const RT& m12 = plane2.c();
+  const RT& b1  = - plane2.d();
+  const RT& m20 = plane3.a();
+  const RT& m21 = plane3.b();
+  const RT& m22 = plane3.c();
+  const RT& b2  = - plane3.d();
 
   // Minors common to two determinants
-  const FT minor_0 = m00*m11 - m10*m01;
-  const FT minor_1 = m00*m21 - m20*m01;
-  const FT minor_2 = m10*m21 - m20*m11;
+  const RT minor_0 = m00*m11 - m10*m01;
+  const RT minor_1 = m00*m21 - m20*m01;
+  const RT minor_2 = m10*m21 - m20*m11;
 
-  const FT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M
+  const RT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M
 
   if(is_zero(den)){
     return std::nullopt;
   }
 
-  const FT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0;  // determinant of M with M[x:2] swapped with [b0,b1,b2]
+  const RT num3 = minor_0*b2 - minor_1*b1 + minor_2*b0;  // determinant of M with M[x:2] swapped with [b0,b1,b2]
 
   // Minors common to two determinants
-  const FT minor_3 = b0*m12 - b1*m02;
-  const FT minor_4 = b0*m22 - b2*m02;
-  const FT minor_5 = b1*m22 - b2*m12;
+  const RT minor_3 = b0*m12 - b1*m02;
+  const RT minor_4 = b0*m22 - b2*m02;
+  const RT minor_5 = b1*m22 - b2*m12;
 
   // num1 has opposite signs because b0 and M[:1] have been swapped
-  const FT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01;  // determinant of M with M[x:0] swapped with [b0,b1,b2]
-  const FT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00;  // determinant of M with M[x:1] swapped with [b0,b1,b2]
+  const RT num1 = - minor_3*m21 + minor_4*m11 - minor_5*m01;  // determinant of M with M[x:0] swapped with [b0,b1,b2]
+  const RT num2 = minor_3*m20 - minor_4*m10 + minor_5*m00;  // determinant of M with M[x:1] swapped with [b0,b1,b2]
 
-  return std::make_optional(typename K::Point_3(num1/den, num2/den, num3/den));
+  return typename K::Point_3(num1, num2, num3, den);
+}
+
+template <class K>
+std::optional<typename K::Point_3>
+intersection_point(const typename K::Plane_3& plane1,
+                   const typename K::Plane_3& plane2,
+                   const typename K::Plane_3& plane3,
+                   const K& k)
+{
+  return intersection_point_RT(plane1, plane2, plane3, k);
 }
 
 template <class K>

Kernel_23/include/CGAL/Kernel/function_objects.h

@@ -2097,25 +2097,38 @@ namespace CommonKernelFunctors {
   template <typename K>
   class Construct_plane_line_intersection_point_3
   {
-    typedef typename K::Plane_3 Plane;
+    typedef typename K::FT FT;
     typedef typename K::Line_3 Line;
     typedef typename K::Point_3 Point;
-    typename K::Construct_plane_3 construct_plane;
-    typename K::Construct_line_3 construct_line;
+    typedef typename K::Plane_3 Plane;
+    typedef typename K::Vector_3 Vector;
 
+    typename K::Construct_cross_product_vector_3 construct_cross_product_vector_3;
+    typename K::Compute_scalar_product_3 scalar_product;
+
+    typename K::Construct_line_3 construct_line;
+    typename K::Construct_vector_3 vector;
+    typename K::Construct_plane_3 construct_plane;
   public:
     Point
     operator()(const Point& p1, const Point& p2, const Point& p3,
                const Point& l1, const Point& l2) const
     {
-      Plane plane = construct_plane(p1, p2, p3);
+      Plane plane = construct_plane( p1, p2, p3 );
       Line line = construct_line( l1, l2 );
 
       const auto res = typename K::Intersect_3()(plane,line);
       CGAL_assertion(res!=std::nullopt);
       const Point* e_pt = std::get_if<Point>(&(*res));
       CGAL_assertion(e_pt!=nullptr);
       return *e_pt;
+
+      // Vector n = construct_cross_product_vector_3( vector(p1, p2), vector(p1, p3));
+      // CGAL_assertion(n!=CGAL::NULL_VECTOR);
+      // Vector l = vector(l1, l2);
+      // CGAL_assertion(l!=CGAL::NULL_VECTOR);
+      // FT t = scalar_product(n, vector(p1,l1)) / scalar_product(n, l);
+      // return l1 + t * l;
     }
 
     Point

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h

@@ -265,6 +265,24 @@ bool close(PolygonMesh& pm, VertexPointMap vpm, typename Traits::Vector_3 plane_
 
 #endif
 
+template <class PolygonMesh,
+          class Visitor>
+void naive_close(PolygonMesh& pm, Visitor& visitor)
+{
+  using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
+
+  std::vector< halfedge_descriptor > border_cycles;
+  extract_boundary_cycles(pm, std::back_inserter(border_cycles));
+  std::vector< Bbox_3 > bboxes;
+
+  for (halfedge_descriptor h : border_cycles)
+  {
+    visitor.before_face_copy(boost::graph_traits<PolygonMesh>::null_face(), pm, pm);
+    Euler::fill_hole(h, pm);
+    visitor.after_face_copy(boost::graph_traits<PolygonMesh>::null_face(), pm, face(h, pm), pm);
+  }
+}
+
 
 template <class Plane_3,
           class TriangleMesh,
@@ -994,7 +1012,7 @@ bool clip(PolygonMesh& pm,
 {
   using parameters::choose_parameter;
   using parameters::get_parameter;
-  using parameters::get_parameter_reference     ;
+  using parameters::get_parameter_reference;
 
   using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
 
@@ -1026,11 +1044,12 @@ bool clip(PolygonMesh& pm,
   const bool allow_self_intersections =
     choose_parameter(get_parameter(np, internal_np::allow_self_intersections), false);
   bool triangulate = !choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), false);
-
+  constexpr bool traits_supports_cdt2 = !internal::Has_member_Does_not_support_CDT2<GT>::value;
+  const bool used_for_kernel = choose_parameter(get_parameter(np, internal_np::used_for_kernel), false);
   auto vos = get(dynamic_vertex_property_t<Oriented_side>(), pm);
   auto ecm = get(dynamic_edge_property_t<bool>(), pm, false);
 
-  if (triangulate && !is_triangle_mesh(pm))
+  if (traits_supports_cdt2 && triangulate && !is_triangle_mesh(pm))
     triangulate = false;
 
   refine_with_plane(pm, plane, parameters::vertex_oriented_side_map(vos)
@@ -1082,15 +1101,26 @@ bool clip(PolygonMesh& pm,
 
   remove_connected_components(pm, ccs_to_remove, fcc);
 
-  if (clip_volume)
+  if constexpr (traits_supports_cdt2)
   {
-    //TODO: add in the traits construct_orthogonal_vector
-    if (triangulate)
-      internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
-    else
-      if (!internal::close<GT>(pm, vpm, plane.orthogonal_vector(), visitor))
-        internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
+    if (clip_volume)
+    {
+      if (!used_for_kernel)
+      {
+        //TODO: add in the traits construct_orthogonal_vector
+        if (triangulate)
+          internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
+        else
+          if (!internal::close<GT>(pm, vpm, plane.orthogonal_vector(), visitor))
+            internal::close_and_triangulate<GT>(pm, vpm, plane.orthogonal_vector(), visitor);
+      }
+      else
+        internal::naive_close(pm, visitor);
+    }
   }
+  else
+    if (clip_volume)
+      internal::naive_close(pm, visitor);
 
   return true;
 }