diff --git a/include/geometrycentral/surface/signpost_intrinsic_triangulation.h b/include/geometrycentral/surface/signpost_intrinsic_triangulation.h index 2653a1e7..8bdd9e58 100644 --- a/include/geometrycentral/surface/signpost_intrinsic_triangulation.h +++ b/include/geometrycentral/surface/signpost_intrinsic_triangulation.h @@ -64,10 +64,10 @@ class SignpostIntrinsicTriangulation : public IntrinsicGeometryInterface { // ====================================================== // Given a point on the input triangulation, returns the corresponding point on the intrinsic triangulation - SurfacePoint equivalentPointOnIntrinsic(const SurfacePoint& pointOnInput); + SurfacePoint equivalentPointOnIntrinsic(SurfacePoint pointOnInput); // Given a point on the intrinsic triangulation, returns the corresponding point on the input triangulation - SurfacePoint equivalentPointOnInput(const SurfacePoint& pointOnIntrinsic); + SurfacePoint equivalentPointOnInput(SurfacePoint pointOnIntrinsic); // Trace out the edges of the intrinsic triangulation along the surface of the input mesh. // Each path is ordered along edge.halfedge(), and includes both the start and end points diff --git a/src/surface/signpost_intrinsic_triangulation.cpp b/src/surface/signpost_intrinsic_triangulation.cpp index 3e41fc0f..22eebbf2 100644 --- a/src/surface/signpost_intrinsic_triangulation.cpp +++ b/src/surface/signpost_intrinsic_triangulation.cpp @@ -32,6 +32,7 @@ SignpostIntrinsicTriangulation::SignpostIntrinsicTriangulation(ManifoldSurfaceMe // == Initialize geometric data inputGeom.requireEdgeLengths(); inputGeom.requireHalfedgeVectorsInVertex(); + inputGeom.requireHalfedgeVectorsInFace(); inputGeom.requireVertexAngleSums(); inputGeom.requireCornerAngles(); @@ -91,6 +92,115 @@ void SignpostIntrinsicTriangulation::setMarkedEdges(const EdgeData& marked markedEdges.setDefault(false); } +SurfacePoint SignpostIntrinsicTriangulation::equivalentPointOnIntrinsic(SurfacePoint pointOnInput) { + pointOnInput = pointOnInput.reduced(); + + // Vertex on inputMesh is preserved on intrinsicMesh + if (pointOnInput.type == SurfacePointType::Vertex) { + return SurfacePoint(intrinsicMesh->vertex(pointOnInput.vertex.getIndex())); + } + + // If edge on inputMesh is preserved, simply return it. Otherwise treat it as a face point. + if (pointOnInput.type == SurfacePointType::Edge) { + if (edgeIsOriginal[pointOnInput.edge.getIndex()]) { + return SurfacePoint(intrinsicMesh->edge(pointOnInput.edge.getIndex()), pointOnInput.tEdge); + } + pointOnInput = pointOnInput.inSomeFace(); + } + + // Examine all possible tracings from the three vertices of pointOnInput.face + std::array startP; + std::array traceVecAngle; + std::array traceVecLen; + for (int inputHE_offset = 0; inputHE_offset < 3; ++inputHE_offset) { + Halfedge inputHE = pointOnInput.face.halfedge(); + for (int i = 0; i < inputHE_offset; ++i) { + inputHE = inputHE.next(); + } + + Vertex inputV = inputHE.vertex(); + startP[inputHE_offset] = intrinsicMesh->vertex(inputV.getIndex()); + + // Get tracing vector from inputV in a local coordinate frame defined by inputHE + Vector2 inputHE_unitVecInFace = inputGeom.halfedgeVectorsInFace[inputHE].normalize(); + std::array vertCoords = {{ + {0., 0.}, + inputGeom.halfedgeVectorsInFace[inputHE] / inputHE_unitVecInFace, // Rotate halfedgeVectorsInFace such that inputHE becomes the X axis + -inputGeom.halfedgeVectorsInFace[inputHE.next().next()] / inputHE_unitVecInFace + }}; + Vector2 traceVec_local = pointOnInput.faceCoords[(1 + inputHE_offset) % 3] * vertCoords[1] + pointOnInput.faceCoords[(2 + inputHE_offset) % 3] * vertCoords[2]; + + // Adjust tracing angle by rescaling and offsetting + traceVecAngle[inputHE_offset] = traceVec_local.arg(); + traceVecAngle[inputHE_offset] *= (inputV.isBoundary() ? 1. : 2.) * M_PI / inputGeom.vertexAngleSums[inputV]; + traceVecAngle[inputHE_offset] += inputGeom.halfedgeVectorsInVertex[inputHE].arg(); + traceVecLen[inputHE_offset] = traceVec_local.norm(); + } + + // Select traceVec whose length is the smallest + int selected = + traceVecLen[0] < traceVecLen[1] && traceVecLen[0] < traceVecLen[2] ? 0 : + traceVecLen[1] < traceVecLen[2] ? 1 : + 2; + Vector2 traceVec = Vector2::fromAngle(traceVecAngle[selected]) * traceVecLen[selected]; + TraceGeodesicResult intrinsicTraceResult = traceGeodesic(*this, startP[selected], traceVec); + return intrinsicTraceResult.endPoint; +} + +SurfacePoint SignpostIntrinsicTriangulation::equivalentPointOnInput(SurfacePoint pointOnIntrinsic) { + pointOnIntrinsic = pointOnIntrinsic.reduced(); + + // We already know where each intrinsicMesh vertex is located on inputMesh + if (pointOnIntrinsic.type == SurfacePointType::Vertex) { + return vertexLocations[pointOnIntrinsic.vertex]; + } + + // If intrinsicMesh edge is preserved, simply return it. Otherwise treat it as a face point. + if (pointOnIntrinsic.type == SurfacePointType::Edge) { + if (edgeIsOriginal[pointOnIntrinsic.edge]) { + return SurfacePoint(inputMesh.edge(pointOnIntrinsic.edge.getIndex()), pointOnIntrinsic.tEdge); + } + pointOnIntrinsic = pointOnIntrinsic.inSomeFace(); + } + + // Examine all possible tracings from the three vertices of pointOnIntrinsic.face + std::array startP; + std::array traceVecAngle; + std::array traceVecLen; + for (int intrinsicHE_offset = 0; intrinsicHE_offset < 3; ++intrinsicHE_offset) { + Halfedge intrinsicHE = pointOnIntrinsic.face.halfedge(); + for (int i = 0; i < intrinsicHE_offset; ++i) { + intrinsicHE = intrinsicHE.next(); + } + + Vertex intrinsicV = intrinsicHE.vertex(); + startP[intrinsicHE_offset] = vertexLocations[intrinsicV]; + + // Get tracing vector from intrinsicV in a local coordinate frame defined by intrinsicHE + Vector2 intrinsicHE_unitVecInFace = halfedgeVectorsInFace[intrinsicHE].normalize(); + std::array vertCoords = {{ + {0., 0.}, + halfedgeVectorsInFace[intrinsicHE] / intrinsicHE_unitVecInFace, // Rotate halfedgeVectorsInFace such that intrinsicHE becomes the X axis + -halfedgeVectorsInFace[intrinsicHE.next().next()] / intrinsicHE_unitVecInFace + }}; + Vector2 traceVec_local = pointOnIntrinsic.faceCoords[(1 + intrinsicHE_offset) % 3] * vertCoords[1] + pointOnIntrinsic.faceCoords[(2 + intrinsicHE_offset) % 3] * vertCoords[2]; + + // Adjust tracing angle by rescaling and offsetting + traceVecAngle[intrinsicHE_offset] = traceVec_local.arg(); + traceVecAngle[intrinsicHE_offset] *= 1.0 / vertexAngleScaling(intrinsicV); + traceVecAngle[intrinsicHE_offset] += halfedgeVector(intrinsicHE).arg(); + traceVecLen[intrinsicHE_offset] = traceVec_local.norm(); + } + + // Select traceVec whose length is the smallest + int selected = + traceVecLen[0] < traceVecLen[1] && traceVecLen[0] < traceVecLen[2] ? 0 : + traceVecLen[1] < traceVecLen[2] ? 1 : + 2; + Vector2 traceVec = Vector2::fromAngle(traceVecAngle[selected]) * traceVecLen[selected]; + TraceGeodesicResult inputTraceResult = traceGeodesic(inputGeom, startP[selected], traceVec); + return inputTraceResult.endPoint; +} std::vector SignpostIntrinsicTriangulation::traceHalfedge(Halfedge he, bool trimEnd) { diff --git a/src/surface/trace_geodesic.cpp b/src/surface/trace_geodesic.cpp index 2cbaa905..6cf905e2 100644 --- a/src/surface/trace_geodesic.cpp +++ b/src/surface/trace_geodesic.cpp @@ -699,6 +699,7 @@ TraceGeodesicResult traceGeodesic(IntrinsicGeometryInterface& geom, SurfacePoint geom.unrequireHalfedgeVectorsInVertex(); geom.unrequireHalfedgeVectorsInFace(); + result.endPoint = startP; result.endingDir = Vector2::zero(); // probably want to ensure we still return a point in a face... diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e8966344..ea2b8865 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -98,6 +98,7 @@ set(TEST_SRCS src/halfedge_geometry_test.cpp src/surface_misc_test.cpp src/linear_algebra_test.cpp + src/signpost_intrinsic_triangulation_test.cpp ) add_executable(geometry-central-test "${TEST_SRCS}") diff --git a/test/src/signpost_intrinsic_triangulation_test.cpp b/test/src/signpost_intrinsic_triangulation_test.cpp new file mode 100644 index 00000000..e64cf5f2 --- /dev/null +++ b/test/src/signpost_intrinsic_triangulation_test.cpp @@ -0,0 +1,194 @@ +#include "geometrycentral/surface/signpost_intrinsic_triangulation.h" + +#include "load_test_meshes.h" + +using namespace geometrycentral; +using namespace geometrycentral::surface; + +// helpers +namespace { + +std::mt19937 mt(42); + +double unitRand() { + std::uniform_real_distribution dist(0.0, 1.0); + return dist(mt); +} + +void insertVerticesAtRandom(SignpostIntrinsicTriangulation& signpostTri) { + for (Face f : signpostTri.intrinsicMesh->faces()) { + if (unitRand() > 0.5) { + signpostTri.insertVertex(SurfacePoint{f, {1/3., 1/3., 1/3.}}); + } + } +} + +void flipEdgesAtRandom(SignpostIntrinsicTriangulation& signpostTri) { + for (Edge e : signpostTri.intrinsicMesh->edges()) { + if (unitRand() > 0.5) { + signpostTri.flipEdgeIfPossible(e); + } + } +} + +void mutateIntrinsicMesh(SignpostIntrinsicTriangulation& signpostTri) { + insertVerticesAtRandom(signpostTri); + flipEdgesAtRandom(signpostTri); +} + +double getBoundingBoxDiagonal(VertexPositionGeometry& geometry) { + Vector3 minP = Vector3::constant(std::numeric_limits::infinity()); + Vector3 maxP = Vector3::constant(-std::numeric_limits::infinity()); + for (Vertex v : geometry.mesh.vertices()) { + Vector3 p = geometry.vertexPositions[v]; + minP = componentwiseMin(minP, p); + maxP = componentwiseMax(maxP, p); + } + return norm(minP - maxP); +} + +const int N = 7; +const double EPS = 1e-3; + +} // namespace + +class SignpostIntrinsicTriangulationSuite : public MeshAssetSuite {}; + +// ============================================================ +// =============== Equivalent point tests +// ============================================================ + +TEST_F(SignpostIntrinsicTriangulationSuite, EquivalentPointOnIntrinsic_FacePoint) { + for (auto& asset : {getAsset("bob_small.ply", true), getAsset("lego.ply", true), getAsset("sphere_small.ply", true), getAsset("spot.ply", true)}) { + asset.printThyName(); + SignpostIntrinsicTriangulation signpostTri(*asset.manifoldMesh, *asset.geometry); + + mutateIntrinsicMesh(signpostTri); + + double boundingBoxDiagonal = getBoundingBoxDiagonal(*asset.geometry); + + for (Face f : signpostTri.inputMesh.faces()) { + Vector3 i; + for (i[0] = 1; i[0] < N; ++i[0]) { + for (i[1] = 1; i[1] < N - i[0]; ++i[1]) { + i[2] = N - i[0] - i[1]; + Vector3 faceCoords = i / N; + SurfacePoint pointOnInput_before{f, faceCoords}; + SurfacePoint pointOnIntrinsic = signpostTri.equivalentPointOnIntrinsic(pointOnInput_before); + SurfacePoint pointOnInput_after = signpostTri.equivalentPointOnInput(pointOnIntrinsic); + double error = (pointOnInput_before.interpolate(asset.geometry->inputVertexPositions) - pointOnInput_after.interpolate(asset.geometry->inputVertexPositions)).norm(); + EXPECT_LT(error, boundingBoxDiagonal * EPS); + } + } + } + } +} + +TEST_F(SignpostIntrinsicTriangulationSuite, EquivalentPointOnInput_FacePoint) { + for (auto& asset : {getAsset("bob_small.ply", true), getAsset("lego.ply", true), getAsset("sphere_small.ply", true), getAsset("spot.ply", true)}) { + asset.printThyName(); + SignpostIntrinsicTriangulation signpostTri(*asset.manifoldMesh, *asset.geometry); + + mutateIntrinsicMesh(signpostTri); + + double boundingBoxDiagonal = getBoundingBoxDiagonal(*asset.geometry); + + for (Face f : signpostTri.intrinsicMesh->faces()) { + Vector3 i; + for (i[0] = 1; i[0] < N; ++i[0]) { + for (i[1] = 1; i[1] < N - i[0]; ++i[1]) { + i[2] = N - i[0] - i[1]; + Vector3 faceCoords = i / N; + SurfacePoint pointOnIntrinsic_before{f, faceCoords}; + SurfacePoint pointOnInput_before = signpostTri.equivalentPointOnInput(pointOnIntrinsic_before); + SurfacePoint pointOnIntrinsic_after = signpostTri.equivalentPointOnIntrinsic(pointOnInput_before); + SurfacePoint pointOnInput_after = signpostTri.equivalentPointOnInput(pointOnIntrinsic_after); + double error = (pointOnInput_before.interpolate(asset.geometry->inputVertexPositions) - pointOnInput_after.interpolate(asset.geometry->inputVertexPositions)).norm(); + EXPECT_LT(error, boundingBoxDiagonal * EPS); + } + } + } + } +} + +TEST_F(SignpostIntrinsicTriangulationSuite, EquivalentPointOnIntrinsic_EdgePoint) { + for (auto& asset : {getAsset("bob_small.ply", true), getAsset("lego.ply", true), getAsset("sphere_small.ply", true), getAsset("spot.ply", true)}) { + asset.printThyName(); + SignpostIntrinsicTriangulation signpostTri(*asset.manifoldMesh, *asset.geometry); + + mutateIntrinsicMesh(signpostTri); + + double boundingBoxDiagonal = getBoundingBoxDiagonal(*asset.geometry); + + for (Edge e : signpostTri.inputMesh.edges()) { + for (double i = 1; i < N; ++i) { + SurfacePoint pointOnInput_before{e, i / N}; + SurfacePoint pointOnIntrinsic = signpostTri.equivalentPointOnIntrinsic(pointOnInput_before); + SurfacePoint pointOnInput_after = signpostTri.equivalentPointOnInput(pointOnIntrinsic); + double error = (pointOnInput_before.interpolate(asset.geometry->inputVertexPositions) - pointOnInput_after.interpolate(asset.geometry->inputVertexPositions)).norm(); + EXPECT_LT(error, boundingBoxDiagonal * EPS); + } + } + } +} + +TEST_F(SignpostIntrinsicTriangulationSuite, EquivalentPointOnInput_EdgePoint) { + for (auto& asset : {getAsset("bob_small.ply", true), getAsset("lego.ply", true), getAsset("sphere_small.ply", true), getAsset("spot.ply", true)}) { + asset.printThyName(); + SignpostIntrinsicTriangulation signpostTri(*asset.manifoldMesh, *asset.geometry); + + mutateIntrinsicMesh(signpostTri); + + double boundingBoxDiagonal = getBoundingBoxDiagonal(*asset.geometry); + + for (Edge e : signpostTri.intrinsicMesh->edges()) { + for (double i = 1; i < N; ++i) { + SurfacePoint pointOnIntrinsic_before{e, i / N}; + SurfacePoint pointOnInput_before = signpostTri.equivalentPointOnInput(pointOnIntrinsic_before); + SurfacePoint pointOnIntrinsic_after = signpostTri.equivalentPointOnIntrinsic(pointOnInput_before); + SurfacePoint pointOnInput_after = signpostTri.equivalentPointOnInput(pointOnIntrinsic_after); + double error = (pointOnInput_before.interpolate(asset.geometry->inputVertexPositions) - pointOnInput_after.interpolate(asset.geometry->inputVertexPositions)).norm(); + EXPECT_LT(error, boundingBoxDiagonal * EPS); + } + } + } +} + +TEST_F(SignpostIntrinsicTriangulationSuite, EquivalentPointOnIntrinsic_VertexPoint) { + for (auto& asset : {getAsset("bob_small.ply", true), getAsset("lego.ply", true), getAsset("sphere_small.ply", true), getAsset("spot.ply", true)}) { + asset.printThyName(); + SignpostIntrinsicTriangulation signpostTri(*asset.manifoldMesh, *asset.geometry); + + mutateIntrinsicMesh(signpostTri); + + double boundingBoxDiagonal = getBoundingBoxDiagonal(*asset.geometry); + + for (Vertex v : signpostTri.inputMesh.vertices()) { + SurfacePoint pointOnInput_before{v}; + SurfacePoint pointOnIntrinsic = signpostTri.equivalentPointOnIntrinsic(pointOnInput_before); + SurfacePoint pointOnInput_after = signpostTri.equivalentPointOnInput(pointOnIntrinsic); + double error = (pointOnInput_before.interpolate(asset.geometry->inputVertexPositions) - pointOnInput_after.interpolate(asset.geometry->inputVertexPositions)).norm(); + EXPECT_LT(error, boundingBoxDiagonal * EPS); + } + } +} + +TEST_F(SignpostIntrinsicTriangulationSuite, EquivalentPointOnInput_VertexPoint) { + for (auto& asset : {getAsset("bob_small.ply", true), getAsset("lego.ply", true), getAsset("sphere_small.ply", true), getAsset("spot.ply", true)}) { + asset.printThyName(); + SignpostIntrinsicTriangulation signpostTri(*asset.manifoldMesh, *asset.geometry); + + mutateIntrinsicMesh(signpostTri); + + double boundingBoxDiagonal = getBoundingBoxDiagonal(*asset.geometry); + + for (Vertex v : signpostTri.intrinsicMesh->vertices()) { + SurfacePoint pointOnIntrinsic_before{v}; + SurfacePoint pointOnInput_before = signpostTri.equivalentPointOnInput(pointOnIntrinsic_before); + SurfacePoint pointOnIntrinsic_after = signpostTri.equivalentPointOnIntrinsic(pointOnInput_before); + SurfacePoint pointOnInput_after = signpostTri.equivalentPointOnInput(pointOnIntrinsic_after); + double error = (pointOnInput_before.interpolate(asset.geometry->inputVertexPositions) - pointOnInput_after.interpolate(asset.geometry->inputVertexPositions)).norm(); + EXPECT_LT(error, boundingBoxDiagonal * EPS); + } + } +}