fix bug with flipEdgeNetwork created by piecewise Dijkstra

include/geometrycentral/surface/flip_geodesics.h

@@ -185,7 +185,6 @@ class FlipEdgeNetwork {
   SegmentAngleType locallyShortestTest(Halfedge hePrev, Halfedge heNext);
   double minWedgeAngle(Halfedge hePrev, Halfedge heNext);
   double minWedgeAngle(const FlipPathSegment& segment);
-  bool isStraight(double angleThresh = 1e-4);
   double minAngle();        // minimum over all angles
   double minAngleIsotopy(); // minimum over all angles, excluding those blocked by an path endpoint
   struct ShortestReturnBoth {

src/surface/flip_geodesics.cpp

@@ -308,7 +308,25 @@ std::unique_ptr<FlipEdgeNetwork> FlipEdgeNetwork::constructFromPiecewiseDijkstra
     halfedges.insert(halfedges.end(), dijkstraPath.begin(), dijkstraPath.end());
   }
 
-  return std::unique_ptr<FlipEdgeNetwork>(new FlipEdgeNetwork(mesh_, geom, {halfedges}, extraMark));
+  if (markInterior || halfedges.size() == 0) {
+    return std::unique_ptr<FlipEdgeNetwork>(new FlipEdgeNetwork(mesh_, geom, {halfedges}, extraMark));
+  }
+
+  // Sometimes Dijkstra segments to an intermediate point create a back-and-forth that the FlipOut algorithm
+  // cannot remove, because of an implementation limitation of how we pass it to the stacked representation. This is a
+  // quick filtering pass to collapse and remove such segments when they trivially exist
+  //
+  // NOTE: we can only do this when markInterior==false, because the path will no longer touch the marked vertex.
+  std::vector<Halfedge> heClean;
+  for (Halfedge he : halfedges) {
+    if (!heClean.empty() && heClean.back() == he.twin()) {
+      heClean.pop_back();
+    } else {
+      heClean.push_back(he);
+    }
+  }
+
+  return std::unique_ptr<FlipEdgeNetwork>(new FlipEdgeNetwork(mesh_, geom, {heClean}, extraMark));
 }
 
 

test/CMakeLists.txt

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

test/src/flip_geodesic_test.cpp

