diff --git a/CMakeLists.txt b/CMakeLists.txt index 1ebc7b0..b341616 100644 --- a/CMakeLists.txt +++ b/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 ) diff --git a/README.md b/README.md index c331f68..d38f0d3 100644 --- a/README.md +++ b/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 diff --git a/deps/geometry-central b/deps/geometry-central index c37a5b5..5adae5b 160000 --- a/deps/geometry-central +++ b/deps/geometry-central @@ -1 +1 @@ -Subproject commit c37a5b54e22ba29298217209da04b178612ee239 +Subproject commit 5adae5b40079897867727461f74b5c498201f969 diff --git a/src/cpp/heat_helpers.cpp b/src/cpp/heat_helpers.cpp new file mode 100644 index 0000000..5a020a6 --- /dev/null +++ b/src/cpp/heat_helpers.cpp @@ -0,0 +1,112 @@ +#include "heat_helpers.h" + +SurfacePoint toSurfacePoint(SurfaceMesh& mesh, const std::pair>& p) { + + size_t elementIndex = p.first; + const std::vector& 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 toSignedCurves(SurfaceMesh& mesh, + const std::vector>>>& pythonCurves, + const std::vector& isSigned) { + + std::vector 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 toSignedCurves(SurfaceMesh& mesh, const std::vector>& pythonCurves, + const std::vector& isSigned) { + + std::vector 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> toSignedVertices(SurfaceMesh& mesh, + const std::vector>& pythonCurves) { + std::vector> 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 toSurfacePoints(SurfaceMesh& mesh, + const std::vector>>& pythonPoints) { + + std::vector points; + + for (const auto& point : pythonPoints) { + points.push_back(toSurfacePoint(mesh, point)); + } + + return points; +} + +std::vector toSurfacePoints(SurfaceMesh& mesh, const std::vector& pythonPoints) { + + std::vector 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; +} \ No newline at end of file diff --git a/src/cpp/heat_helpers.h b/src/cpp/heat_helpers.h new file mode 100644 index 0000000..7f4fc81 --- /dev/null +++ b/src/cpp/heat_helpers.h @@ -0,0 +1,28 @@ +#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. + +SurfacePoint toSurfacePoint(SurfaceMesh& mesh, const std::pair>& p); + +std::vector toSignedCurves(SurfaceMesh& mesh, + const std::vector>>>& pythonCurves, + const std::vector& isSigned); + +std::vector toSignedCurves(SurfaceMesh& mesh, const std::vector>& pythonCurves, + const std::vector& isSigned); + +std::vector> toSignedVertices(SurfaceMesh& mesh, + const std::vector>& pythonCurves); + +std::vector toSurfacePoints(SurfaceMesh& mesh, + const std::vector>>& pythonPoints); + +std::vector toSurfacePoints(SurfaceMesh& mesh, const std::vector& pythonPoints); + +SignedHeatOptions toSignedHeatOptions(bool preserveSourceNormals = false, + const std::string& levelSetConstraint = "ZeroSet", + double softLevelSetWeight = -1.); \ No newline at end of file diff --git a/src/cpp/io.cpp b/src/cpp/io.cpp index 61ea105..610f03b 100644 --- a/src/cpp/io.cpp +++ b/src/cpp/io.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include "Eigen/Dense" @@ -57,6 +58,33 @@ std::tuple, DenseMatrix> read_mesh(std::string file return std::make_tuple(V, F); } +std::tuple, std::vector>> 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 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> 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& verts, const DenseMatrix& faces, const DenseMatrix& corner_UVs) { @@ -178,6 +206,7 @@ void write_point_cloud(DenseMatrix 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")); diff --git a/src/cpp/mesh.cpp b/src/cpp/mesh.cpp index 1a98eb6..8c9f4d1 100644 --- a/src/cpp/mesh.cpp +++ b/src/cpp/mesh.cpp @@ -1,9 +1,12 @@ #include "geometrycentral/numerical/linear_algebra_utilities.h" #include "geometrycentral/surface/edge_length_geometry.h" +#include "geometrycentral/surface/fast_marching_method.h" #include "geometrycentral/surface/flip_geodesics.h" #include "geometrycentral/surface/heat_method_distance.h" #include "geometrycentral/surface/manifold_surface_mesh.h" +#include "geometrycentral/surface/marching_triangles.h" #include "geometrycentral/surface/mesh_graph_algorithms.h" +#include "geometrycentral/surface/polygon_mesh_heat_solver.h" #include "geometrycentral/surface/simple_polygon_mesh.h" #include "geometrycentral/surface/surface_mesh.h" #include "geometrycentral/surface/surface_mesh_factories.h" @@ -19,6 +22,8 @@ #include "Eigen/Dense" +#include "heat_helpers.h" + namespace py = pybind11; using namespace geometrycentral; @@ -29,6 +34,116 @@ using namespace geometrycentral::surface; template using overload_cast_ = pybind11::detail::overload_cast_impl; +// Return the E x 2 integer-valued matrix representing the edges of a mesh as internal to geometry-central. +DenseMatrix edges(DenseMatrix verts, DenseMatrix faces) { + + std::unique_ptr mesh(new SurfaceMesh(faces)); + DenseMatrix E(mesh->nEdges(), 2); + for (size_t i = 0; i < mesh->nEdges(); i++) { + E(i, 0) = mesh->edge(i).firstVertex().getIndex(); + E(i, 1) = mesh->edge(i).secondVertex().getIndex(); + } + return E; +} + +// Wrapper for FMMDistance that constructs an internal mesh and geometry +class FastMarchingDistanceEigen { + +public: + FastMarchingDistanceEigen(DenseMatrix verts, DenseMatrix faces) { + // Construct the internal mesh and geometry + mesh.reset(new SurfaceMesh(faces)); + geom.reset(new VertexPositionGeometry(*mesh)); + for (size_t i = 0; i < mesh->nVertices(); i++) { + for (size_t j = 0; j < 3; j++) { + geom->inputVertexPositions[i][j] = verts(i, j); + } + } + geom->requireVertexIndices(); + geom->requireEdgeIndices(); + } + + Vector compute_distance(std::vector>>> curves, + std::vector> distances, bool sign = false) { + + size_t nCurves = curves.size(); + std::vector>> initialDistances(nCurves); + std::vector> initDists(nCurves); + for (size_t i = 0; i < nCurves; i++) initDists[i] = std::vector(curves[i].size(), 0.); + + for (size_t i = 0; i < std::min(distances.size(), nCurves); i++) { + for (size_t j = 0; j < std::min(distances[i].size(), curves[i].size()); j++) { + initDists[i][j] = distances[i][j]; + } + } + for (size_t i = 0; i < nCurves; i++) { + for (size_t j = 0; j < curves[i].size(); j++) { + initialDistances[i].emplace_back(toSurfacePoint(*mesh, curves[i][j]), initDists[i][j]); + } + } + VertexData phi = FMMDistance(*geom, initialDistances, sign); + return phi.toVector(); + } + +private: + std::unique_ptr mesh; + std::unique_ptr geom; +}; + +// Wrapper for marching triangles that constructs an internal mesh and geometry +class MarchingTrianglesEigen { + +public: + MarchingTrianglesEigen(DenseMatrix verts, DenseMatrix faces) { + // Construct the internal mesh and geometry + mesh.reset(new SurfaceMesh(faces)); + geom.reset(new VertexPositionGeometry(*mesh)); + for (size_t i = 0; i < mesh->nVertices(); i++) { + for (size_t j = 0; j < 3; j++) { + geom->inputVertexPositions[i][j] = verts(i, j); + } + } + + geom->requireVertexIndices(); + geom->requireEdgeIndices(); + geom->requireFaceIndices(); + } + + std::vector>>> marching_triangles(DenseMatrix u, + double isoval) { + + VertexData f(*mesh, u); + std::vector> curves = marchingTriangles(*geom, f, isoval); + size_t nCurves = curves.size(); + std::vector>>> contour(nCurves); + for (size_t i = 0; i < nCurves; i++) { + std::pair> pt; + for (size_t j = 0; j < curves[i].size(); j++) { + const SurfacePoint& p = curves[i][j]; + switch (p.type) { + case (SurfacePointType::Vertex): { + pt.first = geom->vertexIndices[p.vertex]; + break; + } + case (SurfacePointType::Edge): { + pt.first = geom->edgeIndices[p.edge]; + pt.second = {p.tEdge}; + break; + } + default: { + break; + } + } + contour[i].push_back(pt); + } + } + return contour; + } + +private: + std::unique_ptr mesh; + std::unique_ptr geom; +}; // A wrapper class for the heat method solver, which exposes Eigen in/out class HeatMethodDistanceEigen { @@ -165,6 +280,131 @@ class VectorHeatMethodEigen { std::unique_ptr solver; }; +// A wrapper class for the signed heat method solver, which exposes Eigen in/out +class SignedHeatMethodEigen { + + // TODO use intrinsic triangulations here + +public: + SignedHeatMethodEigen(DenseMatrix verts, DenseMatrix faces, double tCoef = 1.0) { + + // Construct the internal mesh and geometry + mesh.reset(new SurfaceMesh(faces)); + geom.reset(new VertexPositionGeometry(*mesh)); + for (size_t i = 0; i < mesh->nVertices(); i++) { + for (size_t j = 0; j < 3; j++) { + geom->inputVertexPositions[i][j] = verts(i, j); + } + } + + // Build the solver + solver.reset(new SignedHeatSolver(*geom, tCoef)); + } + + Vector compute_distance(const std::vector>>>& curves, + const std::vector& isSigned, + const std::vector>>& points, + bool preserveSourceNormals = false, const std::string& levelSetConstraint = "ZeroSet", + double softLevelSetWeight = -1.0) { + + std::vector sourceCurves = toSignedCurves(*mesh, curves, isSigned); + std::vector sourcePoints = toSurfacePoints(*mesh, points); + SignedHeatOptions solveOptions = toSignedHeatOptions(preserveSourceNormals, levelSetConstraint, softLevelSetWeight); + VertexData phi = solver->computeDistance(sourceCurves, sourcePoints, solveOptions); + return phi.toVector(); + } + +private: + std::unique_ptr mesh; + std::unique_ptr geom; + std::unique_ptr solver; +}; + +// A wrapper class for the polygon mesh heat method solver, which exposes Eigen in/out +class PolygonMeshHeatMethodEigen { + +public: + PolygonMeshHeatMethodEigen(const DenseMatrix& verts, const std::vector>& faces, + double tCoef = 1.0) { + + // Construct the internal mesh and geometry + mesh.reset(new SurfaceMesh(faces)); + geom.reset(new VertexPositionGeometry(*mesh)); + for (size_t i = 0; i < mesh->nVertices(); i++) { + for (size_t j = 0; j < 3; j++) { + geom->inputVertexPositions[i][j] = verts(i, j); + } + } + + // Build the solver + solver.reset(new PolygonMeshHeatSolver(*geom, tCoef)); + } + + Vector compute_distance(Vector sourceVerts) { + std::vector sources; + for (size_t i = 0; i < sourceVerts.rows(); i++) { + sources.push_back(mesh->vertex(sourceVerts(i))); + } + VertexData dist = solver->computeDistance(sources); + return dist.toVector(); + } + + Vector extend_scalar(Vector sourceVerts, Vector values) { + std::vector> sources; + for (size_t i = 0; i < sourceVerts.rows(); i++) { + sources.emplace_back(mesh->vertex(sourceVerts(i)), values(i)); + } + VertexData ext = solver->extendScalars(sources); + return ext.toVector(); + } + + // Returns an extrinsic representation of the tangent frame being used internally, as X/Y/N vectors. + // This is a direct copy of what's already in the VectorHeatMethodEigen class. It's not ideal to repeat code, but this + // also seems like the simplest solution without making breaking changes. + std::tuple, DenseMatrix, DenseMatrix> get_tangent_frames() { + + // Just in case we don't already have it + geom->requireVertexNormals(); + geom->requireVertexTangentBasis(); + + // unpack + VertexData basisX(*mesh); + VertexData basisY(*mesh); + for (Vertex v : mesh->vertices()) { + basisX[v] = geom->vertexTangentBasis[v][0]; + basisY[v] = geom->vertexTangentBasis[v][1]; + } + + return std::tuple, DenseMatrix, DenseMatrix>( + EigenMap(basisX), EigenMap(basisY), EigenMap(geom->vertexNormals)); + } + + DenseMatrix transport_tangent_vectors(Vector sourceVerts, DenseMatrix values) { + std::vector> sources; + for (size_t i = 0; i < sourceVerts.rows(); i++) { + sources.emplace_back(mesh->vertex(sourceVerts(i)), Vector2{values(i, 0), values(i, 1)}); + } + VertexData ext = solver->transportTangentVectors(sources); + + return EigenMap(ext); + } + + Vector compute_signed_distance(const std::vector>& curves, + const std::string& levelSetConstraint = "ZeroSet") { + + std::vector> sourceCurves = toSignedVertices(*mesh, curves); + SignedHeatOptions solveOptions = toSignedHeatOptions(false, levelSetConstraint, -1.0); + VertexData dist = solver->computeSignedDistance(sourceCurves, solveOptions.levelSetConstraint); + return dist.toVector(); + } + + +private: + std::unique_ptr mesh; + std::unique_ptr geom; + std::unique_ptr solver; +}; + // A wrapper class for flip-based geodesics class EdgeFlipGeodesicsManager { @@ -405,6 +645,16 @@ class GeodesicTracer { // clang-format off void bind_mesh(py::module& m) { + m.def("edges", &edges, py::arg("verts"), py::arg("faces")); + + py::class_(m, "MeshFastMarchingDistance") + .def(py::init, DenseMatrix>()) + .def("compute_distance", &FastMarchingDistanceEigen::compute_distance, py::arg("curves"), py::arg("distances"), py::arg("sign")); + + py::class_(m, "MeshMarchingTriangles") + .def(py::init, DenseMatrix>()) + .def("marching_triangles", &MarchingTrianglesEigen::marching_triangles, py::arg("u"), py::arg("isoval")); + py::class_(m, "MeshHeatMethodDistance") .def(py::init, DenseMatrix, double, bool>()) .def("compute_distance", &HeatMethodDistanceEigen::compute_distance, py::arg("source_vert")) @@ -420,6 +670,25 @@ void bind_mesh(py::module& m) { .def("transport_tangent_vectors", &VectorHeatMethodEigen::transport_tangent_vectors, py::arg("source_verts"), py::arg("vectors")) .def("compute_log_map", &VectorHeatMethodEigen::compute_log_map, py::arg("source_vert")); + py::class_(m, "MeshSignedHeatMethod") + .def(py::init, DenseMatrix, double>()) + .def("compute_distance", &SignedHeatMethodEigen::compute_distance, + py::arg("curves") = std::vector>>>(), + py::arg("is_signed") = std::vector(), + py::arg("points") = std::vector>>(), + py::arg("preserve_source_normals") = false, + py::arg("level_set_constraint") = "ZeroSet", + py::arg("soft_level_set_weight") = -1.0); + + py::class_(m, "PolygonMeshHeatSolver") + .def(py::init, std::vector>, double>()) + .def("compute_distance", &PolygonMeshHeatMethodEigen::compute_distance, py::arg("source_verts")) + .def("extend_scalar", &PolygonMeshHeatMethodEigen::extend_scalar, py::arg("source_verts"), py::arg("values")) + .def("get_tangent_frames", &PolygonMeshHeatMethodEigen::get_tangent_frames) + .def("transport_tangent_vectors", &PolygonMeshHeatMethodEigen::transport_tangent_vectors, py::arg("source_verts"), py::arg("vectors")) + .def("compute_signed_distance", &PolygonMeshHeatMethodEigen::compute_signed_distance, + py::arg("curves") = std::vector>(), + py::arg("level_set_constraint") = "ZeroSet"); py::class_(m, "EdgeFlipGeodesicsManager") .def(py::init, DenseMatrix>()) diff --git a/src/cpp/point_cloud.cpp b/src/cpp/point_cloud.cpp index a0f3389..3c9579b 100644 --- a/src/cpp/point_cloud.cpp +++ b/src/cpp/point_cloud.cpp @@ -1,17 +1,19 @@ #include "geometrycentral/pointcloud/point_cloud.h" #include "geometrycentral/numerical/linear_algebra_utilities.h" +#include "geometrycentral/pointcloud/local_triangulation.h" #include "geometrycentral/pointcloud/point_cloud_heat_solver.h" #include "geometrycentral/pointcloud/point_cloud_io.h" #include "geometrycentral/pointcloud/point_position_geometry.h" -#include "geometrycentral/pointcloud/local_triangulation.h" #include "geometrycentral/utilities/eigen_interop_helpers.h" #include #include #include +#include #include "Eigen/Dense" +#include "heat_helpers.h" namespace py = pybind11; @@ -60,6 +62,31 @@ class PointCloudHeatSolverEigen { return dist.toVector(); } + // Solve for signed distance from a collection of curves + Vector compute_signed_distance(const std::vector>& curves, + const DenseMatrix& cloudNormals, bool preserveSourceNormals = false, + const std::string& levelSetConstraint = "ZeroSet", + double softLevelSetWeight = -1.) { + + if (cloudNormals.rows() != cloud->nPoints()) { + throw std::runtime_error("`cloudNormals` should have as many entries as points in the point cloud."); + } + std::vector> sourceCurves; + for (size_t i = 0; i < curves.size(); i++) { + sourceCurves.emplace_back(); + for (const auto& pIdx : curves[i]) { + sourceCurves.back().push_back(cloud->point(pIdx)); + } + } + + PointData pointNormals(*cloud); + for (size_t i = 0; i < cloudNormals.rows(); i++) { + for (int j = 0; j < 3; j++) pointNormals[i][j] = cloudNormals(i, j); + } + SignedHeatOptions options = toSignedHeatOptions(preserveSourceNormals, levelSetConstraint, softLevelSetWeight); + PointData dist = solver->computeSignedDistance(sourceCurves, pointNormals, options); + return dist.toVector(); + } Vector extend_scalar(Vector sourcePoints, Vector values) { std::vector> sources; @@ -69,7 +96,7 @@ class PointCloudHeatSolverEigen { PointData ext = solver->extendScalars(sources); return ext.toVector(); } - + // Returns an extrinsic representation of the tangent frame being used internally, as X/Y/N vectors. std::tuple, DenseMatrix, DenseMatrix> get_tangent_frames() { @@ -125,7 +152,8 @@ class PointCloudHeatSolverEigen { // A class that exposes the local pointcloud trinagulation to python class PointCloudLocalTriangulation { public: - PointCloudLocalTriangulation(DenseMatrix points, bool withDegeneracyHeuristic) : withDegeneracyHeuristic(withDegeneracyHeuristic) { + PointCloudLocalTriangulation(DenseMatrix points, bool withDegeneracyHeuristic) + : withDegeneracyHeuristic(withDegeneracyHeuristic) { // Construct the internal cloud and geometry cloud.reset(new PointCloud(points.rows())); @@ -139,7 +167,8 @@ class PointCloudLocalTriangulation { Eigen::Matrix get_local_triangulation() { - PointData>> local_triangulation = buildLocalTriangulations(*cloud, *geom, withDegeneracyHeuristic); + PointData>> local_triangulation = + buildLocalTriangulations(*cloud, *geom, withDegeneracyHeuristic); int max_neigh = 0; @@ -157,10 +186,10 @@ class PointCloudLocalTriangulation { for (Point v : cloud->points()) { int i = 0; - for (auto const &neighs : local_triangulation[v]) { - out(v.getIndex(), i + 0) = (int) neighs[0].getIndex(); - out(v.getIndex(), i + 1) = (int) neighs[1].getIndex(); - out(v.getIndex(), i + 2) = (int) neighs[2].getIndex(); + for (auto const& neighs : local_triangulation[v]) { + out(v.getIndex(), i + 0) = (int)neighs[0].getIndex(); + out(v.getIndex(), i + 1) = (int)neighs[1].getIndex(); + out(v.getIndex(), i + 2) = (int)neighs[2].getIndex(); i += 3; } } @@ -188,7 +217,13 @@ void bind_point_cloud(py::module& m) { .def("get_tangent_frames", &PointCloudHeatSolverEigen::get_tangent_frames) .def("transport_tangent_vector", &PointCloudHeatSolverEigen::transport_tangent_vector, py::arg("source_point"), py::arg("vector")) .def("transport_tangent_vectors", &PointCloudHeatSolverEigen::transport_tangent_vectors, py::arg("source_points"), py::arg("vectors")) - .def("compute_log_map", &PointCloudHeatSolverEigen::compute_log_map, py::arg("source_point")); + .def("compute_log_map", &PointCloudHeatSolverEigen::compute_log_map, py::arg("source_point")) + .def("compute_signed_distance", &PointCloudHeatSolverEigen::compute_signed_distance, + py::arg("curves") = std::vector>(), + py::arg("cloud_normals") = Eigen::Matrix(), + py::arg("preserve_source_normals") = false, + py::arg("level_set_constraint") = "ZeroSet", + py::arg("soft_level_set_weight") = -1.0); py::class_(m, "PointCloudLocalTriangulation") .def(py::init, bool>()) diff --git a/src/potpourri3d/core.py b/src/potpourri3d/core.py index 6fca3a7..19f46f6 100644 --- a/src/potpourri3d/core.py +++ b/src/potpourri3d/core.py @@ -9,14 +9,18 @@ def validate_mesh(V, F, force_triangular=False, test_indices=False): raise ValueError("vertices should be a 2d Nx3 numpy array") n_vert = V.shape[0] - if len(F.shape) != 2 or F.shape[1] < 3: - raise ValueError("faces should be a 2d NxD numpy array, where D >= 3") + if isinstance(F, np.ndarray): + if len(F.shape) != 2 or F.shape[1] < 3: + raise ValueError("faces should be a 2d NxD numpy array, where D >= 3") - if force_triangular and F.shape[1] != 3: - raise ValueError("faces must be triangular; dimensions should be Nx3") + if force_triangular: + if isinstance(F, np.ndarray) and F.shape[1] != 3: + raise ValueError("faces must be triangular; dimensions should be Nx3") + elif not isinstance(F, np.ndarray) and all(len(polygon) == 3 for polygon in F): + raise ValueError("faces must be triangular") if test_indices: - max_elem = np.amin(F) + max_elem = np.amin(F) if isinstance(F, np.ndarray) else max(max(F)) if max_elem >= n_vert: raise ValueError("There is an out-of-bounds face index. Faces should be zero-based array of indices to vertices") diff --git a/src/potpourri3d/io.py b/src/potpourri3d/io.py index 4f13780..a93a1a8 100644 --- a/src/potpourri3d/io.py +++ b/src/potpourri3d/io.py @@ -9,6 +9,11 @@ def read_mesh(filename): F = np.ascontiguousarray(F) return V, F +def read_polygon_mesh(filename): + V, F = pp3db.read_polygon_mesh(filename) + V = np.ascontiguousarray(V) + return V, F + def write_mesh(V, F, filename, UV_coords=None, UV_type=None): # TODO generalize this to take indexed UVs # (the underlying geometry-central writer needs to support it first) diff --git a/src/potpourri3d/mesh.py b/src/potpourri3d/mesh.py index 8972386..bf9bb5c 100644 --- a/src/potpourri3d/mesh.py +++ b/src/potpourri3d/mesh.py @@ -6,6 +6,28 @@ from .core import * +class MeshFastMarchingDistanceSolver(): + + def __init__(self, V, F): + validate_mesh(V, F, force_triangular=True, test_indices=True) + self.bound_solver = pp3db.MeshFastMarchingDistance(V, F) + + def compute_distance(self, curves, distances=[], sign=False): + return self.bound_solver.compute_distance(curves, distances, sign) + +class MeshMarchingTrianglesSolver(): + + def __init__(self, V, F): + validate_mesh(V, F, force_triangular=True, test_indices=True) + self.bound_solver = pp3db.MeshMarchingTriangles(V, F) + + def marching_triangles(self, u, isoval=0.): + return self.bound_solver.marching_triangles(u, isoval) + +def marching_triangles(V, F, u, isoval=0.): + solver = MeshMarchingTrianglesSolver(V, F) + return solver.marching_triangles(u, isoval) + class MeshHeatMethodDistanceSolver(): def __init__(self, V, F, t_coef=1., use_robust=True): @@ -57,6 +79,53 @@ def transport_tangent_vectors(self, v_inds, vectors): def compute_log_map(self, v_ind): return self.bound_solver.compute_log_map(v_ind) + +class MeshSignedHeatSolver(): + + def __init__(self, V, F, t_coef=1.): + validate_mesh(V, F, force_triangular=True, test_indices=True) + self.bound_solver = pp3db.MeshSignedHeatMethod(V, F, t_coef) + + def compute_distance(self, curves, curve_signs=[], points=[], + preserve_source_normals=False, level_set_constraint="ZeroSet", soft_level_set_weight=-1): + + ''' + Args: + curves [list]: List of curves. Each curve is a list of tuples (element_indices, barycentric_coords). + curve_signs [list]: list of bools, indicating whether each curve in `curves' + is signed (true) or unsigned (false). + By default, curves are assumed to be signed (oriented). + ''' + return self.bound_solver.compute_distance(curves, curve_signs, points, + preserve_source_normals, level_set_constraint, soft_level_set_weight) + + +class PolygonMeshHeatSolver(): + + def __init__(self, V, F, t_coef=1.): + validate_mesh(V, F, force_triangular=False, test_indices=True) + self.bound_solver = pp3db.PolygonMeshHeatSolver(V, F, t_coef) + + def compute_distance(self, v_inds): + return self.bound_solver.compute_distance(v_inds) + + def extend_scalar(self, v_inds, values): + if len(v_inds) != len(values): + raise ValueError("source vertex indices and values array should be same shape") + return self.bound_solver.extend_scalar(v_inds, values) + + def get_tangent_frames(self): + return self.bound_solver.get_tangent_frames() + + def transport_tangent_vectors(self, v_inds, vectors): + if len(v_inds) != len(vectors): + raise ValueError("source vertex indices and values array should be same length") + return self.bound_solver.transport_tangent_vectors(v_inds, vectors) + + def compute_signed_distance(self, curves, level_set_constraint="ZeroSet"): + return self.bound_solver.compute_signed_distance(curves, level_set_constraint) + + class EdgeFlipGeodesicSolver(): def __init__(self, V, F, t_coef=1.): @@ -179,3 +248,7 @@ def vertex_areas(V, F): vertex_area /= 3. return vertex_area + +def edges(V, F): + validate_mesh(V, F, force_triangular=False) + return pp3db.edges(V, F) diff --git a/src/potpourri3d/point_cloud.py b/src/potpourri3d/point_cloud.py index fb65874..2d1ceb9 100644 --- a/src/potpourri3d/point_cloud.py +++ b/src/potpourri3d/point_cloud.py @@ -36,6 +36,11 @@ def transport_tangent_vectors(self, p_inds, vectors): def compute_log_map(self, p_ind): return self.bound_solver.compute_log_map(p_ind) + def compute_signed_distance(self, curves, cloud_normals, + preserve_source_normals=False, level_set_constraint="ZeroSet", soft_level_set_weight=-1): + return self.bound_solver.compute_signed_distance(curves, cloud_normals, + preserve_source_normals, level_set_constraint, soft_level_set_weight) + class PointCloudLocalTriangulation(): diff --git a/test/sample.py b/test/sample.py index 6077dcb..361156c 100644 --- a/test/sample.py +++ b/test/sample.py @@ -39,24 +39,82 @@ # Vector heat (transport vector) ext = solver.transport_tangent_vector(1, [6., 6.]) -ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY +ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY ps_mesh.add_vector_quantity("transport vec", ext3D) ext = solver.transport_tangent_vectors([1, 22], [[6., 6.], [3., 4.]]) -ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY +ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY ps_mesh.add_vector_quantity("transport vec2", ext3D) # Vector heat (log map) logmap = solver.compute_log_map(1) ps_mesh.add_parameterization_quantity("logmap", logmap) +# Signed heat +curve_idxs = [809, 808, 830, 783, 843, 181, 696, 781, 302, 937, 793, 26, 1025, 702, 1047, 1088, 487, 981, 1239, 1238, 402, 57, 1189, 161, 471, 529, 560, 316, 766, 833] +broken_curves = [ + [809, 808, 830, 783, 843, 181, 696, 781, 302, 937, 793], + [1088, 487, 981, 1239, 1238, 402, 57, 1189], + [529, 560, 316, 766, 833] + ] +point_idxs = [1085, 1000] + +# Convert to SurfacePoint format +points = [(idx, []) for idx in point_idxs] +curves = [[(idx, []) for idx in curve] for curve in broken_curves] +sp_curves = [ + [(603, [0.3, 0.3, 0.4]), (726, []), (1660, [0.3, 0.3, 0.4]), (718, []), (2058, [0.3, 0.3, 0.4]), (162, []), (583, []), (88, []), (26, []), (1025, []), (600, []), (1394, [0.3, 0.3]), (1101, []), (1517, [0.3, 0.3]), (919, []), (62, [])] + ] + +E = pp3d.edges(V, F) +def interpolate_coordinates(curve): + positions = [] + for point in curve: + (idx, coords) = point + if len(coords) < 1: + positions.append(V[idx,:]) + elif len(coords) == 1: + pos = (1.-coords[0]) * V[E[idx,0],:] + coords[0] * V[E[idx,1],:] + positions.append(pos) + elif len(coords) >= 2: + pos = coords[0]*V[F[idx,0],:] + coords[1]*V[F[idx,1],:] + (1.-coords[0]-coords[1])*V[F[idx,2],:] + positions.append(pos) + return positions + +sp_curves_positions = interpolate_coordinates(sp_curves[0]) + +signed_solver = pp3d.MeshSignedHeatSolver(V, F) +unsigned_dist = signed_solver.compute_distance(curves, [False, False, False], points) +signed_dist = signed_solver.compute_distance(curves, [], [], + preserve_source_normals=True, level_set_constraint="none", soft_level_set_weight=-1) +sp_dist = signed_solver.compute_distance(sp_curves, [], [], + preserve_source_normals=False, level_set_constraint="ZeroSet", soft_level_set_weight=-1) +mixed_dist = signed_solver.compute_distance(curves, [True, True, True], points, + preserve_source_normals=False, level_set_constraint="ZeroSet", soft_level_set_weight=-1) +isocontour = pp3d.marching_triangles(V, F, signed_dist, 0.1) +contour = interpolate_coordinates(isocontour[0]) +ps.register_curve_network("broken curve 0", np.array(V[broken_curves[0]]), edges="line") +ps.register_curve_network("broken curve 1", np.array(V[broken_curves[1]]), edges="line") +ps.register_curve_network("broken curve 2", np.array(V[broken_curves[2]]), edges="line") +ps.register_curve_network("curve (barycentric points)", np.array(sp_curves_positions), edges="line") +ps.register_curve_network("isocontour", np.array(contour), edges="line") +ps_mesh.add_distance_quantity("unsigned distance", unsigned_dist) +ps_mesh.add_distance_quantity("signed distance", signed_dist) +ps_mesh.add_distance_quantity("signed distance (barycentric)", sp_dist) +ps_mesh.add_distance_quantity("mixed distance", mixed_dist) + +# Fast marching +fmm_solver = pp3d.MeshFastMarchingDistanceSolver(V, F) +fmm_dist = fmm_solver.compute_distance([[(idx, []) for idx in curve_idxs]], sign=True) +ps_mesh.add_distance_quantity("fast marching distance", fmm_dist) + # Flip geodesics path_solver = pp3d.EdgeFlipGeodesicSolver(V,F) for k in range(50): - for i in range(5): - path_pts = path_solver.find_geodesic_path(v_start=1, v_end=22+i) - ps.register_curve_network("flip path " + str(i), path_pts, edges='line') - + for i in range(5): + path_pts = path_solver.find_geodesic_path(v_start=1, v_end=22+i) + ps.register_curve_network("flip path " + str(i), path_pts, edges='line') + path_pts = path_solver.find_geodesic_path_poly([1173, 148, 870, 898]) ps.register_curve_network("flip path poly", path_pts, edges='line') @@ -76,6 +134,49 @@ trace_pts = tracer.trace_geodesic_from_face(31, np.array((0.1, 0.4, 0.5)), np.array((0.3, 0.5, 0.4))) ps.register_curve_network("trace face geodesic", trace_pts, edges='line') +# Areas +vert_area = pp3d.vertex_areas(V,F) +ps_mesh.add_scalar_quantity("vert area", vert_area) +face_area = pp3d.face_areas(V,F) +ps_mesh.add_scalar_quantity("face area", face_area, defined_on='faces') + + +# Laplacian +L = pp3d.cotan_laplacian(V,F,denom_eps=1e-6) +M = scipy.sparse.diags(vert_area) +k_eig = 6 +evals, evecs = scipy.sparse.linalg.eigsh(L, k_eig, M, sigma=1e-8) +for i in range(k_eig): + ps_mesh.add_scalar_quantity("evec " + str(i), evecs[:,i]) + +## = Polygon mesh test + +V, F = pp3d.read_polygon_mesh("bunny_poly.ply") +ps_polymesh = ps.register_surface_mesh("polygon mesh", V, F) +solver = pp3d.PolygonMeshHeatSolver(V, F) + +# Distance +dist = solver.compute_distance([100, 1000]) +ps_polymesh.add_distance_quantity("unsigned distance", dist) + +# Vector heat (extend scalar) +ext = solver.extend_scalar([100, 1000], [0., 6.]) +ps_polymesh.add_scalar_quantity("ext", ext) + +# Vector heat (transport vector) +basisX, basisY, basisN = solver.get_tangent_frames() +ext = solver.transport_tangent_vectors([100, 1000], [[6., 6.], [3., 4.]]) +ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY +ps_polymesh.add_vector_quantity("transport vec2", ext3D) + +# Signed heat +curves = [ + [96, 1037, 1238, 1239, 851, 487, 1088, 601, 217, 942, 702, 631, 805, 793, 937, 1237, 302, 781, 185, 923, 908, 843, 909, 910, 889, 809, 833, 766, 317, 267, 296, 345, 289, 484, 465, 985, 1168, 1099, 17, 212] + ] +signed_dist = solver.compute_signed_distance(curves) +ps.register_curve_network("polygon mesh curve", np.array(V[curves[0]]), edges="loop") +ps_polymesh.add_distance_quantity("signed distance", signed_dist) + ## = Point cloud test P = V ps_cloud = ps.register_point_cloud("cloud", P) @@ -101,11 +202,11 @@ # Vector heat (transport vector) ext = solver.transport_tangent_vector(1, [6., 6.]) -ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY +ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY ps_cloud.add_vector_quantity("transport vec", ext3D) ext = solver.transport_tangent_vectors([1, 22], [[6., 6.], [3., 4.]]) -ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY +ext3D = ext[:,0,np.newaxis] * basisX + ext[:,1,np.newaxis] * basisY ps_cloud.add_vector_quantity("transport vec2", ext3D) # Vector heat (log map) @@ -113,19 +214,12 @@ ps_cloud.add_scalar_quantity("logmapX", logmap[:,0]) ps_cloud.add_scalar_quantity("logmapY", logmap[:,1]) -# Areas -vert_area = pp3d.vertex_areas(V,F) -ps_mesh.add_scalar_quantity("vert area", vert_area) -face_area = pp3d.face_areas(V,F) -ps_mesh.add_scalar_quantity("face area", face_area, defined_on='faces') - - -# Laplacian -L = pp3d.cotan_laplacian(V,F,denom_eps=1e-6) -M = scipy.sparse.diags(vert_area) -k_eig = 6 -evals, evecs = scipy.sparse.linalg.eigsh(L, k_eig, M, sigma=1e-8) -for i in range(k_eig): - ps_mesh.add_scalar_quantity("evec " + str(i), evecs[:,i]) +# Signed heat +curves = [ + [809, 808, 830, 783, 843, 181, 696, 781, 302, 937, 793, 26, 1025, 704, 1047, 1088, 487, 981, 1239, 1238, 402, 57, 1189, 161, 471, 529, 560, 316, 766, 832] + ] +signed_dist = solver.compute_signed_distance(curves, basisN, + preserve_source_normals=False, level_set_constraint="ZeroSet", soft_level_set_weight=-1) +ps_cloud.add_scalar_quantity("signed distance", signed_dist) ps.show()