Exact geodesics: honour max_propagation_distance without stop points; require curvatures once

src/surface/exact_geodesics.cpp

@@ -444,6 +444,7 @@ void GeodesicAlgorithmExact::propagate(const std::vector<SurfacePoint>& sources,
 
   IntervalWithStop candidates[2];
   geom.requireCornerAngles();
+  geom.requireVertexGaussianCurvatures();
 
   while (!m_queue.empty()) {
     m_queue_max_size = std::max(m_queue.size(), m_queue_max_size);
@@ -469,9 +470,7 @@ void GeodesicAlgorithmExact::propagate(const std::vector<SurfacePoint>& sources,
     bool const last_interval = min_interval->next() == nullptr;
 
     auto saddleOrBoundary = [&](Vertex v) -> bool {
-      geom.requireVertexGaussianCurvatures();
       bool saddle = geom.vertexGaussianCurvatures[v] < 0;
-      geom.unrequireVertexGaussianCurvatures();
       return saddle || vertexIsBoundary[v];
     };
 
@@ -607,6 +606,7 @@ void GeodesicAlgorithmExact::propagate(const std::vector<SurfacePoint>& sources,
   m_time_consumed = (static_cast<double>(stop) - static_cast<double>(start)) / CLOCKS_PER_SEC;
 
   geom.unrequireCornerAngles();
+  geom.unrequireVertexGaussianCurvatures();
 
   /*	for(unsigned i=0; i<m_edge_interval_lists.size(); ++i)
     {
@@ -631,7 +631,8 @@ bool GeodesicAlgorithmExact::check_stop_conditions(unsigned& index) const {
   }
 
   if (m_stop_vertices.empty()) {
-    return false;
+    // without stop vertices, the distance bound alone stops the propagation
+    return stop_dist < GEODESIC_INF;
   }
 
   while (index < m_stop_vertices.size()) {

test/CMakeLists.txt

@@ -108,6 +108,7 @@ set(TEST_SRCS
   src/point_cloud_test.cpp
   src/poisson_disk_sampler_test.cpp
   src/polygon_operators_test.cpp
+  src/exact_geodesics_test.cpp
   src/flip_geodesic_test.cpp
   src/stl_reader_test.cpp
   src/surface_misc_test.cpp

test/src/exact_geodesics_test.cpp

@@ -0,0 +1,115 @@
+#include "geometrycentral/surface/exact_geodesics.h"
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+#include "geometrycentral/surface/vertex_position_geometry.h"
+
+#include "load_test_meshes.h"
+
+#include "gtest/gtest.h"
+
+#include <algorithm>
+#include <vector>
+
+using namespace geometrycentral;
+using namespace geometrycentral::surface;
+
+class ExactGeodesicSuite : public MeshAssetSuite {};
+
+namespace {
+
+// A distance that about half of the vertices are closer than
+double medianDistance(SurfaceMesh& mesh, const VertexData<double>& dist) {
+  std::vector<double> d;
+  for (Vertex v : mesh.vertices()) d.push_back(dist[v]);
+  std::nth_element(d.begin(), d.begin() + d.size() / 2, d.end());
+  return d[d.size() / 2];
+}
+
+} // namespace
+
+// The distance field does not depend on whether the caller holds the curvatures used by the saddle test
+TEST_F(ExactGeodesicSuite, ResultIndependentOfRequiredQuantities) {
+  for (const MeshAsset& a : {getAsset("spot.ply", true), getAsset("bob_small.ply", true)}) {
+    a.printThyName();
+    SurfaceMesh& mesh = *a.mesh;
+    VertexPositionGeometry& geometry = *a.geometry;
+    Vertex source = mesh.vertex(0);
+
+    GeodesicAlgorithmExact mmp(mesh, geometry);
+    mmp.propagate(source);
+    VertexData<double> dist = mmp.getDistanceFunction();
+
+    geometry.requireVertexGaussianCurvatures();
+    GeodesicAlgorithmExact mmpHeld(mesh, geometry);
+    mmpHeld.propagate(source);
+    VertexData<double> distHeld = mmpHeld.getDistanceFunction();
+    geometry.unrequireVertexGaussianCurvatures();
+
+    for (Vertex v : mesh.vertices()) {
+      EXPECT_EQ(dist[v], distHeld[v]);
+    }
+  }
+}
+
+// With a distance bound and no stop points, propagation stops at the bound: vertices closer than the bound
+// get their exact distance, and some vertices beyond it are left unreached
+TEST_F(ExactGeodesicSuite, MaxPropagationDistanceWithoutStopPoints) {
+  for (const MeshAsset& a : {getAsset("spot.ply", true), getAsset("bob_small.ply", true)}) {
+    a.printThyName();
+    SurfaceMesh& mesh = *a.mesh;
+    VertexPositionGeometry& geometry = *a.geometry;
+    Vertex source = mesh.vertex(0);
+
+    GeodesicAlgorithmExact full(mesh, geometry);
+    full.propagate(source);
+    VertexData<double> ref = full.getDistanceFunction();
+    double bound = medianDistance(mesh, ref);
+
+    GeodesicAlgorithmExact mmp(mesh, geometry);
+    mmp.propagate(source, bound);
+    VertexData<double> dist = mmp.getDistanceFunction();
+
+    size_t unreached = 0;
+    for (Vertex v : mesh.vertices()) {
+      if (ref[v] < bound) {
+        EXPECT_NEAR(dist[v], ref[v], 1e-12 * bound);
+      }
+      if (dist[v] >= GEODESIC_INF) unreached++;
+    }
+    EXPECT_GT(unreached, 0u);
+  }
+}
+
+// With stop points, propagation stops once they are covered; with a distance bound too, once both hold
+TEST_F(ExactGeodesicSuite, MaxPropagationDistanceWithStopPoints) {
+  for (const MeshAsset& a : {getAsset("spot.ply", true), getAsset("bob_small.ply", true)}) {
+    a.printThyName();
+    SurfaceMesh& mesh = *a.mesh;
+    VertexPositionGeometry& geometry = *a.geometry;
+    Vertex source = mesh.vertex(0);
+
+    GeodesicAlgorithmExact full(mesh, geometry);
+    full.propagate(source);
+    VertexData<double> ref = full.getDistanceFunction();
+    double bound = medianDistance(mesh, ref);
+
+    // A stop vertex well inside the bound
+    Vertex stop;
+    for (Vertex v : mesh.vertices()) {
+      if (v != source && ref[v] < 0.5 * bound && (stop == Vertex() || ref[v] > ref[stop])) stop = v;
+    }
+    ASSERT_TRUE(stop != Vertex());
+
+    GeodesicAlgorithmExact mmpStop(mesh, geometry);
+    mmpStop.propagate(source, GEODESIC_INF, {stop});
+    EXPECT_NEAR(mmpStop.getDistance(stop), ref[stop], 1e-12 * bound);
+
+    GeodesicAlgorithmExact mmpBoth(mesh, geometry);
+    mmpBoth.propagate(source, bound, {stop});
+    VertexData<double> dist = mmpBoth.getDistanceFunction();
+    for (Vertex v : mesh.vertices()) {
+      if (ref[v] < bound) {
+        EXPECT_NEAR(dist[v], ref[v], 1e-12 * bound);
+      }
+    }
+  }
+}