fixed bug in geodesic distances, adjusted unit test

src/pmp/algorithms/SurfaceGeodesic.cpp

@@ -16,21 +16,24 @@ namespace pmp {
 
 //=============================================================================
 
-SurfaceGeodesic::SurfaceGeodesic(SurfaceMesh& mesh, std::vector<Vertex> seed,
-                                 Scalar maxdist, bool use_virtual_edges)
+SurfaceGeodesic::SurfaceGeodesic(SurfaceMesh& mesh, 
+                                 bool use_virtual_edges)
     : mesh_(mesh),
-      seed_(seed),
-      maxdist_(maxdist),
       use_virtual_edges_(use_virtual_edges)
 {
-    distance_ = mesh_.add_vertex_property<Scalar>("geodesic:distance");
+    distance_  = mesh_.add_vertex_property<Scalar>("geodesic:distance");
     processed_ = mesh_.add_vertex_property<bool>("geodesic:processed");
-
-    front_ = new PriorityQueue(VertexCmp(distance_));
     find_virtual_edges();
-    init_front();
-    propagate_front();
-    delete front_;
+}
+
+//-----------------------------------------------------------------------------
+
+SurfaceGeodesic::SurfaceGeodesic(SurfaceMesh& mesh, 
+                                 const std::vector<Vertex>& seed,
+                                 Scalar maxdist, bool use_virtual_edges)
+    : SurfaceGeodesic(mesh, use_virtual_edges)
+{
+    compute(seed, maxdist);
 }
 
 //-----------------------------------------------------------------------------
@@ -43,6 +46,19 @@ SurfaceGeodesic::~SurfaceGeodesic()
 
 //-----------------------------------------------------------------------------
 
+void SurfaceGeodesic::compute(const std::vector<Vertex>& seed,
+                              Scalar maxdist)
+{
+    seed_    = seed;
+    maxdist_ = maxdist;
+    front_   = new PriorityQueue(VertexCmp(distance_));
+    init_front();
+    propagate_front();
+    delete front_;
+}
+
+//-----------------------------------------------------------------------------
+
 void SurfaceGeodesic::find_virtual_edges()
 {
     Halfedge hh, hhh;
@@ -167,14 +183,14 @@ void SurfaceGeodesic::init_front()
     for (auto v : mesh_.vertices())
     {
         processed_[v] = false;
-        distance_[v] = FLT_MAX;
+        distance_[v]  = FLT_MAX;
     }
 
     // initialize seed vertices
     for (auto v : seed_)
     {
         processed_[v] = true;
-        distance_[v] = 0.0;
+        distance_[v]  = 0.0;
     }
 
     // initialize seed's one-ring
@@ -186,7 +202,7 @@ void SurfaceGeodesic::init_front()
                 pmp::distance(mesh_.position(v), mesh_.position(vv));
             if (dist < distance_[vv])
             {
-                distance_[vv] = dist;
+                distance_[vv]  = dist;
                 processed_[vv] = true;
             }
         }
@@ -370,21 +386,6 @@ Scalar SurfaceGeodesic::distance(Vertex v0, Vertex v1, Vertex v2, Scalar r0,
     if (c < 0.0)
         return dykstra;
 
-    // Novotni: intersect two circles
-    // use Novotni when distances are not too large
-    const double l = pmp::distance(A, B);
-    if (std::max(TA, TB) / l < 10)
-    {
-        if (valid_triangle(l, a, b) && valid_triangle(l, TA, TB))
-        {
-            const double x2 = 0.5f * (l * l + b * b - a * a) / l;
-            const double y2 = sqrt(b * b - x2 * x2);
-            const double x = 0.5f * (l * l + TA * TA - TB * TB) / l;
-            const double y = -sqrt(TA * TA - x * x);
-            return pmp::distance(dvec2(x2, y2), dvec2(x, y));
-        }
-    }
-
     // Kimmel: solve quadratic equation
     const double u = TB - TA;
     const double aa = a * a + b * b - 2.0 * a * b * c;

src/pmp/algorithms/SurfaceGeodesic.h

@@ -29,16 +29,25 @@ namespace pmp {
 //!
 //! The methods works by a Dykstra-like breadth first traversal from
 //! the seed vertices, implemented by a head structure.
-//! See \cite kimmel_1998_geodesic and \cite novotni_2002_geodesic for details.
+//! See \cite kimmel_1998_geodesic for details.
 class SurfaceGeodesic
 {
 public:
+    //! Empty constructor. Computes virtual edges only (if set to true).
+    //! Call compute() to compute geodesic distances.
+    SurfaceGeodesic(SurfaceMesh& mesh, bool use_virtual_edges = true);
+
     //! Construct with mesh and seed vertices.
     //! Grow around seed up to specified maximum distance.
     //! Set whether to use virtual edges (more computation, more accurate result)
-    SurfaceGeodesic(SurfaceMesh& mesh, std::vector<Vertex> seed,
+    SurfaceGeodesic(SurfaceMesh& mesh, const std::vector<Vertex>& seed,
                     Scalar maxdist = FLT_MAX, bool use_virtual_edges = true);
 
+    //! Compute geodesic distances from specified seed points
+    //! up to the specified maximum distance.
+    void compute(const std::vector<Vertex>& seed,
+                 Scalar maxdist = FLT_MAX);
+
     // destructor
     ~SurfaceGeodesic();
 

tests/SurfaceGeodesicTest.cpp

@@ -30,24 +30,49 @@ TEST(SurfaceGeodesicTest, geodesic)
     Scalar d(0);
     for (auto v: mesh.vertices())
         d = std::max(d, geodist(v));
-    EXPECT_FLOAT_EQ(d, 3.13061);
+    EXPECT_FLOAT_EQ(d, 3.1348989);
+
+    // map distances to texture coordinates
+    geodist.distance_to_texture_coordinates();
+    auto tex = mesh.vertex_property<TexCoord>("v:tex");
+    EXPECT_TRUE(tex);
 }
 
-TEST(SurfaceGeodesicTest, geodesic_to_texture)
+TEST(SurfaceGeodesicTest, geodesic_symmetry)
 {
     // read irregular mesh (to have virtual edges)
     SurfaceMesh mesh;
     EXPECT_TRUE(mesh.read("pmp-data/off/bunny_adaptive.off"));
 
-    // use first vertex as seed
+    SurfaceGeodesic geodist(mesh);
     std::vector<Vertex> seed;
-    seed.push_back(Vertex(0));
+    Vertex v0, v1;
+    Scalar d0, d1;
 
-    // compute geodesic distance
-    SurfaceGeodesic geodist(mesh, seed);
+    // grow from first vector
+    v0 = Vertex(0);
+    seed.clear();
+    seed.push_back(v0);
+    geodist.compute(seed);
 
-    // map distances to texture coordinates
-    geodist.distance_to_texture_coordinates();
-    auto tex = mesh.vertex_property<TexCoord>("v:tex");
-    EXPECT_TRUE(tex);
+    // find maximum geodesic distance
+    d0=0;
+    for (auto v: mesh.vertices())
+    {
+        if (geodist(v) > d0)
+        {
+            d0 = geodist(v);
+            v1 = v;
+        }
+    }
+
+    // grow back from max-dist vertex to vertex 0
+    seed.clear();
+    seed.push_back(v1);
+    geodist.compute(seed);
+    d1 = geodist(v0);
+
+    // expect both distance to be the same
+    Scalar err = fabs(d0-d1) / (0.5*(d0+d1));
+    EXPECT_LT(err, 0.001);
 }