diff --git a/include/geometrycentral/surface/flip_geodesics.h b/include/geometrycentral/surface/flip_geodesics.h index 04de0aeb..2f820d7c 100644 --- a/include/geometrycentral/surface/flip_geodesics.h +++ b/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 { diff --git a/src/surface/flip_geodesics.cpp b/src/surface/flip_geodesics.cpp index a9abbcbb..43c11655 100644 --- a/src/surface/flip_geodesics.cpp +++ b/src/surface/flip_geodesics.cpp @@ -308,7 +308,25 @@ std::unique_ptr FlipEdgeNetwork::constructFromPiecewiseDijkstra halfedges.insert(halfedges.end(), dijkstraPath.begin(), dijkstraPath.end()); } - return std::unique_ptr(new FlipEdgeNetwork(mesh_, geom, {halfedges}, extraMark)); + if (markInterior || halfedges.size() == 0) { + return std::unique_ptr(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 heClean; + for (Halfedge he : halfedges) { + if (!heClean.empty() && heClean.back() == he.twin()) { + heClean.pop_back(); + } else { + heClean.push_back(he); + } + } + + return std::unique_ptr(new FlipEdgeNetwork(mesh_, geom, {heClean}, extraMark)); } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b739fc4b..ddb63090 100644 --- a/test/CMakeLists.txt +++ b/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 ) diff --git a/test/src/flip_geodesic_test.cpp b/test/src/flip_geodesic_test.cpp new file mode 100644 index 00000000..93f19f70 --- /dev/null +++ b/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 + +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 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 waypoints = {mesh.vertex(0), mesh.vertex(nV / 3), mesh.vertex(2 * nV / 3)}; + + std::unique_ptr 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 waypoints = {mesh.vertex(0), mesh.vertex(nV / 3), mesh.vertex(2 * nV / 3)}; + + std::unique_ptr 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 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 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> buildOverlapTestMesh() { + std::vector verts = {{0, 0, 0}, {0, 3, 0}, {0, 3, 1}, {0, 5, 0}, {0, 5, 5}, {0, 3, 2}, {0, 5, 0}}; + std::vector> 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 meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildOverlapTestMesh(); + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexPositionGeometry& geom = *geomPtr; + + std::vector path = {mesh.vertex(0), mesh.vertex(2), mesh.vertex(3)}; + std::unique_ptr 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 meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildOverlapTestMesh(); + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexPositionGeometry& geom = *geomPtr; + + std::vector path = {mesh.vertex(1), mesh.vertex(5), mesh.vertex(3)}; + std::unique_ptr 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 meshPtr; + std::unique_ptr geomPtr; + std::tie(meshPtr, geomPtr) = buildOverlapTestMesh(); + ManifoldSurfaceMesh& mesh = *meshPtr; + VertexPositionGeometry& geom = *geomPtr; + + std::vector path = {mesh.vertex(0), mesh.vertex(2), mesh.vertex(5), mesh.vertex(1)}; + std::unique_ptr network = + FlipEdgeNetwork::constructFromPiecewiseDijkstraPath(mesh, geom, path, false); + + ASSERT_EQ(network->paths.size(), 1); + EXPECT_EQ(network->paths[0]->getHalfedgeList().size(), 1); +}