Enable to fix the coordinates of boundary vertices

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

@@ -227,6 +227,7 @@ CGAL_add_named_parameter(subdivision_ratio_t, subdivision_ratio, subdivision_rat
 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/include/CGAL/Variational_shape_approximation.h

@@ -123,6 +123,7 @@ private:
   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;
@@ -820,6 +821,12 @@ public:
    *     \cgalParamDefault{`true`}
    *   \cgalParamNEnd
    *
+   *   \cgalParamNBegin{optimize_boundary_anchor_location}
+   *     \cgalParamDescription{If `true`, optimize the anchor locations of boundary vertices}
+   *     \cgalParamType{`Boolean`}
+   *     \cgalParamDefault{`true`}
+   *   \cgalParamNEnd
+   *
    *   \cgalParamNBegin{pca_plane}
    *     \cgalParamDescription{If `true`, use PCA plane fitting, otherwise use the default area averaged plane parameters}
    *     \cgalParamType{`Boolean`}
@@ -837,6 +844,7 @@ public:
     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), true);
     const bool pca_plane = choose_parameter(get_parameter(np, internal_np::pca_plane), false);
 
     // compute averaged edge length, used in chord subdivision
@@ -862,7 +870,7 @@ public:
     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);
@@ -1353,7 +1361,7 @@ private:
     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
@@ -1520,7 +1528,7 @@ private:
 
   /*!
    * @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.
@@ -1845,8 +1853,8 @@ private:
    * @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),
+   * @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).
+   * 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
    */
@@ -1865,8 +1873,14 @@ private:
 
     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;
@@ -1895,11 +1909,13 @@ private:
         const FT chord_len = CGAL::approximate_sqrt(chord_vec.squared_length());
         bool degenerate_chord = false;
         if (chord_len > FT(0.0)) {
+          Segment_3 seg(pt_begin, pt_end);
             chord_vec = scale_functor(chord_vec, FT(1.0) / chord_len);
             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)]);
+              //Vector_3 vec = vector_functor(pt_begin, m_vpoint_map[target(*citr, *m_ptm)]);
-                vec = cross_product_functor(chord_vec, vec);
+              //vec = cross_product_functor(chord_vec, vec);
-                const FT dist = CGAL::approximate_sqrt(vec.squared_length());
+              //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;
@@ -2007,9 +2023,15 @@ private:
    * @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)) {
@@ -2018,6 +2040,7 @@ private:
       }
 
       // 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)
 

Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_border_test.cpp

