diff --git a/src/surface/common_subdivision.cpp b/src/surface/common_subdivision.cpp index 0cf727c6..d4961f55 100644 --- a/src/surface/common_subdivision.cpp +++ b/src/surface/common_subdivision.cpp @@ -193,7 +193,12 @@ EdgeData CommonSubdivision::interpolateEdgeLengthsA(const EdgeDataedges()) { SurfacePoint tail = sourcePoints[e.halfedge().tailVertex()]->posA; SurfacePoint tip = sourcePoints[e.halfedge().tipVertex()]->posA; - Face f = sharedFace(tail, tip); + // Resolve the containing source face from an adjacent common subdivision + // face rather than from the endpoint SurfacePoints: on a Delta-complex, + // the endpoints may share several elements (e.g. two vertices joined by + // multiple edges), and sharedFace() could pick a face adjacent to the + // wrong one, producing a wildly wrong length. + Face f = sourceFaceA[e.halfedge().face()]; GC_SAFETY_ASSERT(f != Face(), "common subdivision edges must be contained in a face"); tail = tail.inFace(f); tip = tip.inFace(f); @@ -216,7 +221,10 @@ EdgeData CommonSubdivision::interpolateEdgeLengthsB(const EdgeDataedges()) { SurfacePoint tail = sourcePoints[e.halfedge().tailVertex()]->posB; SurfacePoint tip = sourcePoints[e.halfedge().tipVertex()]->posB; - Face f = sharedFace(tail, tip); + // See the comment in interpolateEdgeLengthsA: resolve the containing + // source face from an adjacent common subdivision face, since the + // endpoint SurfacePoints alone are ambiguous on a Delta-complex. + Face f = sourceFaceB[e.halfedge().face()]; GC_SAFETY_ASSERT(f != Face(), "common subdivision edges must be contained in a face"); tail = tail.inFace(f); tip = tip.inFace(f); diff --git a/test/src/intrinsic_triangulation_test.cpp b/test/src/intrinsic_triangulation_test.cpp index 01128fdc..4672bf05 100644 --- a/test/src/intrinsic_triangulation_test.cpp +++ b/test/src/intrinsic_triangulation_test.cpp @@ -370,6 +370,41 @@ TEST_F(IntrinsicTriangulationSuite, CommonSubdivisionCompareIntegerSignpost) { // TODO: also test with signposts? +// A single edge flip on a tetrahedron creates a Delta-complex in which two +// vertices are joined by two distinct edges. interpolateEdgeLengthsA/B used +// to resolve each common subdivision edge's containing face from its +// endpoint SurfacePoints alone, which is ambiguous in this configuration +// and produced lengths of the wrong edge. +TEST_F(IntrinsicTriangulationSuite, CommonSubdivisionEdgeLengthsDeltaComplex) { + for (size_t iE = 0; iE < 6; iE++) { + auto a = getAsset("tet.obj", true); + ManifoldSurfaceMesh& mesh = *a.manifoldMesh; + VertexPositionGeometry& origGeometry = *a.geometry; + + IntegerCoordinatesIntrinsicTriangulation tri(mesh, origGeometry); + + if (!tri.flipEdgeIfPossible(tri.intrinsicMesh->edge(iE))) continue; + + CommonSubdivision& cs = tri.getCommonSubdivision(); + cs.constructMesh(); + + // Ground truth: every common subdivision edge lies in a single face of + // both meshes, so extrinsic positions interpolated across mesh A give + // exact lengths + VertexPositionGeometry csGeo(*cs.mesh, cs.interpolateAcrossA(origGeometry.vertexPositions)); + csGeo.requireEdgeLengths(); + + EdgeData lengthsFromLenB = cs.interpolateEdgeLengthsB(tri.edgeLengths); + origGeometry.requireEdgeLengths(); + EdgeData lengthsFromLenA = cs.interpolateEdgeLengthsA(origGeometry.edgeLengths); + + for (Edge e : cs.mesh->edges()) { + EXPECT_NEAR(csGeo.edgeLengths[e], lengthsFromLenB[e], 1e-5); + EXPECT_NEAR(csGeo.edgeLengths[e], lengthsFromLenA[e], 1e-5); + } + } +} + TEST_F(IntrinsicTriangulationSuite, FunctionTransfer) { for (const MeshAsset& a : {getAsset("fox.ply", true)}) { a.printThyName();