Improving projection API to converge to MeshDomain

Author: fprotais

Mesh_smoothing_3/doc/Mesh_smoothing_3/Concepts/C3t3Projector.h

@@ -72,7 +72,7 @@ using Curve_edge = `std::pair<Curve_index, Edge>`;
 /*!
 Return the plane the patch face should align to. 
 */
-`std::pair<Point_3, Vector_3>` patch_projection_plane(Patch_face patch_face, Point_3 face_center, double face_radius) const;
+std::pair<Point_3, Vector_3> patch_projection_plane(Patch_face patch_face, std::vector<Point_3> face_points) const;
 
 /*!
 Return if a patch face should be projected or not. 
@@ -83,7 +83,7 @@ bool project_patch_face(Patch_face patch_face) const;
 /*!
 Return the line the curve edge should align to. 
 */
-std::pair<Point_3, Vector_3> curve_projection_tangent(Curve_edge curve_edge, Point_3 edge_center, double segment_size) const;
+std::pair<Point_3, Vector_3> curve_projection_tangent(Curve_edge curve_edge, std::array<Point_3,2> edge_points) const;
 
 /*!
 Return if a curve edge should be projected or not. 

Mesh_smoothing_3/include/CGAL/Mesh_smoothing_3/boundary_aware_mesh_smoothing.h

@@ -21,6 +21,7 @@
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits_3.h>
 #include <CGAL/AABB_primitive.h>
+#include <CGAL/centroid.h>
 
 #include <map>
 #include <vector>
@@ -121,6 +122,16 @@ get(const Edge_to_point_property_map<C3t3>&, const typename C3t3::Edge& e)
 } // namespace internal
 
 
+
+template <typename Point3,
+          typename Vector3>
+struct Projection_data {
+    Point3 origin = Point3();
+    Vector3 vector = Vector3();
+    bool activated = true;
+    double weight = 1.;
+};
+
 /*!
 * \ingroup pkgMeshSmoothing3Projection
 *
@@ -143,36 +154,39 @@ class C3t3_mesh_projector {
 public:
     using Point_3 = typename K::Point_3;
     using Vector_3 = typename K::Vector_3;
+    using Projection = Projection_data<Point_3, Vector_3>;
 
     using Surface_patch_index = typename C3t3::Surface_patch_index;
     using Facet = typename C3t3::Facet;
     using Patch_face = std::pair<Surface_patch_index, Facet>;
 
-    std::pair<Point_3, Vector_3> patch_projection_plane(Patch_face patch_face, Point_3 face_center, double /*face_radius*/) const {
+    Projection patch_projection_plane(Patch_face patch_face, std::vector<Point_3> const &face_points) const {
+        Projection projection;
+        Point_3 face_center = CGAL::centroid(face_points.begin(), face_points.end());
         auto res = _facet_trees.at(patch_face.first).closest_point_and_primitive(face_center);
-        Point_3 closest_point = res.first;
+        projection.origin = res.first;
         const auto triangle = _c3t3->triangulation().triangle(res.second);
-        Vector_3 normal = CGAL::unit_normal(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2));
-        return {closest_point, normal};
+        projection.vector = CGAL::unit_normal(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2));
+        return projection;
     }
-    bool project_patch_face(Patch_face) const { return true; }
 
     using Curve_index = typename C3t3::Curve_index;
     using Edge = typename C3t3::Edge;
     using Curve_edge = std::pair<Curve_index, Edge>;
 
