Add marching triangles

docs/docs/surface/utilities/barycentric_vector.md

@@ -25,7 +25,7 @@ which indicates what kind of vector it is.
 
 ??? func "`#!cpp BarycentricVector::BarycentricVector(Face f, Vector3 faceCoords)`"
 
-    Construct a barycentric vector that lies along the given face `f`, with the given barycentric coordinates `faceCoords`.
+    Construct a barycentric vector that lies in the given face `f`, with the given barycentric coordinates `faceCoords`.
 
 ??? func "`#!cpp BarycentricVector::BarycentricVector(Edge e, Vector2 edgeCoords)`"
 
@@ -35,6 +35,10 @@ which indicates what kind of vector it is.
 
     Construct a (zero-length) barycentric vector that lies on the given vertex `v`.
 
+??? func "`#!cpp BarycentricVector::BarycentricVector(Halfedge he, Face f)`"
+
+    Construct a barycentric vector from the given halfedge, that lies in the given face `f`.
+
 ??? func "`#!cpp BarycentricVector::BarycentricVector(SurfacePoint pA, SurfacePoint pB)`"
 
     Construct a barycentric vector from the given [surface points](../surface_point), i.e. barycentric points. The vector direction goes from `pA` to `pB`.

include/geometrycentral/surface/barycentric_vector.h

