Adding Signed Heat Method (& a bit more!)

CMakeLists.txt

@@ -15,6 +15,7 @@ add_subdirectory(deps/geometry-central)
 pybind11_add_module(potpourri3d_bindings
   src/cpp/core.cpp
   src/cpp/io.cpp
+  src/cpp/heat_helpers.cpp
   src/cpp/mesh.cpp
   src/cpp/point_cloud.cpp
 )

README.md

@@ -1,12 +1,12 @@
 # potpourri3d
 
-A Python library of various algorithms and utilities for 3D triangle meshes and point clouds. Managed by [Nicholas Sharp](https://nmwsharp.com), with new tools added lazily as needed. Currently, mainly bindings to C++ tools from [geometry-central](http://geometry-central.net/).
+A Python library of various algorithms and utilities for 3D triangle meshes, polygon meshes, and point clouds. Managed by [Nicholas Sharp](https://nmwsharp.com), with new tools added lazily as needed. Currently, mainly bindings to C++ tools from [geometry-central](http://geometry-central.net/).
 
 `pip install potpourri3d`
 
 The blend includes:
 - Mesh and point cloud reading/writing to a few file formats
-- Use **heat methods** to compute distance, parallel transport, logarithmic maps, and more
+- Use **heat methods** to compute unsigned and signed distances, parallel transport, logarithmic maps, and more
 - Computing geodesic polylines along surface via edge flips
 - More!
 
@@ -29,9 +29,12 @@ python -m pip install potpourri3d --no-binary potpourri3d
 - [Input / Output](#input--output)
 - [Mesh basic utilities](#mesh-basic-utilities)
 - [Mesh Distance](#mesh-distance)
+- [Mesh Signed Distance](#mesh-signed-distance)
+- [Mesh Fast Marching Distance](#mesh-fast-marching-distance)
 - [Mesh Vector Heat](#mesh-vector-heat)
 - [Mesh Geodesic Paths](#mesh-geodesic-paths)
 - [Mesh Geodesic Tracing](#mesh-geodesic-tracing)
+- [Polygon Mesh Distance & Transport](#polygon-mesh-distance--transport)
 - [Point Cloud Distance & Vector Heat](#point-cloud-distance--vector-heat)
 - [Other Point Cloud Routines](#other-point-cloud-routines)
 
@@ -40,6 +43,7 @@ python -m pip install potpourri3d --no-binary potpourri3d
 Read/write meshes and point clouds from some common formats.
 
 - `read_mesh(filename)` Reads a mesh from file. Returns numpy matrices `V, F`, a Nx3 real numpy array of vertices and a Mx3 integer numpy array of 0-based face indices (or Mx4 for a quad mesh, etc).
+- `read_polygon_mesh(filename)` Reads a mesh from file. Returns numpy matrices `V, F`, where `V` is a Nx3 real numpy array of vertices, and a `polygons` is a nested list of integers; each sub-list represents a polygon face with 0-based face indices.
   - `filename` the path to read the file from. Currently supports the same file types as [geometry-central](http://geometry-central.net/surface/utilities/io/#supported-file-types). The file type is inferred automatically from the path extension.
 
 - `write_mesh(V, F, filename, UV_coords=None, UV_type=None)` Write a mesh to file, optionally with UV coords.
@@ -61,6 +65,14 @@ Read/write meshes and point clouds from some common formats.
 - `face_areas(V, F)` computes a length-F real numpy array of face areas for a triangular mesh
 - `vertex_areas(V, F)` computes a length-V real numpy array of vertex areas for a triangular mesh (equal to 1/3 the sum of the incident face areas)
 - `cotan_laplacian(V, F, denom_eps=0.)` computes the cotan-Laplace matrix as a VxV real sparse csr scipy matrix. Optionally, set `denom_eps` to a small value like `1e-6` to get some additional stability in the presence of degenerate faces.
+- `edges(V, F)` returns the Ex2 integer-valued matrix representing the edges of the given surface mesh, as constructed internally. The i-th row gives the indices of the i-th edge's two endpoint vertices.
+- _Barycentric points_ are used as input and output to some algorithms below, specified as 2-tuples of the form `(element_index, barycentric_coordinates)`:
+  - Vertices are specified as `(vertex_index, )`
+  - Edges are specified as `(edge_index, [t])` where t ∈ [0,1] is the parameter along the edge
+  - Faces are specified as `(face_index, [tA, tB])` where `tA`, `tB` (and optionally, `tC`) are barycentric coordinates in the face. If `tC` is not specified, then `tC` is inferred to be `1 - tA - tB`.
+- `MarchingTrianglesSolver(V, F)` construct an instance of a solver class for contouring scalar functions on triangle meshes using the marching triangles algorithm. 
+  - `MarchingTrianglesSolver.marching_triangles(u, isoval=0.)` takes as input a vector `u` representing a scalar function defined on mesh vertices, and an isovalue; returns a list of lists of barycentric points, where each sublist represents a single connected curve component.
+  - `marching_triangles(V, F, u, isoval=0.)` is similar to the above, but one-off instead of stateful. Returns a list of lists of barycentric points.
 
 ### Mesh Distance
 
@@ -91,6 +103,82 @@ The heat method works by solving a sequence of linear PDEs on the surface of you
 - `compute_distance(V, F, v_ind)` Similar to above, but one-off instead of stateful. Returns an array of distances.
 - `compute_distance_multisource(V, F, v_ind_list)` Similar to above, but one-off instead of stateful. Returns an array of distances.
 
+### Mesh Signed Distance
+
+Use the [signed heat method](https://nzfeng.github.io/research/SignedHeatMethod/index.html) to compute signed distance on meshes, robust to holes and noise. Repeated solves are fast after initial setup.
+
+```python
+import potpourri3d as pp3d
+
+V, F = # your mesh
+solver = pp3d.MeshSignedHeatSolver(V, F)
+
+# Specify a curve as a sequence of barycentric points
+curves = [
+           [
+             (61, [0.3, 0.3]), # face
+             (7, []), # vertex
+             (16, [0.3, 0.3, 0.4]), # face
+             (11, [0.4]), # edge
+             (71, []), # vertex
+             (20, [0.3, 0.3, 0.4]), # face
+             (13, []), # vertex
+             (58, []) # vertex
+             ]
+         ]
+
+# Compute a distance field combining signed distance to curve sources, and unsigned distance to point sources.
+dist = solver.compute_distance(curves, [], points) 
+```
+
+- `MeshSignedHeatSolver.compute_distance(curves, is_signed, points, preserve_source_normals=False, level_set_constraint="ZeroSet", soft_level_set_weight=-1)`
+  - `curves` a list of lists of source points; each point is specified via barycentric coordinates.
+  - `is_signed` a list of bools, one for each curve in `curves`, indicating whether one should compute signed distance (`True`) or unsigned distance (`False`) to a curve. All `True` by default.
+  - `points` a list of source vertex indices
+  - `preserve_source_normals` whether to additionally constrain the normals of the curve. Generally not necessary.
+  - `level_set_constraint` whether to apply level set constraints, with options "ZeroSet", "None", "Multiple". Generally set to "ZeroSet" (set by default).
+  - `soft_level_set_weight` float; if positive, gives the weight with which the given level set constraint is "softly" enforced (negative by default). Generally not necessary.
+
+### Mesh Fast Marching Distance
+
+```python
+import potpourri3d as pp3d
+
+V, F = # your mesh
+solver = pp3d.MeshFastMarchingDistanceSolver(V, F)
+
+# Specify each curve as a sequence of barycentric points
+curves = [
+           [
+             (61, [0.3, 0.3]), # face
+             (7, []), # vertex
+             (16, [0.3, 0.3, 0.4]), # face
+             (11, [0.4]), # edge
+             (71, []), # vertex
+             (20, [0.3, 0.3, 0.4]), # face
+             (13, []), # vertex
+             (58, []) # vertex
+             ]
+         ]
+
+# Compute a signed distance field to a set of closed curves.
+signed_dist = solver.compute_distance(curves, sign=True) 
+
+# Compute unsigned to a set of points.
+points = [
+          [
+            (71, []), # vertex
+            (18, [0.5]) # edge
+          ]
+         ]
+unsigned_dist = solver.compute_distance(points, sign=False) 
+```
+
+- `MeshFastMarchingDistanceSolver.compute_distance(curves, distances=[], sign=False)`
+  - `curves` a list of lists of source points; each point is specified via barycentric coordinates.
+  - `distances` a list of lists of initial distances. Default initial distances are 0.
+  - `sign` if False, compute unsigned distance; if True, compute signed distance. (When initial distances are not 0, "signed" means that the gradient of distance is continuous across the source curves.)
+
 ### Mesh Vector Heat
 
 Use the [vector heat method](https://nmwsharp.com/research/vector-heat-method/) to compute various interpolation & vector-based quantities on meshes. Repeated solves are fast after initial setup.
@@ -187,10 +275,65 @@ trace_pts = tracer.trace_geodesic_from_vertex(22, np.array((0.3, 0.5, 0.4)))
 - `GeodesicTracer.trace_geodesic_from_face(start_face, bary_coords, direction_xyz, max_iterations=None)` similar to above, but from a point in a face. `bary_coords` is a length-3 vector of barycentric coordinates giving the location within the face to start from.
 
 Set `max_iterations` to terminate early after tracing the path through some number of faces/edges (default: no limit).
+
+### Polygon Mesh Distance & Transport
+
+Use the [heat method for unsigned geodesic distance](https://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/), the [signed heat method](https://nzfeng.github.io/research/SignedHeatMethod/index.html) to compute signed distance, and the [vector heat method](https://nmwsharp.com/research/vector-heat-method/) to compute various interpolation & vector-based quantities on general polygon meshes (including mixed-degree meshes, such as tri-quad meshes). Repeated solves are fast after initial setup.
+
+```python
+import potpourri3d as pp3d
+
+V, polygons = # your polygon mesh
+solver = pp3d.PolygonMeshHeatSolver(V, F)
+
+# Compute unsigned geodesic distance to vertices 12 and 17
+dist = solver.compute_distance([12, 17])
+
+# Extend the value `0.` from vertex 12 and `1.` from vertex 17. Any vertex 
+# geodesically closer to 12. will take the value 0., and vice versa 
+# (plus some slight smoothing)
+ext = solver.extend_scalar([12, 17], [0.,1.])
+
+# Get the tangent frames which are used by the solver to define tangent data
+# at each vertex
+basisX, basisY, basisN = solver.get_tangent_frames()
+
+# Parallel transport a vector along the surface
+# (and map it to a vector in 3D)
+sourceV = 22
+ext = solver.transport_tangent_vector(sourceV, [6., 6.])
+ext3D = ext[:,0,np.newaxis] * basisX +  ext[:,1,np.newaxis] * basisY
+
+# Compute signed distance to the oriented curve(s) denoted by a vertex sequence.
+curves = [
+           [9, 10, 12, 13, 51, 48], 
+           [79, 93, 12, 30, 78, 18, 92], 
+           [90, 84, 19, 91, 82, 81, 83]
+         ]
+signed_dist = solver.compute_signed_distance(curves)
+```
+
+- `PolygonMeshHeatSolver(V, polygons, t_coef=1.)` construct an instance of the solver class.
+  - `V` a Nx3 real numpy array of vertices 
+  - `polygons` a list of lists; each sub-list represents a polygon face with 0-based face indices (integers).
+  - `t_coef` set the time used for short-time heat flow. Generally don't change this. If necessary, larger values may make the solution more stable at the cost of smoothing it out.
+- `PolygonMeshHeatSolver.extend_scalar(v_inds, values)` nearest-geodesic-neighbor interpolate values defined at vertices. Vertices will take the value from the closest source vertex (plus some slight smoothing)
+  - `v_inds` a list of source vertices
+  - `values` a list of scalar values, one for each source vertex
+- `PolygonMeshHeatSolver.get_tangent_frames()` get the coordinate frames used to define tangent data at each vertex. Returned as a tuple of basis-X, basis-Y, and normal axes, each as an Nx3 array. May be necessary for change-of-basis into or out of tangent vector convention.
+- `PolygonMeshHeatSolver.transport_tangent_vectors(v_inds, vectors)` parallel transports a collection of vectors across a surface, such that each vertex takes the vector from its nearest-geodesic-neighbor.
+  - `v_inds` a list of source vertices
+  - `vectors` a list of 2D tangent vectors, one for each source vertex
+- `PolygonMeshHeatSolver.compute_distance(v_inds)` 
+  - `v_inds` a list of source vertices
+- `PolygonMeshHeatSolver.compute_signed_distance(curves, level_set_constraint="ZeroSet")`
+  - `curves` a list of lists of source vertices
+  - `level_set_constraint` whether to apply level set constraints, with options "ZeroSet", "None", "Multiple". Generally set to "ZeroSet" (set by default).
+
 
 ### Point Cloud Distance & Vector Heat
 
-Use the [heat method for geodesic distance](https://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/) and [vector heat method](https://nmwsharp.com/research/vector-heat-method/) to compute various interpolation & vector-based quantities on point clouds. Repeated solves are fast after initial setup.
+Use the [heat method for unsigned geodesic distance](https://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/), the [signed heat method](https://nzfeng.github.io/research/SignedHeatMethod/index.html) to compute signed distance, and the [vector heat method](https://nmwsharp.com/research/vector-heat-method/) to compute various interpolation & vector-based quantities on point clouds. Repeated solves are fast after initial setup.
 
 ![point cloud vector heat examples](https://github.com/nmwsharp/potpourri3d/blob/master/media/point_heat_solvers.jpg)
 
@@ -221,6 +364,14 @@ ext3D = ext[:,0,np.newaxis] * basisX +  ext[:,1,np.newaxis] * basisY
 
 # Compute the logarithmic map
 logmap = solver.compute_log_map(sourceP)
+
+# Signed distance to the oriented curve(s) denoted by a point sequence.
+curves = [
+           [9, 10, 12, 13, 51, 48], 
+           [79, 93, 12, 30, 78, 18, 92], 
+           [90, 84, 19, 91, 82, 81, 83]
+         ]
+signed_dist = solver.compute_signed_distance(curves, basisN)
 ```
 
 - `PointCloudHeatSolver(P, t_coef=1.)` construct an instance of the solver class.
@@ -238,7 +389,12 @@ logmap = solver.compute_log_map(sourceP)
   - `vectors` a list of 2D tangent vectors, one for each source point
 - `PointCloudHeatSolver.compute_log_map(p_ind)` compute the logarithmic map centered at the given source point
   - `p_ind` index of the source point
-
+- `PointCloudHeatSolver.compute_signed_distance(curves, cloud_normals, preserve_source_normals=False, level_set_constraint="ZeroSet", soft_level_set_weight=-1)`
+  - `curves` a list of lists of source point indices
+  - `cloud_normals` a list of 3D normal vectors, one for each point in the point cloud
+  - `preserve_source_normals` whether to additionally constrain the normals of the curve. Generally not necessary.
+  - `level_set_constraint` whether to apply level set constraints, with options "ZeroSet", "None", "Multiple". Generally set to "ZeroSet" (set by default).
+  - `soft_level_set_weight` float; if positive, gives the weight with which the given level set constraint is "softly" enforced (negative by default). Generally not necessary.
 
 
 ### Other Point Cloud Routines

src/cpp/heat_helpers.cpp

@@ -0,0 +1,112 @@
+#include "heat_helpers.h"
+
+SurfacePoint toSurfacePoint(SurfaceMesh& mesh, const std::pair<int64_t, std::vector<double>>& p) {
+
+  size_t elementIndex = p.first;
+  const std::vector<double>& baryCoords = p.second;
+
+  if (baryCoords.size() < 1) {
+    return SurfacePoint(mesh.vertex(elementIndex));
+  } else if (baryCoords.size() == 1) {
+    return SurfacePoint(mesh.edge(elementIndex), baryCoords[0]);
+  } else if (baryCoords.size() >= 2) {
+    double tC = (baryCoords.size() == 3) ? baryCoords[2] : 1.0 - baryCoords[0] - baryCoords[1];
+    Vector3 faceCoords = {baryCoords[0], baryCoords[1], tC};
+    return SurfacePoint(mesh.face(elementIndex), faceCoords);
+  } else {
+    throw std::runtime_error("Invalid barycentric coordinates for surface point.");
+  }
+}
+
+std::vector<Curve> toSignedCurves(SurfaceMesh& mesh,
+                                  const std::vector<std::vector<std::pair<int64_t, std::vector<double>>>>& pythonCurves,
+                                  const std::vector<bool>& isSigned) {
+
+  std::vector<Curve> curves;
+
+  for (size_t i = 0; i < pythonCurves.size(); i++) {
+    Curve curve;
+    curve.isSigned = (i < isSigned.size()) ? isSigned[i] : true; // default to signed if not specified
+
+    for (const auto& node : pythonCurves[i]) {
+      curve.nodes.push_back(toSurfacePoint(mesh, node));
+    }
+    curves.push_back(curve);
+  }
+
+  return curves;
+}
+
+std::vector<Curve> toSignedCurves(SurfaceMesh& mesh, const std::vector<std::vector<int64_t>>& pythonCurves,
+                                  const std::vector<bool>& isSigned) {
+
+  std::vector<Curve> curves;
+
+  for (size_t i = 0; i < pythonCurves.size(); i++) {
+    Curve curve;
+    curve.isSigned = (i < isSigned.size()) ? isSigned[i] : true; // default to signed if not specified
+    curve.nodes = toSurfacePoints(mesh, pythonCurves[i]);
+    curves.push_back(curve);
+  }
+
+  return curves;
+}
+
+std::vector<std::vector<Vertex>> toSignedVertices(SurfaceMesh& mesh,
+                                                  const std::vector<std::vector<int64_t>>& pythonCurves) {
+  std::vector<std::vector<Vertex>> curves;
+  for (size_t i = 0; i < pythonCurves.size(); i++) {
+    curves.emplace_back();
+    for (const auto& vIdx : pythonCurves[i]) {
+      curves.back().push_back(mesh.vertex(vIdx));
+    }
+  }
+  return curves;
+}
+
+std::vector<SurfacePoint> toSurfacePoints(SurfaceMesh& mesh,
+                                          const std::vector<std::pair<int64_t, std::vector<double>>>& pythonPoints) {
+
+  std::vector<SurfacePoint> points;
+
+  for (const auto& point : pythonPoints) {
+    points.push_back(toSurfacePoint(mesh, point));
+  }
+
+  return points;
+}
+
+std::vector<SurfacePoint> toSurfacePoints(SurfaceMesh& mesh, const std::vector<int64_t>& pythonPoints) {
+
+  std::vector<SurfacePoint> points;
+
+  for (const auto& vIdx : pythonPoints) {
+    points.emplace_back(mesh.vertex(vIdx));
+  }
+
+  return points;
+}
+
+SignedHeatOptions toSignedHeatOptions(bool preserveSourceNormals, const std::string& levelSetConstraint,
+                                      double softLevelSetWeight) {
+
+  auto toLower = [&](const std::string& s) -> std::string {
+    std::string t = s;
+    std::transform(t.begin(), t.end(), t.begin(), [](unsigned char c) { return std::tolower(c); });
+    return t;
+  };
+
+  SignedHeatOptions options;
+  options.preserveSourceNormals = preserveSourceNormals;
+  if (toLower(levelSetConstraint) == "none") {
+    options.levelSetConstraint = LevelSetConstraint::None;
+  }
+  if (toLower(levelSetConstraint) == "zeroset") {
+    options.levelSetConstraint = LevelSetConstraint::ZeroSet;
+  }
+  if (toLower(levelSetConstraint) == "multiple") {
+    options.levelSetConstraint = LevelSetConstraint::Multiple;
+  }
+  options.softLevelSetWeight = softLevelSetWeight;
+  return options;
+}