Add autorefine for triangle soup

Installation/CHANGES.md

@@ -35,6 +35,12 @@ Release date: October 2023
 -   Removed the class templates `Gray_image_mesh_domain_3`, `Implicit_mesh_domain_3`, and `Labeled_image_mesh_domain_3`
     which are deprecated since CGAL-4.13.
 
+### [Polygon Mesh Processing](https://doc.cgal.org/6.0/Manual/packages.html#PkgPolygonMeshProcessing)
+
+-   Added the function
+    `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` that refines a soup of triangles so that no pair of triangles intersects
+     in their interiors. Also added, the function `autorefine()` operating directly on a triangle mesh and updating it
+     using the aforementioned function on a triangle soup.
 
 [Release 5.6](https://github.com/CGAL/cgal/releases/tag/v5.6)
 -----------

Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection.cpp

@@ -96,6 +96,13 @@ void test_coplanar_triangles(){
   assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
+  // TK10 case C'
+  t1=Triangle(Point(88.7921, 89.0007, 1.25), Point(88.1912, 88.3997, 1.25), Point(89.8224, 90.031, 1.25));
+  t2=Triangle(Point(88.0497, 88.2583, 1.25), Point(82.9292, 81.8747, 1.25), Point(91.1726, 91.3812, 1.25));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
 //Intersection is a point
   //edges  are collinear, one vertex in common
   t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
@@ -153,6 +160,13 @@ void test_coplanar_triangles(){
   assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
+  // TK10 case D
+  t1=Triangle(Point(-34.893700000000003, -16.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -18.5351, 3.1334899999999998e-12), Point(-42.393700000000003, -16.0351, 3.1334899999999998e-12));
+  t2=Triangle(Point(-34.893700000000003, -32.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -9.7851400000000002, 3.1334899999999998e-12), Point(-31.643699999999999, -17.201799999999999, 3.1334899999999998e-12));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
 //Intersection is a polygon
   //David's star
   t1=Triangle( Point(0,0,0),Point(1,0,0),Point(0.5,1.5,0) );
@@ -181,6 +195,51 @@ void test_coplanar_triangles(){
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
   assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case A
+  t1=Triangle(Point(3.74861, 12.4822, 14.0112), Point(5.40582, 12.4822, 15.6895), Point(5.37748, 12.4822, 15.7206));
+  t2=Triangle(Point(5.49972, 12.4822, 13.491), Point(5.27627, 12.4822, 15.8106), Point(5.32119, 12.4822, 15.8126));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case C
+  t1=Triangle(Point(5, -94.6659, 3.85175), Point(5, -94.5682, 3.08638), Point(5, -94.8182, 3.08638));
+  t2=Triangle(Point(5, -94.4317, 3.76399), Point(5, -97.6182, 3.08638), Point(5, -94.5659, 2.99682));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case E
+  t1=Triangle(Point(-955.858, -45.032, -0.016), Point(-955.856, -45.032, -0.004), Point(-955.856, -45.032, -0.002));
+  t2=Triangle(Point(-955.856, -45.032, 0.006), Point(-955.854, -45.032, -0.002), Point(-955.876, -45.032, -0.034));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case F
+  t1=Triangle(Point(141.172, 20.576, 155.764), Point(141.172, 20.588, 155.766), Point(141.172, 20.59, 155.766));
+  t2=Triangle(Point(141.172, 20.602, 155.768), Point(141.172, 20.594, 155.766), Point(141.172, 20.574, 155.764));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case D
+  t1=Triangle(Point(152.864, 126.324, 0.950001), Point(152.77, 126.483, 0.950001), Point(153.072, 125.973, 0.950001));
+  t2=Triangle(Point(153.322, 125.551, 0.950001), Point(152.218, 127.415, 0.950001), Point(153.66, 124.768, 0.950001));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
 //Intersection is empty
   t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
   t2=Triangle( Point(-0.1,-0.1,0),Point(-0.1,-0.9,0),Point(-1,-0.1,0) );

Kernel_23/include/CGAL/Kernel/function_objects.h

@@ -2192,6 +2192,107 @@ namespace CommonKernelFunctors {
     }
   };
 
+  template <typename K>
+  class Construct_planes_intersection_point_3
+  {
+    typedef typename K::Plane_3 Plane;
+    typedef typename K::Point_3 Point;
+    typename K::Construct_plane_3 construct_plane;
+  public:
+    typedef Point result_type;
+
+    Point
+    operator()(const Point& p1, const Point& q1, const Point& r1,
+               const Point& p2, const Point& q2, const Point& r2,
+               const Point& p3, const Point& q3, const Point& r3) const
+    {
+      Plane plane1 = construct_plane(p1, q1, r1);
+      Plane plane2 = construct_plane(p2, q2, r2);
+      Plane plane3 = construct_plane(p3, q3, r3);
+
+      const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+
+    Point
+    operator()(const Plane& plane1, const Plane& plane2, const Plane& plane3) const
+    {
+      const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+  };
+
+  template <typename K>
+  class Construct_coplanar_segments_intersection_point_3
+  {
+    typedef typename K::Segment_3 Segment;
+    typedef typename K::Point_3 Point;
+    typename K::Construct_segment_3 construct_segment;
+  public:
+    typedef Point result_type;
+
+    Point
+    operator()(const Point& p1, const Point& q1,
+               const Point& p2, const Point& q2) const
+    {
+      Segment s1 = construct_segment(p1, q1);
+      Segment s2 = construct_segment(p2, q2);
+
+      const auto res = typename K::Intersect_3()(s1, s2);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+
+    Point
+    operator()(const Segment& s1, const Segment& s2) const
+    {
+      const auto res = typename K::Intersect_3()(s1, s2);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+  };
+
+  template <typename K>
+  class Compute_alpha_for_coplanar_triangle_intersection_3
+  {
+    typedef typename K::Point_3 Point_3;
+    typedef typename K::Vector_3 Vector_3;
+  public:
+    typedef typename K::FT  result_type;
+    result_type
+    operator()(const Point_3& p1, const Point_3& p2,       // segment 1
+               const Point_3& p3, const Point_3& p4) const // segment 2
+    {
+      typename K::Construct_vector_3 vector = K().construct_vector_3_object();
+      typename K::Construct_cross_product_vector_3 cross_product =
+        K().construct_cross_product_vector_3_object();
+
+      const Vector_3 v1 = vector(p1, p2);
+      const Vector_3 v2 = vector(p3, p4);
+
+      CGAL_assertion(K().coplanar_3_object()(p1,p2,p3,p4));
+
+      const Vector_3 v3 = vector(p1, p3);
+      const Vector_3 v3v2 = cross_product(v3,v2);
+      const Vector_3 v1v2 = cross_product(v1,v2);
+      const typename K::FT sl = K().compute_squared_length_3_object()(v1v2);
+      CGAL_assertion(!certainly(is_zero(sl)));
+
+      const typename K::FT t = ((v3v2.x()*v1v2.x()) + (v3v2.y()*v1v2.y()) + (v3v2.z()*v1v2.z())) / sl;
+      return t; // p1 + (p2-p1) * t
+    }
+  };
+
   template <typename K>
   class Construct_point_on_2
   {

Kernel_23/include/CGAL/Kernel/interface_macros.h

@@ -398,6 +398,12 @@ CGAL_Kernel_cons(Construct_plane_3,
                  construct_plane_3_object)
 CGAL_Kernel_cons(Construct_plane_line_intersection_point_3,
                  construct_plane_line_intersection_point_3_object)
+CGAL_Kernel_cons(Construct_planes_intersection_point_3,
+                 construct_planes_intersection_point_3_object)
+CGAL_Kernel_cons(Construct_coplanar_segments_intersection_point_3,
+                 construct_coplanar_segments_intersection_point_3_object)
+CGAL_Kernel_cons(Compute_alpha_for_coplanar_triangle_intersection_3,
+                 compute_alpha_for_coplanar_triangle_intersection_3_object)
 CGAL_Kernel_cons(Construct_point_on_2,
                  construct_point_on_2_object)
 CGAL_Kernel_cons(Construct_point_on_3,

Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPAutorefinementVisitor.h

@@ -0,0 +1,25 @@
+/// \ingroup PkgPolygonMeshProcessingConcepts
+/// \cgalConcept
+///
+/// The concept `PMPAutorefinementVisitor` defines the requirements for the visitor
+/// used in `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` to track
+/// the creation of new triangles.
+///
+/// \cgalRefines{CopyConstructible}
+/// \cgalHasModel `CGAL::Polygon_mesh_processing::Autorefinement::Default_visitor`.
+
+class PMPAutorefinementVisitor{
+public:
+
+/// @name Functions called only if at least one intersection has been found
+/// @{
+  /// called when the final number of output triangles is known, `nbt` being the total number of triangles in the output.
+  void number_of_output_triangles(std::size_t nbt);
+  /// called for triangle with no intersection, `tgt_id` is the position in the triangle container after calling
+  /// `autorefine_triangle_soup()`, while `src_id` was its position before calling the function.
+  void verbatim_triangle_copy(std::size_t tgt_id, std::size_t src_id);
+  /// called for each subtriangle created from a triangle with intersection, `tgt_id` is the position in the triangle container after calling
+  /// `autorefine_triangle_soup()` of the subtriangle, while `src_id` was the position of the original support triangle before calling the function.
+  void new_subtriangle(std::size_t tgt_id, std::size_t src_id);
+/// @}
+};

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -109,6 +109,8 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - `CGAL::Polygon_mesh_processing::surface_intersection()`
 - `CGAL::Polygon_mesh_processing::clip()`
 - `CGAL::Polygon_mesh_processing::split()`
+- `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()`
+- `CGAL::Polygon_mesh_processing::autorefine()`
 
 \cgalCRPSection{Meshing Functions}
 - \link PMP_meshing_grp `CGAL::Polygon_mesh_processing::isotropic_remeshing()` \endlink

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -51,6 +51,8 @@ create_single_source_cgal_program("match_faces.cpp")
 create_single_source_cgal_program("cc_compatible_orientations.cpp")
 create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp")
 create_single_source_cgal_program("hausdorff_bounded_error_distance_example.cpp")
+create_single_source_cgal_program("triangle_mesh_autorefinement.cpp")
+create_single_source_cgal_program("soup_autorefinement.cpp")
 
 find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)
 include(CGAL_Eigen3_support)
@@ -121,6 +123,7 @@ if(TARGET CGAL::TBB_support)
   target_link_libraries(self_intersections_example PUBLIC CGAL::TBB_support)
   target_link_libraries(hausdorff_distance_remeshing_example PUBLIC CGAL::TBB_support)
   target_link_libraries(hausdorff_bounded_error_distance_example PUBLIC CGAL::TBB_support)
+  target_link_libraries(soup_autorefinement PUBLIC CGAL::TBB_support)
 
   create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp")
   target_link_libraries(corefinement_parallel_union_meshes PUBLIC CGAL::TBB_support)

Polygon_mesh_processing/examples/Polygon_mesh_processing/soup_autorefinement.cpp

@@ -0,0 +1,42 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/Real_timer.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3                                     Point;
+
+typedef CGAL::Surface_mesh<Point>                           Mesh;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
+                                         : std::string(argv[1]);
+
+  std::vector<Point> input_points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+  CGAL::IO::read_polygon_soup(filename, input_points, input_triangles);
+  PMP::repair_polygon_soup(input_points, input_triangles);
+  PMP::triangulate_polygons(input_points, input_triangles);
+
+  CGAL::Real_timer t;
+  t.start();
+  PMP::autorefine_triangle_soup(input_points, input_triangles,
+                                CGAL::parameters::concurrency_tag(CGAL::Parallel_if_available_tag()));
+  t.stop();
+  std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " in " << t.time() << " sec." << std::endl;
+  CGAL::IO::write_polygon_soup("autorefined.off", input_points, input_triangles, CGAL::parameters::stream_precision(17));
+
+  return 0;
+}

Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp

@@ -0,0 +1,29 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3                                     Point;
+
+typedef CGAL::Surface_mesh<Point>                           Mesh;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  Mesh mesh;
+
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
+                                         : std::string(argv[1]);
+  CGAL::IO::read_polygon_mesh(filename, mesh);
+
+  PMP::autorefine(mesh);
+
+  CGAL::IO::write_polygon_mesh("autorefined.off", mesh, CGAL::parameters::stream_precision(17));
+
+  return 0;
+}

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h

@@ -28,8 +28,8 @@ struct Non_manifold_feature_map
   typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
   typedef dynamic_edge_property_t<std::size_t> Edge_to_id_tag;
   typedef dynamic_vertex_property_t<std::size_t> Vertex_to_id_tag;
-  typedef typename boost::property_map<PolygonMesh, Edge_to_id_tag>::type Edge_to_nm_id;
-  typedef typename boost::property_map<PolygonMesh, Vertex_to_id_tag>::type Vertex_to_nm_id;
+  typedef typename boost::property_map<PolygonMesh, Edge_to_id_tag>::const_type Edge_to_nm_id;
+  typedef typename boost::property_map<PolygonMesh, Vertex_to_id_tag>::const_type Vertex_to_nm_id;
   Edge_to_nm_id e_nm_id;
   Vertex_to_nm_id v_nm_id;
   std::vector< std::vector<edge_descriptor> > non_manifold_edges;
@@ -39,7 +39,7 @@ struct Non_manifold_feature_map
   {}
 
   template <class Vpm>
-  Non_manifold_feature_map(PolygonMesh& pm, Vpm vpm)
+  Non_manifold_feature_map(const PolygonMesh& pm, Vpm vpm)
     : e_nm_id(get(Edge_to_id_tag(), pm))
     , v_nm_id(get(Vertex_to_id_tag(), pm))
   {
@@ -99,6 +99,14 @@ struct Non_manifold_feature_map
       }
     }
   }
+
+  void clear()
+  {
+    non_manifold_edges.clear();
+    non_manifold_vertices.clear();
+    e_nm_id = Edge_to_nm_id();
+    v_nm_id = Vertex_to_nm_id();
+  }
 };
 
 } } // end of CGAL::Polygon_mesh_processing