@@ -20,6 +20,7 @@ struct BarycentricVector {
   BarycentricVector(Vertex v);
   BarycentricVector(Edge e, Vector2 edgeCoords);
   BarycentricVector(Face f, Vector3 faceCoords);
+  BarycentricVector(Halfedge he, Face f);
 
   BarycentricVectorType type = BarycentricVectorType::Face;
 

include/geometrycentral/surface/barycentric_vector.ipp

@@ -44,6 +44,23 @@ inline BarycentricVector::BarycentricVector(SurfacePoint pA, SurfacePoint pB) {
   faceCoords = pB_in_face.faceCoords - pA_in_face.faceCoords;
 }
 
+inline BarycentricVector::BarycentricVector(Halfedge he_, Face f) : type(BarycentricVectorType::Face), face(f) {
+
+  int eIdx = 0;
+  double sign = 0.;
+  for (Halfedge he : f.adjacentHalfedges()) {
+    if (he.edge() == he_.edge()) {
+      sign = (he.tailVertex() == he_.tailVertex() && he.tipVertex() == he_.tipVertex()) ? 1. : -1.;
+      break;
+    }
+    eIdx++;
+  }
+  faceCoords = {0, 0, 0};
+  faceCoords[(eIdx + 1) % 3] = 1;
+  faceCoords[eIdx] = -1;
+  faceCoords *= sign;
+}
+
 // == Methods
 
 inline BarycentricVector BarycentricVector::inSomeFace() const {

include/geometrycentral/surface/marching_triangles.h

@@ -0,0 +1,12 @@
+#pragma once
+
+#include "geometrycentral/surface/barycentric_vector.h"
+
+namespace geometrycentral {
+namespace surface {
+
+std::vector<std::vector<SurfacePoint>> marchingTriangles(IntrinsicGeometryInterface& geom, const VertexData<double>& u,
+                                                         double isoval = 0.);
+
+}
+} // namespace geometrycentral

include/geometrycentral/surface/signed_heat_method.h

@@ -80,7 +80,6 @@ class SignedHeatSolver {
   double lengthOfSegment(const SurfacePoint& pA, const SurfacePoint& pB) const;
   SurfacePoint midSegmentSurfacePoint(const SurfacePoint& pA, const SurfacePoint& pB) const;
   std::complex<double> projectedNormal(const SurfacePoint& pA, const SurfacePoint& pB, const Edge& e) const;
-  BarycentricVector barycentricVectorInFace(const Halfedge& he, const Face& f) const;
   FaceData<BarycentricVector> sampleAtFaceBarycenters(const Vector<std::complex<double>>& Xt);
   Vector<double> integrateWithZeroSetConstraint(const Vector<double>& rhs, const std::vector<Curve>& curves,
                                                 const std::vector<SurfacePoint>& points,

include/geometrycentral/surface/surface_point.ipp

@@ -151,9 +151,9 @@ inline SurfacePoint SurfacePoint::inEdge(Edge targetEdge) const {
   switch (type) {
   case SurfacePointType::Vertex:
     if (vertex == targetEdge.halfedge().tailVertex()) {
-      return SurfacePoint(targetEdge, 0);
-    } else if (vertex == targetEdge.halfedge().tipVertex()) {
       return SurfacePoint(targetEdge, 1);
+    } else if (vertex == targetEdge.halfedge().tipVertex()) {
+      return SurfacePoint(targetEdge, 0);
     }
     break;
   case SurfacePointType::Edge:

src/CMakeLists.txt

@@ -26,6 +26,7 @@ SET(SRCS
   surface/embedded_geometry_interface.cpp
   surface/edge_length_geometry.cpp
   surface/vertex_position_geometry.cpp
+  surface/marching_triangles.cpp
   surface/mutation_manager.cpp
   surface/mesh_graph_algorithms.cpp
   surface/direction_fields.cpp

src/surface/fast_marching_method.cpp

@@ -92,17 +92,17 @@ VertexData<double> FMMDistance(IntrinsicGeometryInterface& geometry,
           // Assign +/- signs to the "third" vertices of each face straddling this edge.
           // These vertices might themselves lie on the curve, in which case we overwrite them below.
           Halfedge he = commonEdge.halfedge();
-          signs[he.next().tipVertex()] = (he.vertex() == pA.vertex) ? 1 : -1;
+          signs[he.next().tipVertex()] = (he.vertex() == pA.vertex) ? -1 : 1;
           signs[he.twin().next().tipVertex()] = -signs[he.next().tipVertex()];
         } else {
           Face commonFace = sharedFace(pA, pB);
           if (commonFace == Face()) {
             throw std::logic_error("For signed fast marching distance, each curve segment must share a common face.");
           }
           BarycentricVector tangent(pA, pB);
-          BarycentricVector normal = tangent.rotate90(geometry);
+          BarycentricVector normal = -tangent.rotate90(geometry);
           for (Vertex v : commonFace.adjacentVertices()) {
-            BarycentricVector u(SurfacePoint(v), pA);
+            BarycentricVector u(pA, SurfacePoint(v));
             signs[v] = (dot(geometry, normal, u) > 0) ? 1 : -1;
           }
         }

src/surface/marching_triangles.cpp

@@ -0,0 +1,129 @@
+#include "geometrycentral/surface/marching_triangles.h"
+
+namespace geometrycentral {
+namespace surface {
+
+namespace {
+
+std::vector<std::vector<std::array<size_t, 2>>>
+getCurveComponents(SurfaceMesh& mesh, const std::vector<SurfacePoint>& curveNodes,
+                   const std::vector<std::array<size_t, 2>>& curveEdges) {
+
+  // Note: This function assumes that `curveNodes` has been de-duplicated.
+  std::vector<std::array<size_t, 2>> edgesToAdd = curveEdges;
+  std::vector<std::vector<std::array<size_t, 2>>> curves;
+  size_t nSegs = curveEdges.size();
+  while (edgesToAdd.size() > 0) {
+    std::array<size_t, 2> startSeg = edgesToAdd.back();
+    edgesToAdd.pop_back();
+    curves.emplace_back();
+    std::vector<std::array<size_t, 2>>& currCurve = curves.back();
+    currCurve.push_back(startSeg);
+
+    // Add segs to the front end until we can't.
+    std::array<size_t, 2> currSeg = startSeg;
+    while (true) {
+      const SurfacePoint& front = curveNodes[currSeg[1]];
+      bool didWeFindOne = false;
+      for (size_t i = 0; i < edgesToAdd.size(); i++) {
+        std::array<size_t, 2> otherSeg = edgesToAdd[i];
+        if (curveNodes[otherSeg[0]] == front) {
+          currSeg = otherSeg;
+          currCurve.push_back(otherSeg);
+          edgesToAdd.erase(edgesToAdd.begin() + i);
+          didWeFindOne = true;
+          break;
+        }
+      }
+      if (!didWeFindOne) break;
+    }
+
+    // Add segs to the back end until we can't.
+    currSeg = startSeg;
+    while (true) {
+      const SurfacePoint& back = curveNodes[currSeg[0]];
+      bool didWeFindOne = false;
+      for (size_t i = 0; i < edgesToAdd.size(); i++) {
+        std::array<size_t, 2> otherSeg = edgesToAdd[i];
+        if (curveNodes[otherSeg[1]] == back) {
+          currSeg = otherSeg;
+          currCurve.insert(currCurve.begin(), otherSeg);
+          edgesToAdd.erase(edgesToAdd.begin() + i);
+          didWeFindOne = true;
+          break;
+        }
+      }
+      if (!didWeFindOne) break;
+    }
+  }
+  return curves;
+}
+} // namespace
+
+std::vector<std::vector<SurfacePoint>> marchingTriangles(IntrinsicGeometryInterface& geom, const VertexData<double>& u,
+                                                         double isoval) {
+
+  SurfaceMesh& mesh = *u.getMesh();
+  std::vector<SurfacePoint> nodes;
+  std::vector<std::array<size_t, 2>> edges;
+  for (Face f : mesh.faces()) {
+    std::vector<SurfacePoint> hits;
+    BarycentricVector gradient(f);
+    for (Halfedge he : f.adjacentHalfedges()) {
+      // Record edge crossings
+      Edge e = he.edge();
+      Vertex v0 = e.firstVertex();
+      Vertex v1 = e.secondVertex();
+      double u0 = u[v0];
+      double u1 = u[v1];
+      double lB = std::min(u0, u1);
+      double uB = std::max(u0, u1);
+      if (lB == uB && lB == isoval) {
+        hits.clear();
+        hits = {SurfacePoint(v0), SurfacePoint(v1)};
+        break;
+      }
+      double t = (isoval - lB) / (uB - lB);
+      if (u0 > u1) t = 1. - t;
+      if (t <= 1. && t >= 0.) hits.emplace_back(e, t);
+      // Compute gradient of the scalar function.
+      BarycentricVector heVec(he.next(), f);
+      BarycentricVector ePerp = heVec.rotate90(geom);
+      gradient += ePerp * u[he.vertex()];
+    }
+    if (hits.size() != 2) continue;
+
+    // Orient segments so that smaller values are always on the "inside" of the curve.
+    std::array<size_t, 2> seg;
+    for (int i = 0; i < 2; i++) {
+      auto iter = std::find(nodes.begin(), nodes.end(), hits[i]);
+      if (iter != nodes.end()) {
+        seg[i] = iter - nodes.begin();
+      } else {
+        nodes.push_back(hits[i]);
+        seg[i] = nodes.size() - 1;
+      }
+    }
+    BarycentricVector segTangent(hits[0], hits[1]);
+    BarycentricVector segNormal = -segTangent.inFace(f).rotate90(geom);
+    if (dot(geom, segNormal, gradient) > 0.) {
+      edges.push_back(seg);
+    } else {
+      edges.push_back({seg[1], seg[0]});
+    }
+  }
+  std::vector<std::vector<std::array<size_t, 2>>> components = getCurveComponents(mesh, nodes, edges);
+  size_t nCurves = components.size();
+  std::vector<std::vector<SurfacePoint>> curves(nCurves);
+  for (size_t i = 0; i < nCurves; i++) {
+    size_t nEdges = components[i].size();
+    for (size_t j = 0; j < nEdges; j++) {
+      curves[i].push_back(nodes[components[i][j][0]]);
+    }
+    curves[i].push_back(nodes[components[i][nEdges - 1][1]]);
+  }
+  return curves;
+}
+
+} // namespace surface
+} // namespace geometrycentral

src/surface/signed_heat_method.cpp

@@ -113,7 +113,7 @@ Vector<std::complex<double>> SignedHeatSolver::integrateVectorHeatFlow(const std
           for (Edge e : f.adjacentEdges()) {
             size_t eIdx = geom.edgeIndices[e];
             double w = scalarCrouzeixRaviart(b, e);
-            BarycentricVector heVec = barycentricVectorInFace(e.halfedge(), f);
+            BarycentricVector heVec(e.halfedge(), f);
             heVec /= heVec.norm(geom);
             BarycentricVector heVecN = heVec.rotate90(geom);
             triplets.emplace_back(m, eIdx, w * dot(geom, heVec, segTangent));
@@ -162,8 +162,8 @@ VertexData<double> SignedHeatSolver::integrateVectorField(const Vector<std::comp
       Halfedge heA = c.halfedge();
       Halfedge heB = heA.next().next();
       BarycentricVector Yj = Y[f];
-      BarycentricVector e1 = barycentricVectorInFace(heA, f);
-      BarycentricVector e2 = -barycentricVectorInFace(heB, f);
+      BarycentricVector e1(heA, f);
+      BarycentricVector e2 = -BarycentricVector(heB, f);
       double cotTheta1 = geom.halfedgeCotanWeights[heA];
       double cotTheta2 = geom.halfedgeCotanWeights[heB];
       double w1 = cotTheta1 * dot(geom, e1, Yj);
@@ -252,11 +252,11 @@ void SignedHeatSolver::buildUnsignedCurveSource(const Curve& curve, Vector<std::
       // Ordinarily, "double-sided" vector information would cancel out along the edge. So in each face, "smear"
       // the info out a bit by parallel-transporting the initial vectors to other edges in adjacent faces.
       Face f = he.face();
-      BarycentricVector segment = barycentricVectorInFace(he, f);
+      BarycentricVector segment(he, f);
       BarycentricVector segNormal = segment.rotate90(geom);
       for (Edge e : f.adjacentEdges()) {
         if (e == commonEdge) continue;
-        BarycentricVector edgeVec = barycentricVectorInFace(e.halfedge(), f);
+        BarycentricVector edgeVec(e.halfedge(), f);
         edgeVec /= edgeVec.norm(geom);
         double sinTheta = dot(geom, segment, edgeVec);
         double cosTheta = dot(geom, segNormal, edgeVec);
@@ -698,23 +698,6 @@ std::complex<double> SignedHeatSolver::projectedNormal(const SurfacePoint& pA, c
   return normal;
 }
 
-BarycentricVector SignedHeatSolver::barycentricVectorInFace(const Halfedge& he_, const Face& f) const {
-
-  int eIdx = 0;
-  double sign = 0.;
-  for (Halfedge he : f.adjacentHalfedges()) {
-    if (he.edge() == he_.edge()) {
-      sign = (he.tailVertex() == he_.tailVertex() && he.tipVertex() == he_.tipVertex()) ? 1. : -1.;
-      break;
-    }
-    eIdx++;
-  }
-  Vector3 faceCoords = {0, 0, 0};
-  faceCoords[(eIdx + 1) % 3] = 1;
-  faceCoords[eIdx] = -1;
-  return BarycentricVector(f, sign * faceCoords);
-}
-
 FaceData<BarycentricVector> SignedHeatSolver::sampleAtFaceBarycenters(const Vector<std::complex<double>>& Xt) {
 
   geom.requireEdgeIndices();
@@ -724,7 +707,7 @@ FaceData<BarycentricVector> SignedHeatSolver::sampleAtFaceBarycenters(const Vect
     Vector3 faceCoords = {0, 0, 0};
     for (Halfedge he : f.adjacentHalfedges()) {
       size_t eIdx = geom.edgeIndices[he.edge()];
-      BarycentricVector e1 = barycentricVectorInFace(he, f);
+      BarycentricVector e1(he, f);
       if (!he.orientation()) e1 *= -1;
       BarycentricVector e2 = e1.rotate90(geom);
       e1 /= e1.norm(geom);