fix bug with flipEdgeNetwork created by piecewise Dijkstra

include/geometrycentral/surface/flip_geodesics.h

@@ -107,7 +107,7 @@ class FlipEdgeNetwork {
   // Run Dijkstra between i'th and (i+1)'th point to initialize path
   static std::unique_ptr<FlipEdgeNetwork>
   constructFromPiecewiseDijkstraPath(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom,
-                                     std::vector<Vertex> points, bool closed = false, bool markInterior = false);
+                                     std::vector<Vertex> points, bool closed = false, bool markInterior = false, bool removeOverlap = true);
 
   // Consturct path(s) from marked edges, heuristically inferring endpoints, loopiness, etc
   static std::unique_ptr<FlipEdgeNetwork> constructFromEdgeSet(ManifoldSurfaceMesh& mesh,

src/surface/flip_geodesics.cpp

@@ -283,7 +283,7 @@ std::unique_ptr<FlipEdgeNetwork> FlipEdgeNetwork::constructFromDijkstraPath(Mani
 std::unique_ptr<FlipEdgeNetwork> FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(ManifoldSurfaceMesh& mesh_,
                                                                                      IntrinsicGeometryInterface& geom,
                                                                                      std::vector<Vertex> points,
-                                                                                     bool closed, bool markInterior) {
+                                                                                     bool closed, bool markInterior, bool removeOverlap) {
 
 
   std::vector<Halfedge> halfedges;
@@ -307,8 +307,33 @@ 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(!removeOverlap || halfedges.size() == 0) {
+  	return std::unique_ptr<FlipEdgeNetwork>(new FlipEdgeNetwork(mesh_, geom, {halfedges}, extraMark));
+	}
+
+	// 2 Adjacent Dijkstra paths may share edge(s) 
+	// This creates a 'back-and-forth' that FlipOut cannot solve
+	std::vector<Halfedge> heClean;
+	Halfedge hePrev, heNext;
+	size_t i, behind;
+	behind = 1;
+	heClean.push_back(halfedges[0]);
+	for(i = 1; i < halfedges.size(); i++) {
+		hePrev = halfedges[i - behind];
+		heNext = halfedges[i];
+		heClean.push_back(halfedges[i]);
+		if(hePrev.edge() == heNext.twin().edge()) {
+			/* Backtrack in case >2 edges overlap */
+			behind = (behind >= i) ? 1 : behind + 2; 
+			heClean.pop_back();
+			heClean.pop_back();
+		} else {
+			behind = 1;
+		}
+	}
+
+  return std::unique_ptr<FlipEdgeNetwork>(new FlipEdgeNetwork(mesh_, geom, {heClean}, extraMark));
 }
 
 

test/CMakeLists.txt

@@ -109,6 +109,7 @@ set(TEST_SRCS
   src/polygon_operators_test.cpp
   src/stl_reader_test.cpp
   src/surface_misc_test.cpp
+  src/piecewise_edge_network_test.cpp
 )
 
 add_executable(geometry-central-test "${TEST_SRCS}")

test/assets/dijkstra.obj

@@ -0,0 +1,12 @@
+v 0.000000 0.000000 0.000000
+v 0.000000 3.000000 0.000000
+v 0.000000 3.000000 1.000000
+v 0.000000 5.000000 0.000000
+v 0.000000 5.000000 5.000000
+v 0.000000 3.000000 2.000000
+v 0.000000 5.000000 0.000000
+f 1 2 7
+f 2 4 5
+f 3 2 5
+f 3 5 6
+f 2 3 7

test/src/load_test_meshes.cpp

@@ -87,6 +87,7 @@ void MeshAssetSuite::SetUpTestSuite() {
   allMeshAssets.emplace_back("platonic_shelf.obj", true);
   allMeshAssets.emplace_back("fox.ply", true);
   allMeshAssets.emplace_back("cat_head.obj", true);
+  allMeshAssets.emplace_back("dijkstra.obj", true);
 
   // Load general surface mesh variants
   allMeshAssets.emplace_back("tet.obj", false);

test/src/piecewise_edge_network_test.cpp

@@ -0,0 +1,93 @@
+#include "geometrycentral/surface/surface_mesh.h"
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+#include <geometrycentral/surface/vertex_position_geometry.h>
+#include <geometrycentral/surface/meshio.h>
+#include <geometrycentral/utilities/vector3.h>
+#include "geometrycentral/surface/flip_geodesics.h"
+
+#include <Eigen/Dense>
+#include <Eigen/src/Core/Matrix.h>
+
+#include "load_test_meshes.h"
+
+#include "gtest/gtest.h"
+
+#include <string>
+#include <iostream>
+#include <memory>
+
+using namespace geometrycentral;
+using namespace geometrycentral::surface;
+using std::cout;
+using std::endl;
+
+class pieceWiseEdgeNetworkTest : public MeshAssetSuite {};
+
+TEST_F(pieceWiseEdgeNetworkTest, deletesSingleDetour) {
+	for(auto& asset : {getAsset("dijkstra.obj", true)}) {
+		Eigen::MatrixXd out;
+		ManifoldSurfaceMesh& mesh = *asset.manifoldMesh; 
+		std::unique_ptr<geometrycentral::surface::FlipEdgeNetwork> network;
+    VertexPositionGeometry& origGeometry = *asset.geometry;
+		std::vector<Vertex> path;
+		std::vector<int> points = {0, 2, 3};
+
+		for(size_t i = 0; i < points.size(); i++) {
+			path.push_back(mesh.vertex(points[i]));
+		}
+
+		network = FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, origGeometry, path, false);
+		for(auto &eptr : network->paths){
+			auto x = eptr->getHalfedgeList();
+			for(auto he : x) {
+				std::cout << he.vertex() << std::endl;
+			}
+			EXPECT_EQ(x.size(), 2);
+		}
+
+	}
+
+}
+
+TEST_F(pieceWiseEdgeNetworkTest, deleteMultipleAtStart) {
+	for(auto& asset : {getAsset("dijkstra.obj", true)}) {
+		Eigen::MatrixXd out;
+		ManifoldSurfaceMesh& mesh = *asset.manifoldMesh; 
+		std::unique_ptr<geometrycentral::surface::FlipEdgeNetwork> network;
+    VertexPositionGeometry& origGeometry = *asset.geometry;
+		std::vector<Vertex> path;
+		std::vector<int> points = {1, 5, 3};
+
+		for(size_t i = 0; i < points.size(); i++) {
+			path.push_back(mesh.vertex(points[i]));
+		}
+
+		network = FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, origGeometry, path, false);
+		for(auto &eptr : network->paths){
+			EXPECT_EQ(eptr->getHalfedgeList().size(), 1);
+		}
+
+	}
+
+}
+
+TEST_F(pieceWiseEdgeNetworkTest, deleteMultipleAtEnd) {
+	for(auto& asset : {getAsset("dijkstra.obj", true)}) {
+		Eigen::MatrixXd out;
+		ManifoldSurfaceMesh& mesh = *asset.manifoldMesh; 
+		std::unique_ptr<geometrycentral::surface::FlipEdgeNetwork> network;
+    VertexPositionGeometry& origGeometry = *asset.geometry;
+		std::vector<Vertex> path;
+		std::vector<int> points = {0, 2, 5, 1};
+
+		for(size_t i = 0; i < points.size(); i++) {
+			path.push_back(mesh.vertex(points[i]));
+		}
+
+		network = FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, origGeometry, path, false);
+		for(auto &eptr : network->paths){
+			EXPECT_EQ(eptr->getHalfedgeList().size(), 1);
+		}
+	}
+
+}