-    std::pair<Point_3, Vector_3> curve_projection_tangent(Curve_edge curve_edge, Point_3 edge_center, double /*segment_size*/) const {
+    Projection curve_projection_tangent(Curve_edge curve_edge, std::array<Point_3,2> const &edge_points) const {
+        Projection projection;
+        Point_3 edge_center = CGAL::midpoint(edge_points[0], edge_points[1]);
         auto res = _edge_trees.at(curve_edge.first).closest_point_and_primitive(edge_center);
-        Point_3 closest_point = res.first;
+        projection.origin =  res.first;
         const auto segment = _c3t3->triangulation().segment(res.second);
-        Vector_3 direction = (segment.target() - segment.source());
-        if (direction.squared_length() > 1e-8) {
-            direction /= CGAL::sqrt(direction.squared_length());
+        projection.vector = (segment.target() - segment.source());
+        if (projection.vector.squared_length() > 1e-8) {
+            projection.vector /= CGAL::sqrt(projection.vector.squared_length());
         }
-        return {closest_point, direction};
+        return projection;
     }
 
-    bool project_curve_edge(Curve_edge) const { return true; }
 
 public:
 
@@ -264,26 +278,82 @@ class C3t3_no_projection {
 public:
     using Point_3 = typename K::Point_3;
     using Vector_3 = typename K::Vector_3;
+    using Projection = Projection_data<Point_3, Vector_3>;
 
     using Surface_patch_index = typename C3t3::Surface_patch_index;
     using Facet = typename C3t3::Facet;
     using Patch_face = std::pair<Surface_patch_index, Facet>;
 
-    std::pair<Point_3, Vector_3> patch_projection_plane(Patch_face, Point_3, double) const {
-        return {Point_3(), Vector_3()};
+    Projection patch_projection_plane(Patch_face, std::vector<Point_3> const &) const {
+        Point_3 face_center = CGAL::centroid(face_points.begin(), face_points.end());
+        Vector_3 normal = CGAL::unit_normal(face_points[0], face_points[1], face_points[2]); // only works for triangles, but let's start with that
+        Ray_3 ray_positive(face_center,  normal);
+        Ray_3 ray_negative(face_center, -normal);
+        return Projection{Point_3(), Vector_3(), false};
     }
-    bool project_patch_face(Patch_face) const { return false; }
 
     using Curve_index = typename C3t3::Curve_index;
     using Edge = typename C3t3::Edge;
     using Curve_edge = std::pair<Curve_index, Edge>;
 
-    std::pair<Point_3, Vector_3> curve_projection_tangent(Curve_edge, Point_3, double) const {
-        return {Point_3(), Vector_3()};
+    Projection curve_projection_tangent(Curve_edge, std::array<Point_3,2> const &) const {
+        return Projection{Point_3(), Vector_3(), false};
     }
-    bool project_curve_edge(Curve_edge) const { return false; }
 
 };
+
+
+/*
+* \ingroup pkgMeshSmoothing3Projection
+*
+* \brief provides projection to a Mesh Domain
+*
+* @tparam C3t3 model of `MeshDomain_3`
+*
+* \cgalModels{C3t3Projector}
+*
+\sa `CGAL::boundary_aware_mesh_smoothing`
+*
+*/
+template<typename C3t3, typename MeshDomain>
+class Mesh_domain_projection {
+    using K = typename C3t3::Triangulation::Geom_traits;
+public:
+    using Point_3 = typename K::Point_3;
+    using Vector_3 = typename K::Vector_3;
+    using Projection = Projection_data<Point_3, Vector_3>;
+
+    using Surface_patch_index = typename C3t3::Surface_patch_index;
+    using Facet = typename C3t3::Facet;
+    using Patch_face = std::pair<Surface_patch_index, Facet>;
+
+    Projection patch_projection_plane(Patch_face, std::vector<Point_3> const &) const {
+        return Projection();
+    }
+
+    using Curve_index = typename C3t3::Curve_index;
+    using Edge = typename C3t3::Edge;
+    using Curve_edge = std::pair<Curve_index, Edge>;
+
+    Projection curve_projection_tangent(Curve_edge, std::array<Point_3,2> const &) const {
+        return Projection();
+    }
+
+public: 
+
+    /*
+        Class constructor
+
+        \param mesh_domain contains the domain used for projection
+
+    */
+    Mesh_domain_projection(MeshDomain const& mesh_domain)
+    : _domain(mesh_domain)
+    {}
+
+    MeshDomain const& _domain;
+};
+
 } // namespace Mesh_smoothing_3
 
 /*!
@@ -415,16 +485,16 @@ void boundary_aware_mesh_smoothing  (
     using Vector_3 = typename  C3t3::Triangulation::Geom_traits::Vector_3;
 
     // converting Projector to Mesh_smoothing_3 queries
-    smoother.set_boundary_query([&](const Point_3& pt, std::pair<typename C3t3::Surface_patch_index, typename C3t3::Facet> patch_face, double radius) {
-        if (!projector.project_patch_face(patch_face)) return std::make_tuple(Point_3(0.,0.,0.), Vector_3{1.,0.,0.}, 0.);
-        auto [point, normal] = projector.patch_projection_plane(patch_face, pt, radius);
-        return std::make_tuple(point, normal, 1.0);
+    smoother.set_boundary_query([&](std::vector<Point_3> const& pts, std::pair<typename C3t3::Surface_patch_index, typename C3t3::Facet> patch_face) {
+        Mesh_smoothing_3::Projection_data<Point_3, Vector_3> proj = projector.patch_projection_plane(patch_face, pts);
+        if (!proj.activated) proj.weight = 0.;
+        return std::make_tuple(proj.origin, proj.vector, proj.weight);
     });
 
-    smoother.set_curves_query([&](const Point_3& pt, std::pair<typename C3t3::Curve_index, typename C3t3::Edge> curve_edge, double segment_size) {
-        if (!projector.project_curve_edge(curve_edge)) return std::make_tuple(Point_3(0.,0.,0.), Vector_3{1.,0.,0.}, 0.);
-        auto [point, tangent] = projector.curve_projection_tangent(curve_edge, pt, segment_size);
-        return std::make_tuple(point, tangent, 1.0);
+    smoother.set_curves_query([&](std::array<Point_3, 2> const& pts, std::pair<typename C3t3::Curve_index, typename C3t3::Edge> curve_edge) {
+        Mesh_smoothing_3::Projection_data<Point_3, Vector_3> proj = projector.curve_projection_tangent(curve_edge, pts);
+        if (!proj.activated) proj.weight = 0.;
+        return std::make_tuple(proj.origin, proj.vector, proj.weight);
     });
 
     smoother.set_verbose(verbose);

Mesh_smoothing_3/include/CGAL/Mesh_smoothing_3/internal/Mesh_smoothing_3_impl.h

@@ -272,6 +272,9 @@ void Mesh_smoother<TetrahedralMesh, BoundaryMesh, EdgeNetwork>::create_compress_
         }
         return (*res.first).second;
     };
+
+    // initialize volume
+
     _tetrahedra.reserve(_mesh.nb_cells());
     _tetrahedron_refs.reserve(_mesh.nb_cells());
     for (auto cell : _mesh.cell_range()) {
@@ -300,7 +303,6 @@ void Mesh_smoother<TetrahedralMesh, BoundaryMesh, EdgeNetwork>::create_compress_
 
     _tetrahedra.shrink_to_fit();
     _tetrahedron_refs.shrink_to_fit();
-    _compressed_coords = Eigen::Map<Eigen::VectorXd>(temp_coordinates_storage.data(), static_cast<Eigen::Index>(temp_coordinates_storage.size()));
 
     // computing vertex -> tet data structure
     std::vector<unsigned> curr_id(nb_points,0);
@@ -314,8 +316,61 @@ void Mesh_smoother<TetrahedralMesh, BoundaryMesh, EdgeNetwork>::create_compress_
             _vert2tet_corner[v][curr_id[v]++] = 4*t+i;
         }
     }
-    initialize_boundary();
-    initialize_curve_network();
+
+    // initialize boundary
+    
+    std::vector<unsigned> currFace(10);
+    _bnd_faces.reserve(_boundary.nb_faces());
+    _face_surface_id.reserve(_boundary.nb_faces());
+    std::unordered_map<unsigned, std::vector<std::array<unsigned, 2>>> vert2faces; // todo: maybe use a map from boost
+    for (auto face : _boundary.face_range()) {
+        bool hasUnlocked = false;
+        for (auto vertex_descriptor : _boundary.face_vertices(face)) {
+            if (!is_vert_locked(vertex_descriptor)) {
+                hasUnlocked = true;
+                break;
+            }
+        }
+        if (!hasUnlocked) continue;
+        currFace.clear();
+        for (auto vertex_descriptor : _boundary.face_vertices(face)) {
+            currFace.push_back(get_compressed_point_id(vertex_descriptor));
+        }
+        for (unsigned i = 0; i < currFace.size(); ++i) {
+            if (_lock_boundary) {
+                for (unsigned d = 0; d < 3; ++d) {
+                    _compressed_locks[3*currFace[i]+d] = true;
+                }
+            }
+            std::array<unsigned, 2> location_pair = {static_cast<unsigned>(_bnd_faces.size()), i};
+            auto res = vert2faces.emplace(currFace[i], std::vector<std::array<unsigned, 2>>{location_pair});
+            if (!res.second) (*res.first).second.push_back(location_pair);
+        }
+        _bnd_faces.push_back(currFace);
+        _face_surface_id.push_back(_boundary.patch_id(face));
+    }
+    // not ordered, do a sort? then I lose memory alignment? check if becomes bottleneck
+    _vert_and_face_corners.assign(vert2faces.begin(), vert2faces.end());
+
+
+    // initialize curve network
+
+    _curve_edges.reserve(_edge_network.nb_edges());
+    _curve_ids.reserve(_edge_network.nb_edges());
+    for (auto edge : _edge_network.edge_range()) {
+        Vertex_descriptor v0 = _edge_network.edge_vertex(edge, 0);
+        Vertex_descriptor v1 = _edge_network.edge_vertex(edge, 1);
+        if (is_vert_locked(v0) && is_vert_locked(v1)) continue;
+        _curve_edges.push_back({get_compressed_point_id(v0), get_compressed_point_id(v1)});
+        _curve_ids.push_back(_edge_network.curve_id(edge));
+    }
+
+
+
+    // initialize rest of data
+    _point_targets.reserve(_vertex_target_positions.size());
+
+    _compressed_coords = Eigen::Map<Eigen::VectorXd>(temp_coordinates_storage.data(), static_cast<Eigen::Index>(temp_coordinates_storage.size()));
 
     if (_update_validator != nullptr) {
         _current_coords_to_check.reserve(_mesh.nb_vertices());