signpost_intrinsic_triangulation: implement equivalentPointOn{Intrins…

…ic,Input}
nmwsharp · Sep 10, 2020 · c96d703 · c96d703
1 parent 3423ff4
commit c96d703
Show file tree

Hide file tree

Showing 5 changed files with 308 additions and 2 deletions.
diff --git 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
@@ -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<char>& 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<SurfacePoint, 3> startP;
+  std::array<double, 3> traceVecAngle;
+  std::array<double, 3> 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<Vector2, 3> 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<SurfacePoint, 3> startP;
+  std::array<double, 3> traceVecAngle;
+  std::array<double, 3> 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<Vector2, 3> 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<SurfacePoint> SignpostIntrinsicTriangulation::traceHalfedge(Halfedge he, bool trimEnd) {
 

diff --git 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
@@ -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
@@ -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<double> 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<double>::infinity());
+  Vector3 maxP = Vector3::constant(-std::numeric_limits<double>::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);
+    }
+  }
+}