Preprocess curves for the signed heat method (#229)

* First commit for preprocessing curves

* Preprocess curves by partitioning them into well-sampled components

Author: nzfeng

docs/docs/surface/algorithms/signed_heat_method.md

@@ -64,6 +64,8 @@ The diffusion time can also be changed using the following function.
 
     Re-computes the time used for short time heat flow `tCoef * h^2`, where `h` is the mean distance between nodes of the mesh.
 
+Note: On triangle meshes, each mesh edge crossed by an input curve should be explicitly represented by a SurfacePoint (otherwise the curve's geometry is not fully specified). If this is not the case, the code internally partitions input curves so that this condition is satisfied; the robustness of the heat method should fill in small gaps. Depending on how sparsely your curve is sampled, this may produce curve components with fewer than two nodes, which will be ignored; if many of your curves have fewer than two nodes, consider [FlipOut geodesics](../flip_geodesics) for shortest-path completion of curves.
+
 ## Helper Types
 
 ### Curves

include/geometrycentral/surface/signed_heat_method.h

@@ -61,6 +61,7 @@ class SignedHeatSolver {
   SparseMatrix<double> massMat, doubleMassMat, doubleConnectionLaplacian, doubleVectorOp;
 
   // Helpers
+  std::vector<Curve> preprocessCurves(const std::vector<Curve>& curves) const;
   Vector<std::complex<double>> integrateVectorHeatFlow(const std::vector<Curve>& curves,
                                                        const std::vector<SurfacePoint>& points,
                                                        const SignedHeatOptions& options);

src/surface/signed_heat_method.cpp

@@ -32,8 +32,9 @@ VertexData<double> SignedHeatSolver::computeDistance(const std::vector<Curve>& c
                                                      const std::vector<SurfacePoint>& points,
                                                      const SignedHeatOptions& options) {
 
-  Vector<std::complex<double>> Xt = integrateVectorHeatFlow(curves, points, options);
-  return integrateVectorField(Xt, curves, points, options);
+  std::vector<Curve> processedCurves = preprocessCurves(curves);
+  Vector<std::complex<double>> Xt = integrateVectorHeatFlow(processedCurves, points, options);
+  return integrateVectorField(Xt, processedCurves, points, options);
 }
 
 VertexData<double> SignedHeatSolver::computeDistance(const std::vector<Curve>& curves,
@@ -56,6 +57,34 @@ void SignedHeatSolver::setDiffusionTimeCoefficient(double tCoef_) {
   doubleVectorOp = doubleMassMat + shortTime * doubleConnectionLaplacian;
 }
 
+std::vector<Curve> SignedHeatSolver::preprocessCurves(const std::vector<Curve>& curves) const {
+  // If there are gaps in the curves (i.e. curves are not sampled densely along the mesh), then there is an ambiguity of
+  // how to connect up the points into curves.
+  // Rather than, for example, automatically connecting up curves using a geodesic path -- which would impose curves
+  // that the user might not have meant -- we simply break up the input curves into components that do reflect how
+  // they're sampled. (This may, however, create curves with fewer than 2 nodes, which would get ignored.)
+  std::vector<Curve> newCurves;
+  for (const Curve& curve : curves) {
+    newCurves.emplace_back();
+    newCurves.back().isSigned = curve.isSigned;
+    size_t nNodes = curve.nodes.size();
+    for (size_t i = 0; i < nNodes - 1; i++) {
+      const SurfacePoint& pA = curve.nodes[i];
+      const SurfacePoint& pB = curve.nodes[i + 1];
+      newCurves.back().nodes.push_back(pA);
+      Face commonFace = sharedFace(pA, pB);
+      if (commonFace == Face()) {
+        // Create new curve
+        newCurves.emplace_back();
+        newCurves.back().isSigned = curve.isSigned;
+      }
+    }
+    // Don't forget the last point
+    newCurves.back().nodes.push_back(curve.nodes[nNodes - 1]);
+  }
+  return newCurves;
+}
+
 Vector<std::complex<double>> SignedHeatSolver::integrateVectorHeatFlow(const std::vector<Curve>& curves,
                                                                        const std::vector<SurfacePoint>& points,
                                                                        const SignedHeatOptions& options) {