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 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,11 @@ 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 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 +42,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.
@@ -91,6 +94,49 @@ 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
+  * Vertices: (vertex_index, )
+  * Edges: (edge_index, [t]) where t ∈ [0,1] is the parameter along the edge
+  * Faces: (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.
+'''
+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 = signed_solver.compute_distance(curves, [], points) 
+
+```
+
+- `MeshSignedHeatMethod.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 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 +233,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 +322,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 +347,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,120 @@
+#include "heat_helpers.h"
+
+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]) {
+      size_t elementIndex = node.first;
+      const std::vector<double>& baryCoords = node.second;
+
+      if (baryCoords.size() < 1) {
+        curve.nodes.emplace_back(mesh.vertex(elementIndex));
+      } else if (baryCoords.size() == 1) {
+        curve.nodes.emplace_back(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};
+        curve.nodes.emplace_back(mesh.face(elementIndex), faceCoords);
+      } else {
+        throw std::runtime_error("Invalid barycentric coordinates for a curve 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) {
+    size_t elementIndex = point.first;
+    const std::vector<double>& baryCoords = point.second;
+
+    if (baryCoords.size() < 1) {
+      points.emplace_back(mesh.vertex(elementIndex));
+    } else if (baryCoords.size() == 1) {
+      points.emplace_back(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};
+      points.emplace_back(mesh.face(elementIndex), faceCoords);
+    } else {
+      throw std::runtime_error("Invalid barycentric coordinates for surface 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;
+}

src/cpp/heat_helpers.h

@@ -0,0 +1,26 @@
+#include "geometrycentral/surface/signed_heat_method.h"
+
+using namespace geometrycentral;
+using namespace geometrycentral::surface;
+
+// Convert from intermediary representation -- used for easier Python bindings -- to the input structs used in
+// geometry-central.
+
+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> toSignedCurves(SurfaceMesh& mesh, const std::vector<std::vector<int64_t>>& pythonCurves,
+                                  const std::vector<bool>& isSigned);
+
+std::vector<std::vector<Vertex>> toSignedVertices(SurfaceMesh& mesh,
+                                                  const std::vector<std::vector<int64_t>>& pythonCurves);
+
+std::vector<SurfacePoint> toSurfacePoints(SurfaceMesh& mesh,
+                                          const std::vector<std::pair<int64_t, std::vector<double>>>& pythonPoints);
+
+std::vector<SurfacePoint> toSurfacePoints(SurfaceMesh& mesh, const std::vector<int64_t>& pythonPoints);
+
+SignedHeatOptions toSignedHeatOptions(bool preserveSourceNormals = false,
+                                      const std::string& levelSetConstraint = "ZeroSet",
+                                      double softLevelSetWeight = -1.);

src/cpp/io.cpp

@@ -14,6 +14,7 @@
 #include <pybind11/eigen.h>
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
 
 #include "Eigen/Dense"
 
@@ -57,6 +58,33 @@ std::tuple<DenseMatrix<double>, DenseMatrix<int64_t>> read_mesh(std::string file
   return std::make_tuple(V, F);
 }
 
+std::tuple<DenseMatrix<double>, std::vector<std::vector<int64_t>>> read_polygon_mesh(std::string filename) {
+
+  // Call the mesh reader
+  SimplePolygonMesh pmesh(filename);
+
+  if (pmesh.nFaces() == 0) throw std::runtime_error("read mesh has no faces");
+
+  // Manually copy the vertex array
+  DenseMatrix<double> V(pmesh.nVertices(), 3);
+  for (size_t i = 0; i < pmesh.nVertices(); i++) {
+    for (size_t j = 0; j < 3; j++) {
+      V(i, j) = pmesh.vertexCoordinates[i][j];
+    }
+  }
+
+  std::vector<std::vector<int64_t>> polygons(pmesh.nFaces());
+  for (size_t i = 0; i < pmesh.nFaces(); i++) {
+    size_t fDegree = pmesh.polygons[i].size();
+    polygons[i].resize(fDegree);
+    for (size_t j = 0; j < fDegree; j++) {
+      polygons[i][j] = pmesh.polygons[i][j];
+    }
+  }
+
+  return std::make_tuple(V, polygons);
+}
+
 namespace { // anonymous helers
 SimplePolygonMesh buildMesh(const DenseMatrix<double>& verts, const DenseMatrix<int64_t>& faces,
                             const DenseMatrix<double>& corner_UVs) {
@@ -178,6 +206,7 @@ void write_point_cloud(DenseMatrix<double> points, std::string filename) {
 void bind_io(py::module& m) {
 
   m.def("read_mesh", &read_mesh, "Read a mesh from file.", py::arg("filename"));
+  m.def("read_polygon_mesh", &read_polygon_mesh, "Read a polygon mesh from file.", py::arg("filename"));
 
   m.def("write_mesh", &write_mesh, "Write a mesh to file.", py::arg("verts"), py::arg("faces"), py::arg("filename"));
   m.def("write_mesh_pervertex_uv", &write_mesh_pervertex_uv, "Write a mesh to file.", py::arg("verts"), py::arg("faces"), py::arg("UVs"), py::arg("filename"));