Add cotangent-based Laplace-Beltrami and mean curvature support to fv_surface_mesh

src/YgorMath.cc

@@ -6284,6 +6284,159 @@ fv_surface_mesh<T,I>::surface_area(int64_t n) const {
     template double fv_surface_mesh<double, uint64_t >::surface_area(int64_t) const;
 #endif
 
+namespace {
+    template <class T>
+    T
+    cotangent_for_opposing_vertex(const vec3<T> &P1, const vec3<T> &P2, const vec3<T> &P_opposite){
+        const auto A = P1 - P_opposite;
+        const auto B = P2 - P_opposite;
+        const auto denom = A.Cross(B).length();
+        if(denom <= std::numeric_limits<T>::epsilon()){
+            throw std::runtime_error("Unable to compute cotangent weight due to degenerate triangle. Cannot continue.");
+        }
+        return A.Dot(B) / denom;
+    }
+}
+
+template <class T, class I>
+std::array<T,2>
+fv_surface_mesh<T,I>::cotangent_weights(I v1, I v2, I v3, I v4) const {
+    const auto &P1 = this->vertices.at(v1);
+    const auto &P2 = this->vertices.at(v2);
+    const auto &P3 = this->vertices.at(v3);
+    const auto &P4 = this->vertices.at(v4);
+
+    return std::array<T,2>{{
+        cotangent_for_opposing_vertex(P1, P2, P3),
+        cotangent_for_opposing_vertex(P1, P2, P4)
+    }};
+}
+#ifndef YGORMATH_DISABLE_ALL_SPECIALIZATIONS
+    template std::array<float,2>  fv_surface_mesh<float , uint32_t >::cotangent_weights(uint32_t, uint32_t, uint32_t, uint32_t) const;
+    template std::array<float,2>  fv_surface_mesh<float , uint64_t >::cotangent_weights(uint64_t, uint64_t, uint64_t, uint64_t) const;
+
+    template std::array<double,2> fv_surface_mesh<double, uint32_t >::cotangent_weights(uint32_t, uint32_t, uint32_t, uint32_t) const;
+    template std::array<double,2> fv_surface_mesh<double, uint64_t >::cotangent_weights(uint64_t, uint64_t, uint64_t, uint64_t) const;
+#endif
+
+template <class T, class I>
+std::vector<vec3<T>>
+fv_surface_mesh<T,I>::laplace_beltrami_operator() const {
+    std::vector<vec3<T>> lb(this->vertices.size(), vec3<T>(0,0,0));
+    std::vector<T> barycentric_areas(this->vertices.size(), static_cast<T>(0));
+
+    std::map<std::pair<I,I>, std::vector<I>> edge_to_opposing_vertices;
+    for(const auto &fv : this->faces){
+        if(fv.size() < 3) continue;
+        if(fv.size() != 3){
+            throw std::runtime_error("Laplace-Beltrami operator requires a triangular surface mesh. Cannot continue.");
+        }
+        const auto i0 = fv.at(0);
+        const auto i1 = fv.at(1);
+        const auto i2 = fv.at(2);
+
+        if( (i0 == i1) || (i1 == i2) || (i0 == i2) ){
+            throw std::runtime_error("Laplace-Beltrami operator requires faces with three distinct vertices. Cannot continue.");
+        }
+        if( (this->vertices.size() <= static_cast<size_t>(i0))
+        ||  (this->vertices.size() <= static_cast<size_t>(i1))
+        ||  (this->vertices.size() <= static_cast<size_t>(i2)) ){
+            throw std::runtime_error("Face references a non-existent vertex. Cannot continue.");
+        }
+
+        const auto &P0 = this->vertices.at(i0);
+        const auto &P1 = this->vertices.at(i1);
+        const auto &P2 = this->vertices.at(i2);
+        const auto area = (P1 - P0).Cross(P2 - P0).length() / static_cast<T>(2);
+        if(area <= std::numeric_limits<T>::epsilon()){
+            throw std::runtime_error("Laplace-Beltrami operator requires non-degenerate triangles. Cannot continue.");
+        }
+
+        const auto area_contrib = area / static_cast<T>(3);
+        barycentric_areas.at(i0) += area_contrib;
+        barycentric_areas.at(i1) += area_contrib;
+        barycentric_areas.at(i2) += area_contrib;
+
+        const auto key01 = (i0 < i1) ? std::make_pair(i0, i1) : std::make_pair(i1, i0);
+        const auto key12 = (i1 < i2) ? std::make_pair(i1, i2) : std::make_pair(i2, i1);
+        const auto key20 = (i2 < i0) ? std::make_pair(i2, i0) : std::make_pair(i0, i2);
+
+        edge_to_opposing_vertices[key01].push_back(i2);
+        edge_to_opposing_vertices[key12].push_back(i0);
+        edge_to_opposing_vertices[key20].push_back(i1);
+    }
+
+    for(const auto &edge_opps : edge_to_opposing_vertices){
+        const auto i = edge_opps.first.first;
+        const auto j = edge_opps.first.second;
+        const auto &opps = edge_opps.second;
+
+        T w = static_cast<T>(0);
+        if(opps.size() == 1){
+            const auto &P_i = this->vertices.at(i);
+            const auto &P_j = this->vertices.at(j);
+            const auto &P_opp = this->vertices.at(opps.at(0));
+            w = cotangent_for_opposing_vertex(P_i, P_j, P_opp);
+        }else if(opps.size() == 2){
+            const auto cot_w = this->cotangent_weights(i, j, opps.at(0), opps.at(1));
+            w = cot_w.at(0) + cot_w.at(1);
+        }else{
+            throw std::runtime_error("Edge incident to more than two triangles. Laplace-Beltrami operator requires a manifold triangular surface mesh. Cannot continue.");
+        }
+
+        const auto edge_ij = this->vertices.at(j) - this->vertices.at(i);
+        lb.at(i) += edge_ij * w;
+        lb.at(j) -= edge_ij * w;
+    }
+
+    for(size_t i = 0; i < lb.size(); ++i){
+        const auto area = barycentric_areas.at(i);
+        if(area <= std::numeric_limits<T>::epsilon()){
+            throw std::runtime_error("Laplace-Beltrami operator requires every vertex to have positive local area and to be referenced by at least one face. Vertices with zero local area are invalid. Cannot continue.");
+        }
+        lb.at(i) /= (static_cast<T>(2) * area);
+    }
+
+    return lb;
+}
+#ifndef YGORMATH_DISABLE_ALL_SPECIALIZATIONS
+    template std::vector<vec3<float >> fv_surface_mesh<float , uint32_t >::laplace_beltrami_operator() const;
+    template std::vector<vec3<float >> fv_surface_mesh<float , uint64_t >::laplace_beltrami_operator() const;
+
+    template std::vector<vec3<double>> fv_surface_mesh<double, uint32_t >::laplace_beltrami_operator() const;
+    template std::vector<vec3<double>> fv_surface_mesh<double, uint64_t >::laplace_beltrami_operator() const;
+#endif
+
+template <class T, class I>
+std::vector<T>
+fv_surface_mesh<T,I>::mean_curvature(const std::vector<vec3<T>> &laplace_beltrami) const {
+    std::vector<T> out;
+    out.reserve(laplace_beltrami.size());
+    for(const auto &lbi : laplace_beltrami){
+        out.push_back( lbi.length() / static_cast<T>(2) );
+    }
+    return out;
+}
+
+template <class T, class I>
+std::vector<T>
+fv_surface_mesh<T,I>::mean_curvature() const {
+    const auto lb = this->laplace_beltrami_operator();
+    return this->mean_curvature(lb);
+}
+#ifndef YGORMATH_DISABLE_ALL_SPECIALIZATIONS
+    template std::vector<float > fv_surface_mesh<float , uint32_t >::mean_curvature(const std::vector<vec3<float >> &) const;
+    template std::vector<float > fv_surface_mesh<float , uint64_t >::mean_curvature(const std::vector<vec3<float >> &) const;
+
+    template std::vector<double> fv_surface_mesh<double, uint32_t >::mean_curvature(const std::vector<vec3<double>> &) const;
+    template std::vector<double> fv_surface_mesh<double, uint64_t >::mean_curvature(const std::vector<vec3<double>> &) const;
+    template std::vector<float > fv_surface_mesh<float , uint32_t >::mean_curvature() const;
+    template std::vector<float > fv_surface_mesh<float , uint64_t >::mean_curvature() const;
+
+    template std::vector<double> fv_surface_mesh<double, uint32_t >::mean_curvature() const;
+    template std::vector<double> fv_surface_mesh<double, uint64_t >::mean_curvature() const;
+#endif
+
 // Regenerates this->involved_faces using this->vertices and this->faces.
 template <class T, class I>
 void

src/YgorMath.h

@@ -593,6 +593,27 @@ template <class T, class I>   class fv_surface_mesh {
         // all faces.
         T surface_area(int64_t n = -1) const;
 
+        // Computes cotangent weights for the edge (v1, v2) in a triangle mesh.
+        //
+        // Topology convention:
+        //   - (v1, v2) must be an edge shared by two adjacent triangles in the mesh.
+        //   - v3 is the vertex opposite this edge in one triangle, so that (v1, v2, v3) forms a face.
+        //   - v4 is the vertex opposite this edge in the other adjacent triangle, so that (v2, v1, v4) forms a face.
+        //
+        // In other words, v3 and v4 are the third vertices of the two faces incident on the edge (v1, v2). The returned
+        // array contains the cotangents of the angles at v3 and v4, respectively.
+        //
+        // Example of how to identify v3 and v4:
+        //   - Find all faces that contain both v1 and v2.
+        //   - For each such face, the remaining vertex index (neither v1 nor v2) is one of v3 or v4.
+        std::array<T,2> cotangent_weights(I v1, I v2, I v3, I v4) const;
+
+        // Applies the cotangent-weight Laplace-Beltrami operator to all vertices.
+        std::vector<vec3<T>> laplace_beltrami_operator() const;
+
+        // Estimates per-vertex mean curvature magnitudes from the Laplace-Beltrami operator.
+        std::vector<T> mean_curvature() const;
+
         // Regenerates this->involved_faces using this->vertices and this->faces.
         void recreate_involved_face_index(void);
 

tests2/YgorMath/fv_surface_mesh.cc

@@ -122,6 +122,82 @@ TEST_CASE( "fv_surface_mesh class" ){
         REQUIRE( mesh1.surface_area() == 0.5 );
     }
 
+    SUBCASE("cotangent_weights"){
+        fv_surface_mesh<double, uint32_t> mesh2;
+        mesh2.vertices = {{ vec3<double>( 1.0,  1.0,  1.0),
+                            vec3<double>(-1.0, -1.0,  1.0),
+                            vec3<double>(-1.0,  1.0, -1.0),
+                            vec3<double>( 1.0, -1.0, -1.0) }};
+        mesh2.faces = {{ static_cast<uint32_t>(0), static_cast<uint32_t>(1), static_cast<uint32_t>(2) },
+                       { static_cast<uint32_t>(0), static_cast<uint32_t>(3), static_cast<uint32_t>(1) },
+                       { static_cast<uint32_t>(0), static_cast<uint32_t>(2), static_cast<uint32_t>(3) },
+                       { static_cast<uint32_t>(1), static_cast<uint32_t>(3), static_cast<uint32_t>(2) }};
+
+        const auto cots = mesh2.cotangent_weights(0U, 1U, 2U, 3U);
+        const auto expected = 1.0 / std::sqrt(3.0);
+        REQUIRE( std::abs(cots.at(0) - expected) < eps );
+        REQUIRE( std::abs(cots.at(1) - expected) < eps );
+    }
+
+    SUBCASE("laplace_beltrami_operator and mean_curvature"){
+        fv_surface_mesh<double, uint32_t> mesh2;
+        mesh2.vertices = {{ vec3<double>( 1.0,  1.0,  1.0),
+                            vec3<double>(-1.0, -1.0,  1.0),
+                            vec3<double>(-1.0,  1.0, -1.0),
+                            vec3<double>( 1.0, -1.0, -1.0) }};
+        mesh2.faces = {{ static_cast<uint32_t>(0), static_cast<uint32_t>(1), static_cast<uint32_t>(2) },
+                       { static_cast<uint32_t>(0), static_cast<uint32_t>(3), static_cast<uint32_t>(1) },
+                       { static_cast<uint32_t>(0), static_cast<uint32_t>(2), static_cast<uint32_t>(3) },
+                       { static_cast<uint32_t>(1), static_cast<uint32_t>(3), static_cast<uint32_t>(2) }};
+
+        const auto lb = mesh2.laplace_beltrami_operator();
+        REQUIRE( lb.size() == mesh2.vertices.size() );
+        for(size_t i = 0; i < lb.size(); ++i){
+            const auto expected = mesh2.vertices.at(i) * (-2.0 / 3.0);
+            REQUIRE( std::abs(lb.at(i).x - expected.x) < eps );
+            REQUIRE( std::abs(lb.at(i).y - expected.y) < eps );
+            REQUIRE( std::abs(lb.at(i).z - expected.z) < eps );
+        }
+
+        const auto H = mesh2.mean_curvature();
+        REQUIRE( H.size() == mesh2.vertices.size() );
+        const auto expected_H = 1.0 / std::sqrt(3.0);
+        for(const auto &h : H){
+            REQUIRE( std::abs(h - expected_H) < eps );
+        }
+    }
+
+    SUBCASE("laplace_beltrami_operator throws for non-triangular facets"){
+        fv_surface_mesh<double, uint32_t> mesh2;
+        mesh2.vertices = {{ p1, p2, p3, p4 }};
+        mesh2.faces = {{ static_cast<uint32_t>(0), static_cast<uint32_t>(1), static_cast<uint32_t>(2), static_cast<uint32_t>(3) }};
+        REQUIRE_THROWS( mesh2.laplace_beltrami_operator() );
+        REQUIRE_THROWS( mesh2.mean_curvature() );
+    }
+
+    SUBCASE("laplace_beltrami_operator skips faces with fewer than 3 vertices"){
+        // Face with a single vertex
+        fv_surface_mesh<double, uint32_t> mesh_single;
+        mesh_single.vertices = {{ p1 }};
+        mesh_single.faces = {{ static_cast<uint32_t>(0) }};
+        REQUIRE_NOTHROW( mesh_single.laplace_beltrami_operator() );
+        REQUIRE_NOTHROW( mesh_single.mean_curvature() );
+        const auto lb_single = mesh_single.laplace_beltrami_operator();
+        const auto H_single = mesh_single.mean_curvature();
+        REQUIRE( lb_single.size() == mesh_single.vertices.size() );
+        REQUIRE( H_single.size() == mesh_single.vertices.size() );
+
+        // Face with two vertices
+        fv_surface_mesh<double, uint32_t> mesh_two;
+        mesh_two.vertices = {{ p1, p2 }};
+        mesh_two.faces = {{ static_cast<uint32_t>(0), static_cast<uint32_t>(1) }};
+        REQUIRE_NOTHROW( mesh_two.laplace_beltrami_operator() );
+        REQUIRE_NOTHROW( mesh_two.mean_curvature() );
+        const auto lb_two = mesh_two.laplace_beltrami_operator();
+        const auto H_two = mesh_two.mean_curvature();
+        REQUIRE( lb_two.size() == mesh_two.vertices.size() );
+        REQUIRE( H_two.size() == mesh_two.vertices.size() );
+    }
     SUBCASE("recreate_involved_face_index"){
         REQUIRE( mesh1.involved_faces.size() == 0 );
         mesh1.recreate_involved_face_index();
@@ -831,4 +907,3 @@ TEST_CASE( "Convex_Hull" ){
     }
 }
 
-