Surface_mesh_approximation: Deal with boundary edges (#7574)

## Summary of Changes

Add a named parameter to distinguish between the approximation error for
boundary and non-boundary edges of the input mesh.

## Release Management

* Affected package(s): Surface_mesh_approximation
* Feature/Small Feature (if any): 
* Link to compiled documentation (obligatory for small feature) [*wrong
link name to be changed*](httpssss://wrong_URL_to_be_changed/Manual/Pkg)
* License and copyright ownership:  unchanged

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -229,10 +229,12 @@ CGAL_add_named_parameter(number_of_relaxations_t, number_of_relaxations, number_
 CGAL_add_named_parameter(use_convex_hull_t, use_convex_hull, use_convex_hull)
 
 // meshing parameters
+CGAL_add_named_parameter(boundary_subdivision_ratio_t, boundary_subdivision_ratio, boundary_subdivision_ratio)
 CGAL_add_named_parameter(subdivision_ratio_t, subdivision_ratio, subdivision_ratio)
 CGAL_add_named_parameter(relative_to_chord_t, relative_to_chord, relative_to_chord)
 CGAL_add_named_parameter(with_dihedral_angle_t, with_dihedral_angle, with_dihedral_angle)
 CGAL_add_named_parameter(optimize_anchor_location_t, optimize_anchor_location, optimize_anchor_location)
+CGAL_add_named_parameter(optimize_boundary_anchor_location_t, optimize_boundary_anchor_location, optimize_boundary_anchor_location)
 CGAL_add_named_parameter(pca_plane_t, pca_plane, pca_plane)
 
 // tetrahedral remeshing parameters

Surface_mesh_approximation/doc/Surface_mesh_approximation/Surface_mesh_approximation.txt

@@ -121,7 +121,7 @@ In order to approximate complex boundaries well, more anchors are generated by r
 \f[ d = d / input\_mesh\_average\_edge\_length. \f]
 Optionally, \f$ d \f$ can be measured as the ratio of the chord length:
 \f[ d = d / \Vert(\mathbf{a}, \mathbf{b})\Vert. \f]
-Also, we can add a dihedral angle weight \f$ sin(\mathbf{N}_i,\mathbf{N}_j) \f$ to the distance measurement, where \f$ \mathbf{N}_i,\mathbf{N}_j \f$ are the normals of the proxies separated by the chord \f$ (\mathbf{a}, \mathbf{b}) \f$. If the angle between proxy \f$ P_i \f$ and \f$ P_j \f$ is rather small, then a coarse approximation will do as it does not add geometric information on the shape. Trivial chords (less than 4 edges) are not subdivided if they are non-circular. In case of circular chords, additional anchors maybe added to maintain the topology even if they are trivial, as detailed in Section \ref sma_anchors_additional.
+Also, we can add a dihedral angle weight \f$ sin(\mathbf{N}_i,\mathbf{N}_j) \f$ to the distance measurement, where \f$ \mathbf{N}_i,\mathbf{N}_j \f$ are the normals of the proxies separated by the chord \f$ (\mathbf{a}, \mathbf{b}) \f$. If the angle between proxy \f$ P_i \f$ and \f$ P_j \f$ is rather small, then a coarse approximation will do as it does not add geometric information on the shape. Trivial chords (made of a single edge) are not subdivided. In case of circular chords, additional anchors may be added to maintain the topology, as detailed in Section \ref sma_anchors_additional.
 
 \cgalFigureBegin{chord, chord.jpg}
 Varying the chord error. From left to right: clustering partition, and meshing with decreasing absolute chord error 5, 3 and 1 without dihedral angle weight. The boundaries of the partition (red lines) are approximated with increasing accuracy.
@@ -132,7 +132,7 @@ Varying the chord error. From left to right: clustering partition, and meshing w
 For a boundary cycle without any anchor such as the hole depicted Figure \cgalFigureRef{operations}, we first add a starting anchor to the boundary. We then subdivide this circular chord to ensure that every boundary cycle has at least 2 anchors (i.e., every chord is connecting 2 different anchors, Figure \cgalFigureRef{anchors}). Finally, we add additional anchors to ensure that at least three anchor vertices are generated on every boundary cycle.
 
 \cgalFigureBegin{anchors, anchors.jpg}
-Adding anchors. From left to right: starting from a partition (grey) with a hole (while) and two encircled regions (green and blue), we add a starting anchor (orange disk) to the boundary cycle (red dash line) without any anchor (2nd), subdivide the circular chord (3rd, the number indicates the level of recursion) and add anchors to the boundary cycle with less than 2 anchors (4th, red dash lines).
+Adding anchors. From left to right: starting from a partition (grey) with a hole (white) and two encircled regions (green and blue), we add a starting anchor (orange disk) to the boundary cycle (red dash line) without any anchor (2nd), subdivide the circular chord (3rd, the number indicates the level of recursion) and add anchors to the boundary cycle with less than 2 anchors (4th, red dash lines).
 \cgalFigureEnd
 
 \subsubsection sma_triangulation Discrete Triangulation

Surface_mesh_approximation/include/CGAL/Variational_shape_approximation.h

@@ -123,6 +123,7 @@ class Variational_shape_approximation {
   typedef typename Geom_traits::FT FT;
   typedef typename Geom_traits::Point_3 Point_3;
   typedef typename Geom_traits::Vector_3 Vector_3;
+  typedef typename Geom_traits::Segment_3 Segment_3;
   typedef typename Geom_traits::Plane_3 Plane_3;
   typedef typename Geom_traits::Construct_vector_3 Construct_vector_3;
   typedef typename Geom_traits::Construct_point_3 Construct_point_3;
@@ -795,6 +796,12 @@ class Variational_shape_approximation {
    *     \cgalParamDefault{`5.0`}
    *   \cgalParamNEnd
    *
+   *   \cgalParamNBegin{boundary_subdivision_ratio}
+   *     \cgalParamDescription{the chord subdivision ratio threshold to the chord length or average edge length for boundary edges}
+   *     \cgalParamType{`geom_traits::FT`}
+   *     \cgalParamDefault{`subdivision_ratio`}
+   *   \cgalParamNEnd
+   *
    *   \cgalParamNBegin{relative_to_chord}
    *     \cgalParamDescription{If `true`, the `subdivision_ratio` is the ratio of the furthest vertex distance
    *                           to the chord length, otherwise is the average edge length}
@@ -814,6 +821,12 @@ class Variational_shape_approximation {
    *     \cgalParamDefault{`true`}
    *   \cgalParamNEnd
    *
+   *   \cgalParamNBegin{optimize_boundary_anchor_location}
+   *     \cgalParamDescription{If `true`, optimize the anchor locations of boundary vertices}
+   *     \cgalParamType{`Boolean`}
+   *     \cgalParamDefault{`optimize_anchor_location`}
+   *   \cgalParamNEnd
+   *
    *   \cgalParamNBegin{pca_plane}
    *     \cgalParamDescription{If `true`, use PCA plane fitting, otherwise use the default area averaged plane parameters}
    *     \cgalParamType{`Boolean`}
@@ -827,9 +840,11 @@ class Variational_shape_approximation {
     using parameters::choose_parameter;
 
     const FT subdivision_ratio = choose_parameter(get_parameter(np, internal_np::subdivision_ratio), FT(5.0));
+    const FT boundary_subdivision_ratio = choose_parameter(get_parameter(np, internal_np::boundary_subdivision_ratio), subdivision_ratio);
     const bool relative_to_chord = choose_parameter(get_parameter(np, internal_np::relative_to_chord), false);
     const bool with_dihedral_angle = choose_parameter(get_parameter(np, internal_np::with_dihedral_angle), false);
     const bool optimize_anchor_location = choose_parameter(get_parameter(np, internal_np::optimize_anchor_location), true);
+    const bool optimize_boundary_anchor_location = choose_parameter(get_parameter(np, internal_np::optimize_boundary_anchor_location), optimize_anchor_location);
     const bool pca_plane = choose_parameter(get_parameter(np, internal_np::pca_plane), false);
 
     // compute averaged edge length, used in chord subdivision
@@ -848,14 +863,14 @@ class Variational_shape_approximation {
 
     // generate anchors
     find_anchors();
-    find_edges(subdivision_ratio, relative_to_chord, with_dihedral_angle);
+    find_edges(subdivision_ratio, boundary_subdivision_ratio, relative_to_chord, with_dihedral_angle);
     add_anchors();
 
     // discrete constrained Delaunay triangulation
     pseudo_cdt();
 
     if (optimize_anchor_location)
-      this->optimize_anchor_location();
+      this->optimize_anchor_location(optimize_boundary_anchor_location);
 
     // check manifold-oriented
     return Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(m_tris);
@@ -1346,7 +1361,7 @@ class Variational_shape_approximation {
     const Proxy px = m_metric->fit_proxy(fvec, *m_ptm);
     const FT err = m_metric->compute_error(f, *m_ptm, px);
 
-    // original proxy map should always be falid
+    // original proxy map should always be valid
     const std::size_t prev_px_idx = get(m_fproxy_map, f);
     CGAL_assertion(prev_px_idx != CGAL_VSA_INVALID_TAG);
     // update the proxy error and proxy map
@@ -1513,12 +1528,13 @@ class Variational_shape_approximation {
 
   /*!
    * @brief finds and approximates the chord connecting the anchors.
-   * @param subdivision_ratio boundary chord approximation recursive split creterion
+   * @param subdivision_ratio boundary chord approximation recursive split criterion
    * @param relative_to_chord set `true` if the subdivision_ratio is relative to the chord length (relative sense),
    * otherwise it's relative to the average edge length (absolute sense).
    * @param with_dihedral_angle if set to `true`, add dihedral angle weight to the distance.
    */
   void find_edges(const FT subdivision_ratio,
+    const FT boundary_subdivision_ratio,
     const bool relative_to_chord,
     const bool with_dihedral_angle) {
     // collect candidate halfedges in a set
@@ -1550,7 +1566,7 @@ class Variational_shape_approximation {
         Boundary_chord chord;
         walk_to_next_anchor(he_start, chord);
         m_bcycles.back().num_anchors += subdivide_chord(chord.begin(), chord.end(),
-          subdivision_ratio, relative_to_chord, with_dihedral_angle);
+            subdivision_ratio, boundary_subdivision_ratio, relative_to_chord, with_dihedral_angle);
 
 #ifdef CGAL_SURFACE_MESH_APPROXIMATION_DEBUG
         std::cerr << "#chord_anchor " << m_bcycles.back().num_anchors << std::endl;
@@ -1837,15 +1853,16 @@ class Variational_shape_approximation {
    * @param chord_begin begin iterator of the chord
    * @param chord_end end iterator of the chord
    * @param subdivision_ratio the chord recursive split error threshold
-   * @param relative_to_chord set `true` if the subdivision_ratio is relative to the chord length (relative sense),
-   * otherwise it's relative to the average edge length (absolute sense).
+   * @param relative_to_chord set `true` if the `subdivision_ratio` is relative to the chord length (relative sense),
+   * otherwise it is relative to the average edge length (absolute sense).
    * @param with_dihedral_angle if set to `true` add dihedral angle weight to the distance.
    * @return the number of anchors of the chord apart from the first one
    */
   std::size_t subdivide_chord(
     const Boundary_chord_iterator &chord_begin,
     const Boundary_chord_iterator &chord_end,
     const FT subdivision_ratio,
+    const FT boundary_subdivision_ratio,
     const bool relative_to_chord,
     const bool with_dihedral_angle) {
     const std::size_t chord_size = std::distance(chord_begin, chord_end);
@@ -1854,8 +1871,16 @@ class Variational_shape_approximation {
     const std::size_t anchor_first = get(m_vanchor_map, source(he_first, *m_ptm));
     const std::size_t anchor_last = get(m_vanchor_map, target(he_last, *m_ptm));
 
+    bool is_boundary = is_border_edge(he_first, *m_ptm);
+
+    if(is_boundary && boundary_subdivision_ratio == 0){
+      for (Boundary_chord_iterator citr = chord_begin; *citr != he_last; ++citr) {
+        attach_anchor(*citr);
+      }
+    }
+
     // do not subdivide trivial non-circular chord
-    if ((anchor_first != anchor_last) && (chord_size < 4))
+    if ((anchor_first != anchor_last) && (chord_size < 2))
       return 1;
 
     bool if_subdivide = false;
@@ -1884,11 +1909,9 @@ class Variational_shape_approximation {
       const FT chord_len = CGAL::approximate_sqrt(chord_vec.squared_length());
       bool degenerate_chord = false;
       if (chord_len > FT(0.0)) {
-        chord_vec = scale_functor(chord_vec, FT(1.0) / chord_len);
+        Segment_3 seg(pt_begin, pt_end);
         for (Boundary_chord_iterator citr = chord_begin; citr != chord_end; ++citr) {
-          Vector_3 vec = vector_functor(pt_begin, m_vpoint_map[target(*citr, *m_ptm)]);
-          vec = cross_product_functor(chord_vec, vec);
-          const FT dist = CGAL::approximate_sqrt(vec.squared_length());
+          const FT dist = CGAL::approximate_sqrt(CGAL::squared_distance(m_vpoint_map[target(*citr, *m_ptm)], seg));
           if (dist > dist_max) {
             chord_max = citr;
             dist_max = dist;
@@ -1927,19 +1950,23 @@ class Variational_shape_approximation {
         }
         criterion *= norm_sin;
       }
-
-      if (criterion > subdivision_ratio)
-        if_subdivide = true;
+      if (is_boundary) {
+        if (criterion > boundary_subdivision_ratio)
+          if_subdivide = true;
+      }
+      else {
+        if (criterion > subdivision_ratio)
+          if_subdivide = true;
+      }
     }
-
     if (if_subdivide) {
       // subdivide at the most remote vertex
       attach_anchor(*chord_max);
 
       const std::size_t num_left = subdivide_chord(chord_begin, chord_max + 1,
-        subdivision_ratio, relative_to_chord, with_dihedral_angle);
+        subdivision_ratio, boundary_subdivision_ratio, relative_to_chord, with_dihedral_angle);
       const std::size_t num_right = subdivide_chord(chord_max + 1, chord_end,
-        subdivision_ratio, relative_to_chord, with_dihedral_angle);
+        subdivision_ratio, boundary_subdivision_ratio, relative_to_chord, with_dihedral_angle);
 
       return num_left + num_right;
     }
@@ -1992,9 +2019,15 @@ class Variational_shape_approximation {
    * @brief optimizes the anchor location by averaging the projection points of
    * the anchor vertex to the incident proxy plane.
    */
-  void optimize_anchor_location() {
+  void optimize_anchor_location(bool optimize_boundary_anchor_location) {
     for(Anchor& a : m_anchors) {
       const vertex_descriptor v = a.vtx;
+
+      if(! optimize_boundary_anchor_location && is_border(v,*m_ptm)){
+        a.pos  = m_vpoint_map[v];
+        continue;
+      }
+
       // incident proxy set
       std::set<std::size_t> px_set;
       for(halfedge_descriptor h : halfedges_around_target(v, *m_ptm)) {
@@ -2003,6 +2036,7 @@ class Variational_shape_approximation {
       }
 
       // projection
+      // todo: replace averaging by qem/svd ?  Mael?
       FT sum_area(0.0);
       Vector_3 vec = CGAL::NULL_VECTOR;
       const Point_3 vtx_pt = m_vpoint_map[v];

Surface_mesh_approximation/test/Surface_mesh_approximation/CMakeLists.txt

@@ -15,6 +15,9 @@ if(NOT TARGET CGAL::Eigen3_support)
   return()
 endif()
 
+create_single_source_cgal_program("vsa_border_test.cpp")
+target_link_libraries(vsa_border_test PUBLIC CGAL::Eigen3_support)
+
 create_single_source_cgal_program("vsa_class_interface_test.cpp")
 target_link_libraries(vsa_class_interface_test PUBLIC CGAL::Eigen3_support)