@@ -0,0 +1,170 @@
+#include "geometrycentral/surface/flip_geodesics.h"
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+#include "geometrycentral/surface/surface_mesh_factories.h"
+#include "geometrycentral/surface/vertex_position_geometry.h"
+
+#include "load_test_meshes.h"
+
+#include "gtest/gtest.h"
+
+#include <string>
+
+using namespace geometrycentral;
+using namespace geometrycentral::surface;
+
+class FlipGeodesicSuite : public MeshAssetSuite {};
+
+
+// ============================================================
+// =============== Construction
+// ============================================================
+
+TEST_F(FlipGeodesicSuite, ConstructFromDijkstra) {
+  for (const MeshAsset& a : {getAsset("fox.ply", true), getAsset("cat_head.obj", true)}) {
+    a.printThyName();
+    ManifoldSurfaceMesh& mesh = *a.manifoldMesh;
+    VertexPositionGeometry& geom = *a.geometry;
+
+    std::unique_ptr<FlipEdgeNetwork> network =
+        FlipEdgeNetwork::constructFromDijkstraPath(mesh, geom, mesh.vertex(0), mesh.vertex(mesh.nVertices() / 2));
+
+    ASSERT_NE(network, nullptr);
+    ASSERT_EQ(network->paths.size(), 1);
+    EXPECT_GT(network->paths[0]->size(), 0);
+  }
+}
+
+TEST_F(FlipGeodesicSuite, ConstructFromPiecewiseDijkstra) {
+  for (const MeshAsset& a : {getAsset("fox.ply", true), getAsset("cat_head.obj", true)}) {
+    a.printThyName();
+    ManifoldSurfaceMesh& mesh = *a.manifoldMesh;
+    VertexPositionGeometry& geom = *a.geometry;
+
+    size_t nV = mesh.nVertices();
+    std::vector<Vertex> waypoints = {mesh.vertex(0), mesh.vertex(nV / 3), mesh.vertex(2 * nV / 3)};
+
+    std::unique_ptr<FlipEdgeNetwork> network =
+        FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, geom, waypoints, false);
+
+    ASSERT_NE(network, nullptr);
+    ASSERT_EQ(network->paths.size(), 1);
+    EXPECT_GT(network->paths[0]->size(), 0);
+  }
+}
+
+TEST_F(FlipGeodesicSuite, ConstructClosedLoop) {
+  for (const MeshAsset& a : {getAsset("fox.ply", true), getAsset("cat_head.obj", true)}) {
+    a.printThyName();
+    ManifoldSurfaceMesh& mesh = *a.manifoldMesh;
+    VertexPositionGeometry& geom = *a.geometry;
+
+    size_t nV = mesh.nVertices();
+    std::vector<Vertex> waypoints = {mesh.vertex(0), mesh.vertex(nV / 3), mesh.vertex(2 * nV / 3)};
+
+    std::unique_ptr<FlipEdgeNetwork> network =
+        FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, geom, waypoints, /*closed=*/true);
+
+    ASSERT_NE(network, nullptr);
+    ASSERT_EQ(network->paths.size(), 1);
+    EXPECT_TRUE(network->paths[0]->isClosed);
+  }
+}
+
+
+// ============================================================
+// =============== Shortening
+// ============================================================
+
+TEST_F(FlipGeodesicSuite, ShortenedPathIsStraight) {
+  for (const MeshAsset& a : {getAsset("fox.ply", true), getAsset("cat_head.obj", true)}) {
+    a.printThyName();
+    ManifoldSurfaceMesh& mesh = *a.manifoldMesh;
+    VertexPositionGeometry& geom = *a.geometry;
+
+    std::unique_ptr<FlipEdgeNetwork> network =
+        FlipEdgeNetwork::constructFromDijkstraPath(mesh, geom, mesh.vertex(0), mesh.vertex(mesh.nVertices() / 2));
+    ASSERT_NE(network, nullptr);
+
+    network->iterativeShorten();
+
+    EXPECT_GE(network->minAngle(), M_PI - network->EPS_ANGLE);
+  }
+}
+
+TEST_F(FlipGeodesicSuite, ShorteningDoesNotIncreaseLength) {
+  for (const MeshAsset& a : {getAsset("fox.ply", true), getAsset("cat_head.obj", true)}) {
+    a.printThyName();
+    ManifoldSurfaceMesh& mesh = *a.manifoldMesh;
+    VertexPositionGeometry& geom = *a.geometry;
+
+    std::unique_ptr<FlipEdgeNetwork> network =
+        FlipEdgeNetwork::constructFromDijkstraPath(mesh, geom, mesh.vertex(0), mesh.vertex(mesh.nVertices() / 2));
+    ASSERT_NE(network, nullptr);
+
+    double lengthBefore = network->length();
+    network->iterativeShorten();
+    double lengthAfter = network->length();
+
+    EXPECT_LE(lengthAfter, lengthBefore + 1e-6);
+  }
+}
+
+// ============================================================
+// =============== Piecewise Dijkstra overlap removal
+// ============================================================
+
+// See: https://github.com/nmwsharp/geometry-central/pull/237
+
+// Small 5-triangle mesh for overlap-removal tests
+namespace {
+std::tuple<std::unique_ptr<ManifoldSurfaceMesh>, std::unique_ptr<VertexPositionGeometry>> buildOverlapTestMesh() {
+  std::vector<Vector3> verts = {{0, 0, 0}, {0, 3, 0}, {0, 3, 1}, {0, 5, 0}, {0, 5, 5}, {0, 3, 2}, {0, 5, 0}};
+  std::vector<std::vector<size_t>> faces = {{0, 1, 6}, {1, 3, 4}, {2, 1, 4}, {2, 4, 5}, {1, 2, 6}};
+  return makeManifoldSurfaceMeshAndGeometry(faces, verts);
+}
+} // namespace
+
+TEST_F(FlipGeodesicSuite, PiecewiseDijkstraSingleDetour) {
+  std::unique_ptr<ManifoldSurfaceMesh> meshPtr;
+  std::unique_ptr<VertexPositionGeometry> geomPtr;
+  std::tie(meshPtr, geomPtr) = buildOverlapTestMesh();
+  ManifoldSurfaceMesh& mesh = *meshPtr;
+  VertexPositionGeometry& geom = *geomPtr;
+
+  std::vector<Vertex> path = {mesh.vertex(0), mesh.vertex(2), mesh.vertex(3)};
+  std::unique_ptr<FlipEdgeNetwork> network =
+      FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, geom, path, false);
+
+  ASSERT_EQ(network->paths.size(), 1);
+  EXPECT_EQ(network->paths[0]->getHalfedgeList().size(), 2);
+}
+
+TEST_F(FlipGeodesicSuite, PiecewiseDijkstraMultipleOverlapAtStart) {
+  std::unique_ptr<ManifoldSurfaceMesh> meshPtr;
+  std::unique_ptr<VertexPositionGeometry> geomPtr;
+  std::tie(meshPtr, geomPtr) = buildOverlapTestMesh();
+  ManifoldSurfaceMesh& mesh = *meshPtr;
+  VertexPositionGeometry& geom = *geomPtr;
+
+  std::vector<Vertex> path = {mesh.vertex(1), mesh.vertex(5), mesh.vertex(3)};
+  std::unique_ptr<FlipEdgeNetwork> network =
+      FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, geom, path, false);
+
+  ASSERT_EQ(network->paths.size(), 1);
+  EXPECT_EQ(network->paths[0]->getHalfedgeList().size(), 1);
+}
+
+TEST_F(FlipGeodesicSuite, PiecewiseDijkstraMultipleOverlapAtEnd) {
+  std::unique_ptr<ManifoldSurfaceMesh> meshPtr;
+  std::unique_ptr<VertexPositionGeometry> geomPtr;
+  std::tie(meshPtr, geomPtr) = buildOverlapTestMesh();
+  ManifoldSurfaceMesh& mesh = *meshPtr;
+  VertexPositionGeometry& geom = *geomPtr;
+
+  std::vector<Vertex> path = {mesh.vertex(0), mesh.vertex(2), mesh.vertex(5), mesh.vertex(1)};
+  std::unique_ptr<FlipEdgeNetwork> network =
+      FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, geom, path, false);
+
+  ASSERT_EQ(network->paths.size(), 1);
+  EXPECT_EQ(network->paths[0]->getHalfedgeList().size(), 1);
+}