@@ -0,0 +1,131 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Variational_shape_approximation.h>
+#include <boost/functional/hash_fwd.hpp>
+#include <boost/property_map/vector_property_map.hpp>
+
+#include <CGAL/Surface_mesh_approximation/L21_metric_plane_proxy.h>
+#include <CGAL/Variational_shape_approximation.h>
+
+#include <limits>
+
+#include <CGAL/Surface_mesh.h>
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using SurfaceMesh = CGAL::Surface_mesh<K::Point_3>;
+
+
+typedef boost::property_map<SurfaceMesh, boost::vertex_point_t>::type VertexPointMap;
+typedef typename CGAL::Kernel_traits<
+    typename boost::property_traits<VertexPointMap>::value_type>::Kernel GeomTraits;
+typedef CGAL::Surface_mesh_approximation::
+    L21_metric_plane_proxy<SurfaceMesh, VertexPointMap, GeomTraits>
+        L21_Metric;
+typedef CGAL::Variational_shape_approximation<SurfaceMesh, VertexPointMap, L21_Metric>
+    L21_MeshApproximation;
+
+struct PolygonSoup
+{
+    using PointType = K::Point_3;
+    std::vector<PointType> points;
+    std::vector<std::array<std::size_t, 3>> faces;
+};
+
+auto meshArea = [](const SurfaceMesh& mesh) {
+    double area = 0;
+    for (const auto& f : mesh.faces())
+    {
+        std::vector<K::Point_3> pts;
+        for (const auto& v : mesh.vertices_around_face(mesh.halfedge(f)))
+        {
+            pts.push_back(mesh.point(v));
+        }
+        K::Triangle_3 tri(pts[0], pts[1], pts[2]);
+        area += std::sqrt(tri.squared_area());
+    }
+    return area;
+};
+
+auto soupArea = [](const PolygonSoup& soup) {
+    double area = 0;
+    for (const auto& f : soup.faces)
+    {
+        std::vector<K::Point_3> pts = {soup.points[f[0]], soup.points[f[1]], soup.points[f[2]]};
+        K::Triangle_3 tri(pts[0], pts[1], pts[2]);
+        area += std::sqrt(tri.squared_area());
+    }
+    return area;
+};
+auto soupToMesh = [](const PolygonSoup& soup) {
+    SurfaceMesh mesh;
+    std::vector<SurfaceMesh::Vertex_index> vertices;
+    for (const auto& p : soup.points)
+    {
+        vertices.push_back(mesh.add_vertex(p));
+    }
+    for (const auto& f : soup.faces)
+    {
+        mesh.add_face(vertices[f[0]], vertices[f[1]], vertices[f[2]]);
+    }
+    return mesh;
+};
+
+int main(int, char*[])
+{
+    std::array<double,6> subdivisionRatios          = { 0.10 , 2.0 , 2.0 , 10.0 , 10.0 , 10.00 };
+    std::array<double, 6> boundarySubdivisionRatios = { 0.01 , 2.0 , 0.1 , 10.0 ,  1.0 ,  0.01 };
+
+    for (int i = 0; i < subdivisionRatios.size(); ++i)
+    {
+        const auto subdivisionRatio = subdivisionRatios[i];
+        const auto boundarySubdivisionRatio = boundarySubdivisionRatios[i];
+        std::string filename = "./VSA-" + std::to_string(subdivisionRatio) + "-" + std::to_string(boundarySubdivisionRatio);
+        SurfaceMesh mesh;
+
+        std::ifstream in("data/patch.ply");
+        CGAL::IO::read_PLY(in, mesh);
+        std::cout << vertices(mesh).size() << std::endl;
+        const int maxProxies = 10;
+        const int defaultIterations = 20;
+
+        VertexPointMap vpmap = get(boost::vertex_point, mesh);
+        L21_Metric metric(mesh, vpmap);
+        L21_MeshApproximation approx(mesh, vpmap, metric);
+
+        const auto seedingParams = //
+            CGAL::parameters::seeding_method(
+                CGAL::Surface_mesh_approximation::Seeding_method::INCREMENTAL)
+                .max_number_of_proxies(maxProxies)
+                .min_error_drop(0.001)
+                .number_of_relaxations(5);
+        approx.initialize_seeds(seedingParams);
+
+        approx.run(defaultIterations);
+
+        const auto meshParams = CGAL::parameters::subdivision_ratio(subdivisionRatio)
+                                    .boundary_subdivision_ratio(boundarySubdivisionRatio)
+                                    .with_dihedral_angle(false)
+                                    .optimize_boundary_anchor_location(false)
+                                    .optimize_anchor_location(true)
+                                    .pca_plane(false);
+        const bool isOutputManifold = approx.extract_mesh(meshParams);
+        if (!isOutputManifold)
+        {
+            throw std::runtime_error("not manifold");
+        }
+
+        PolygonSoup soup;
+        const auto outputParams = CGAL::parameters::anchors(std::back_inserter(soup.points))
+                                      .triangles(std::back_inserter(soup.faces));
+        approx.output(outputParams);
+
+        std::cerr << " Mesh area = " << meshArea(mesh) << " Soup area = " << soupArea(soup)
+                  << std::endl;
+        {
+            std::ofstream f(filename + "soup.ply");
+            CGAL::IO::write_PLY(f, soupToMesh(soup));
+        }
+    }
+
+    return 0;
+}