diff --git a/docs/docs/surface/algorithms/vector_heat_method.md b/docs/docs/surface/algorithms/vector_heat_method.md index aebfca55..aae57e30 100644 --- a/docs/docs/surface/algorithms/vector_heat_method.md +++ b/docs/docs/surface/algorithms/vector_heat_method.md @@ -2,7 +2,7 @@ This section describes the _Vector Heat Method_ in geometry-central, which compu Note that these quantities all depend on the _intrinsic_ geometry of a surface (via the `IntrinsicGeometryInterface`). Therefore, these routines can be run on abstract geometric domains as well as traditional surfaces in 3D. -These algorithms are described in [The Vector Heat Method](http://www.cs.cmu.edu/~kmcrane/Projects/VectorHeatMethod/paper.pdf). +These algorithms are described in [The Vector Heat Method](http://www.cs.cmu.edu/~kmcrane/Projects/VectorHeatMethod/paper.pdf) and [The Affine Heat Method](https://www.yousufsoliman.com/projects/the-affine-heat-method.html). `#include "geometrycentral/surface/vector_heat_method.h"` @@ -50,11 +50,11 @@ for (/* ... some inputs ... */ ) { VertexData scalarExtension = vhmSolver->extendScalar(points); ``` -??? func "`#!cpp VertexData VectorHeatSolver::extendScalar( const std::vector>& sources)`" +??? func "`#!cpp VertexData VectorHeatSolver::extendScalar(const std::vector>& sources)`" Compute the nearest-neighbor extension of scalar data defined at isolated vertices to the entire domain. The input is a list of vertices and their corresponding values. -??? func "`#!cpp VertexData VectorHeatSolver::extendScalar( const std::vector>& sources)`" +??? func "`#!cpp VertexData VectorHeatSolver::extendScalar(const std::vector>& sources)`" Compute the nearest-neighbor extension of scalar data defined at isolated points to the entire domain. The input is a list of [surface points](../../utilities/surface_point/) and their corresponding values. @@ -83,25 +83,56 @@ The _logarithmic map_ is a very special 2D local parameterization of a surface a ![octopus logmap](/media/octopus_logmap.jpg){: style="height:300px; display: block; margin-left: auto; margin-right: auto;"} -These routines compute the logarithmic map using the vector heat method. +These routines compute the logarithmic map using the vector heat method or the affine heat method. Several strategies are available, specified by the `LogMapStrategy` enum: -??? func "`#!cpp VertexData VectorHeatSolver::computeLogMap(const Vertex& sourceVert)`" +- `LogMapStrategy::VectorHeat`: the original algorithm from "The Vector Heat Method". Fast, but may have some distortion and warping. +- `LogMapStrategy::AffineLocal`: the fast local algorithm from "The Affine Heat Method". Fast & highly accurate near source. Generates extremely regular coordinates, but may have artifacts if you do not use an intrinsic Delaunay triangulation (see below). +- `LogMapStrategy::AffineAdaptive`: the global algorithm from "The Affine Heat Method". Highest quality. Requires factoring a new matrix each time, so it is significantly slower for repeated solves, though roughly similar cost if just solving once. + +??? func "`#!cpp VertexData VectorHeatSolver::computeLogMap(const Vertex& sourceVert, LogMapStrategy strategy = LogMapStrategy::VectorHeat)`" Compute the logarithmic map with respect to the given source vertex. The angular coordinate of the log map will be respect to the tangent space of the source vertex. -??? func "`#!cpp VertexData VectorHeatSolver::computeLogMap(const SurfacePoint& sourceP)`" +??? func "`#!cpp VertexData VectorHeatSolver::computeLogMap(const SurfacePoint& sourceP, LogMapStrategy strategy = LogMapStrategy::VectorHeat)`" Compute the logarithmic map with respect to the given source point, which is a general [surface point](../../utilities/surface_point/). The angular coordinate of the log map will be respect to the tangent space of the source vertex, edge, or face. +!!! note "intrinsic triangulations make logmaps much more accurate" + + Using an [intrinsic triangulation](/surface/intrinsic_triangulations/basics) will make your logarithmic maps dramatically more accurate, especially on meshes with irregular triangles, at little cost. The `AffineLocal` strategy in particular benefits greatly from operating on a Delaunay intrinsic triangulation. + + **Example: log maps with intrinsic triangulations** + ```cpp + #include <#include "geometrycentral/surface/signpost_intrinsic_triangulation.h"> + + // Construct an intrinsic triangulation + SignpostIntrinsicTriangulation signpostTri(*mesh, *geometry); + signpostTri.flipToDelaunay(); // this is the important step, makes it numerically nice! + + // Remap the vertex to the intrinsic triangulation, if using + Vertex origSourceVert = /* your vertex */; + Vertex intrinsicSourceVert = signpostTri->equivalentPointOnIntrinsic(SurfacePoint(origSourceVert)).vertex; + + // Run the algorithm + VectorHeatMethodSolver solver(signpostTri); + VertexData logmapIntrinsic = solver->computeLogMap(sourceVert, logMapStrategy); + + // Copy the logmap back to your original mesh + VertexData logmap = signpostTri->restrictToInput(logmapIntrinsic); + ``` + + By default, the resulting logarithmic map has coordinates defined in the tangent space of the source point on the intrinsic mesh. At vertices these spaces are the same, but at a general `SurfacePoint` inside some face you might want to align tangent spaces; see the [intrinsic triangulation](/surface/intrinsic_triangulations/basics) documentation for details. There may or may not be a builtin method to do the alignment, depending on your setting. + + ## Citation -If these algorithms contribute to academic work, please cite the following paper: +If these algorithms contribute to academic work, please cite the following paper(s): ```bib @article{sharp2019vector, @@ -114,4 +145,13 @@ If these algorithms contribute to academic work, please cite the following paper year={2019}, publisher={ACM} } + +@article{soliman2025affine, + title={The Affine Heat Method}, + author={Yousuf Soliman, Nicholas Sharp, + booktitle={Computer Graphics Forum}, + volume={44}, + number={5}, + year={2025} +} ``` diff --git a/include/geometrycentral/surface/surface_centers.h b/include/geometrycentral/surface/surface_centers.h index 70382627..9457364f 100644 --- a/include/geometrycentral/surface/surface_centers.h +++ b/include/geometrycentral/surface/surface_centers.h @@ -1,21 +1,27 @@ #pragma once #include "geometrycentral/surface/intrinsic_geometry_interface.h" -#include "geometrycentral/surface/vector_heat_method.h" #include "geometrycentral/surface/manifold_surface_mesh.h" +#include "geometrycentral/surface/vector_heat_method.h" namespace geometrycentral { namespace surface { // Find a center of a collection of points at vertices -SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const std::vector& vertexPts, int p = 2); +SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, + const std::vector& vertexPts, int p = 2, + LogMapStrategy strategy = LogMapStrategy::VectorHeat); SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, VectorHeatMethodSolver& solver, - const std::vector& vertexPts, int p = 2); + const std::vector& vertexPts, int p = 2, + LogMapStrategy strategy = LogMapStrategy::VectorHeat); // Find a center of distribution -SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const VertexData& distribution, int p = 2); +SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, + const VertexData& distribution, int p = 2, + LogMapStrategy strategy = LogMapStrategy::VectorHeat); SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, VectorHeatMethodSolver& solver, - const VertexData& distribution, int p = 2); + const VertexData& distribution, int p = 2, + LogMapStrategy strategy = LogMapStrategy::VectorHeat); } // namespace surface } // namespace geometrycentral diff --git a/include/geometrycentral/surface/vector_heat_method.h b/include/geometrycentral/surface/vector_heat_method.h index dc97179d..786c5810 100644 --- a/include/geometrycentral/surface/vector_heat_method.h +++ b/include/geometrycentral/surface/vector_heat_method.h @@ -19,6 +19,13 @@ namespace surface { // Stateful class. Allows efficient repeated solves +enum class LogMapStrategy { + VectorHeat, // the logmap proposed in the original Vector Heat Method paper + AffineLocal, // the logmap from the Affine Heat Method, allows fast prefactoring for repeated solves, more accurate + // near source + AffineAdaptive, // the logmap from the Affine Heat Method, no prefactoring for repeated solves, but most accurate +}; + class VectorHeatMethodSolver { public: @@ -36,10 +43,9 @@ class VectorHeatMethodSolver { VertexData transportTangentVectors(const std::vector>& sources); VertexData transportTangentVectors(const std::vector>& sources); - // === The Logarithmic map - VertexData computeLogMap(const Vertex& sourceVert, double vertexDistanceShift = 0.); - VertexData computeLogMap(const SurfacePoint& sourceP); + VertexData computeLogMap(const Vertex& sourceVert, LogMapStrategy strategy = LogMapStrategy::VectorHeat); + VertexData computeLogMap(const SurfacePoint& sourceP, LogMapStrategy strategy = LogMapStrategy::VectorHeat); // === Options and parameters const double tCoef; // the time parameter used for heat flow, measured as time = tCoef * mean_edge_length^2 @@ -65,14 +71,21 @@ class VectorHeatMethodSolver { // Solvers std::unique_ptr> scalarHeatSolver; std::unique_ptr>> vectorHeatSolver; + std::unique_ptr> affineHeatSolver; std::unique_ptr> poissonSolver; SparseMatrix massMat; // Helpers void ensureHaveScalarHeatSolver(); void ensureHaveVectorHeatSolver(); + void ensureHaveAffineHeatSolver(); void ensureHavePoissonSolver(); + // Log map approaches, see documentation on the strategy enum + VertexData computeLogMap_VectorHeat(const Vertex& sourceVert, double vertexDistanceShift = 0.); + VertexData computeLogMap_AffineLocal(const Vertex& sourceVert); + VertexData computeLogMap_AffineAdaptive(const Vertex& sourceVert); + void addVertexOutwardBall(Vertex v, Vector>& distGradRHS); }; diff --git a/src/surface/surface_centers.cpp b/src/surface/surface_centers.cpp index 70069499..0044f117 100644 --- a/src/surface/surface_centers.cpp +++ b/src/surface/surface_centers.cpp @@ -10,37 +10,39 @@ namespace geometrycentral { namespace surface { -SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const std::vector& vertexPts, int p) { +SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, + const std::vector& vertexPts, int p, LogMapStrategy strategy) { VertexData dist(geom.mesh, 0.); for (Vertex v : vertexPts) { dist[v] += 1.; } // Forward to the general version - return findCenter(mesh, geom, dist, p); + return findCenter(mesh, geom, dist, p, strategy); } SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, VectorHeatMethodSolver& solver, - const std::vector& vertexPts, int p) { + const std::vector& vertexPts, int p, LogMapStrategy strategy) { VertexData dist(geom.mesh, 0.); for (Vertex v : vertexPts) { dist[v] += 1.; } // Forward to the general version - return findCenter(mesh, geom, solver, dist, p); + return findCenter(mesh, geom, solver, dist, p, strategy); } -SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const VertexData& distribution, int p) { +SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, + const VertexData& distribution, int p, LogMapStrategy strategy) { VectorHeatMethodSolver solver(geom); // Forward to the general version - return findCenter(mesh, geom, solver, distribution, p); + return findCenter(mesh, geom, solver, distribution, p, strategy); } SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, VectorHeatMethodSolver& solver, - const VertexData& distribution, int p) { + const VertexData& distribution, int p, LogMapStrategy strategy) { if (p != 1 && p != 2) { throw std::logic_error("only p=1 or p=2 is supported"); @@ -87,7 +89,7 @@ SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& g auto evalEnergyAndUpdate = [&](SurfacePoint aboutPoint) -> std::tuple { // Compute the current log map // Solve at the face point - VertexData logmap = solver.computeLogMap(aboutPoint); + VertexData logmap = solver.computeLogMap(aboutPoint, strategy); // Evaluate energy and update step double thisEnergy = 0.; diff --git a/src/surface/vector_heat_method.cpp b/src/surface/vector_heat_method.cpp index 65249202..101cdad2 100644 --- a/src/surface/vector_heat_method.cpp +++ b/src/surface/vector_heat_method.cpp @@ -1,4 +1,5 @@ #include "geometrycentral/surface/vector_heat_method.h" +#include "geometrycentral/utilities/vector2.h" namespace geometrycentral { namespace surface { @@ -71,6 +72,67 @@ void VectorHeatMethodSolver::ensureHaveVectorHeatSolver() { geom.unrequireVertexConnectionLaplacian(); } +void VectorHeatMethodSolver::ensureHaveAffineHeatSolver() { + if (affineHeatSolver != nullptr) return; + + geom.requireVertexIndices(); + geom.requireEdgeCotanWeights(); + geom.requireTransportVectorsAlongHalfedge(); + geom.requireHalfedgeVectorsInVertex(); + geom.requireVertexDualAreas(); + + // Get the ingredients + std::vector> triplets; + for (Halfedge he : mesh.halfedges()) { + size_t iTail = geom.vertexIndices[he.vertex()]; + size_t iTip = geom.vertexIndices[he.twin().vertex()]; + + double weight = geom.edgeCotanWeights[he.edge()]; + + Vector2 eij = geom.halfedgeVectorsInVertex[he]; + Vector2 eji = geom.halfedgeVectorsInVertex[he.twin()]; + + Vector2 rot = -eij / eji; + Vector2 t = -eij; + + // Rotation from TjM -> TiM, followed by a translation TjM -> TjM + // The order matters here! + DenseMatrix transport(3, 3); + transport << rot.x, -rot.y, t.x, rot.y, rot.x, t.y, 0, 0, 1.; + + for (size_t p = 0; p < 3; p++) { + triplets.emplace_back(3 * iTail + p, 3 * iTail + p, weight); + for (size_t q = 0; q < 3; q++) { + triplets.emplace_back(3 * iTail + p, 3 * iTip + q, -weight * transport(p, q)); + } + } + } + + // Build the affine connection Laplacian + SparseMatrix L_affineConnection = Eigen::SparseMatrix(3 * mesh.nVertices(), 3 * mesh.nVertices()); + L_affineConnection.setFromTriplets(triplets.begin(), triplets.end()); + + // Build the lumped mass matrix (with 3x identical entries for each 3x3 block) + Vector massVec(3 * mesh.nVertices()); + for (Vertex v : mesh.vertices()) { + for (size_t p = 0; p < 3; p++) { + massVec(3 * geom.vertexIndices[v] + p) = geom.vertexDualAreas[v]; + } + } + SparseMatrix massMat(massVec.asDiagonal()); + + // Build (and factor) the short-time heat flow operator + SparseMatrix L_affineHeatOp = massMat + shortTime * L_affineConnection; + affineHeatSolver.reset( + new SquareSolver(L_affineHeatOp)); // not necessarily SPD in the scalar sense, even with Delaunay + + // Bookkeeping + geom.unrequireVertexIndices(); + geom.unrequireEdgeCotanWeights(); + geom.unrequireTransportVectorsAlongHalfedge(); + geom.unrequireHalfedgeVectorsInVertex(); + geom.unrequireVertexDualAreas(); +} void VectorHeatMethodSolver::ensureHavePoissonSolver() { if (poissonSolver != nullptr) return; @@ -267,8 +329,136 @@ VectorHeatMethodSolver::transportTangentVectors(const std::vector VectorHeatMethodSolver::computeLogMap(const Vertex& sourceVert, LogMapStrategy strategy) { + switch (strategy) { + case LogMapStrategy::VectorHeat: + return computeLogMap_VectorHeat(sourceVert, 0.0); + case LogMapStrategy::AffineLocal: + return computeLogMap_AffineLocal(sourceVert); + case LogMapStrategy::AffineAdaptive: + return computeLogMap_AffineAdaptive(sourceVert); + } +} + +VertexData VectorHeatMethodSolver::computeLogMap_AffineLocal(const Vertex& sourceVert) { + ensureHaveVectorHeatSolver(); + ensureHaveAffineHeatSolver(); + + geom.requireVertexIndices(); + + // Construct the parallel section + VertexData parallel_section = transportTangentVector(sourceVert, Vector2{1, 0}); + + // Build rhs + Vector horizontalRHS = Vector::Zero(3 * mesh.nVertices()); + horizontalRHS[3 * geom.vertexIndices[sourceVert] + 2] += 1.0; // (0,0), in homogeneous coordinates + + // Short-time heat flow on of hte affine values + Vector horizontalSol = affineHeatSolver->solve(horizontalRHS); + + // Construct the map from the heat-flow'd values + VertexData logmap(mesh); + for (Vertex v : mesh.vertices()) { + size_t vInd = geom.vertexIndices[v]; + + Vector2 logCoord = Vector2{horizontalSol(3 * vInd + 0), horizontalSol(3 * vInd + 1)}; + logCoord /= horizontalSol(3 * vInd + 2); // homogenous division + + // Change of basis to coordinate frame at source + logmap[v] = logCoord / parallel_section[v]; + } + + // Bookkeeping + geom.unrequireVertexIndices(); + + return logmap; +} + +VertexData VectorHeatMethodSolver::computeLogMap_AffineAdaptive(const Vertex& sourceVert) { + + geom.requireVertexIndices(); + geom.requireFaceAreas(); + geom.requireVertexDualAreas(); + geom.requireEdgeLengths(); + geom.requireEdgeCotanWeights(); + geom.requireHalfedgeVectorsInVertex(); + + // Construct a parallel section using the coordinate from the source + VertexData parallel_section = transportTangentVector(sourceVert, Vector2{1., 0.}); + + // Build the adaptive connection Laplacian in the coordinate frame + // (this depends on the source, so it cannot be prefactored) + std::vector> triplets; + for (Halfedge he : mesh.halfedges()) { + size_t iTail = geom.vertexIndices[he.vertex()]; + size_t iTip = geom.vertexIndices[he.twin().vertex()]; + + double weight = geom.edgeCotanWeights[he.edge()]; + + // Translation in the edge direction in the trivialization induced by the geodesic frame X + // We have two options on how to do it. Turns out it doesn't really matter which we pick. + + // Do the translation using the values at i + // (we could equivalently do it at j) + Vector2 eij = geom.halfedgeVectorsInVertex[he]; + Vector2 Xi = parallel_section[iTail]; + Vector2 t = -eij / Xi; + + // The developing connection + DenseMatrix transport(3, 3); + transport << 1, 0, t.x, 0, 1, t.y, 0, 0, 1.; + + for (size_t p = 0; p < 3; p++) { + triplets.emplace_back(3 * iTail + p, 3 * iTail + p, weight); + for (size_t q = 0; q < 3; q++) { + triplets.emplace_back(3 * iTail + p, 3 * iTip + q, -weight * transport(p, q)); + } + } + } + SparseMatrix L_affineAdaptive = Eigen::SparseMatrix(3 * mesh.nVertices(), 3 * mesh.nVertices()); + L_affineAdaptive.setFromTriplets(triplets.begin(), triplets.end()); + + // Build the corresponding (lumped) mass matrix, with 3x3 vertex area diagonal blocks + Vector massVec(3 * mesh.nVertices()); + for (Vertex v : mesh.vertices()) + for (size_t p = 0; p < 3; p++) massVec(3 * geom.vertexIndices[v] + p) = geom.vertexDualAreas[v]; + SparseMatrix massMat(massVec.asDiagonal()); + + // Form the short-time heat flow operator and solver + SparseMatrix affineOp = massMat + shortTime * L_affineAdaptive; + SquareSolver affineAdaptiveSolver(affineOp); + + // Build rhs, with a homogenous 0 at the source + Vector horizontalRHS = Vector::Zero(3 * mesh.nVertices()); + horizontalRHS[3 * geom.vertexIndices[sourceVert] + 2] += 1.0; + + // Solve + Vector horizontalSol = affineAdaptiveSolver.solve(horizontalRHS); + + // Compose the result + VertexData logmap(mesh); + for (Vertex v : mesh.vertices()) { + size_t vInd = geom.vertexIndices[v]; + + Vector2 logCoord = Vector2{horizontalSol(3 * vInd + 0), horizontalSol(3 * vInd + 1)}; + logCoord /= horizontalSol(3 * vInd + 2); + + logmap[v] = logCoord; + } + + // Bookkeeping + geom.unrequireVertexIndices(); + geom.unrequireFaceAreas(); + geom.unrequireVertexDualAreas(); + geom.unrequireEdgeLengths(); + geom.unrequireEdgeCotanWeights(); + geom.unrequireHalfedgeVectorsInVertex(); + + return logmap; +} -VertexData VectorHeatMethodSolver::computeLogMap(const Vertex& sourceVert, double vertexDistanceShift) { +VertexData VectorHeatMethodSolver::computeLogMap_VectorHeat(const Vertex& sourceVert, + double vertexDistanceShift) { geom.requireFaceAreas(); geom.requireEdgeLengths(); geom.requireCornerAngles(); @@ -429,14 +619,14 @@ void VectorHeatMethodSolver::addVertexOutwardBall(Vertex vert, Vector VectorHeatMethodSolver::computeLogMap(const SurfacePoint& sourceP) { +VertexData VectorHeatMethodSolver::computeLogMap(const SurfacePoint& sourceP, LogMapStrategy strategy) { geom.requireHalfedgeVectorsInVertex(); geom.requireHalfedgeVectorsInFace(); switch (sourceP.type) { case SurfacePointType::Vertex: { - return computeLogMap(sourceP.vertex); + return computeLogMap(sourceP.vertex, strategy); break; } case SurfacePointType::Edge: { @@ -444,8 +634,8 @@ VertexData VectorHeatMethodSolver::computeLogMap(const SurfacePoint& so // Compute logmaps at both adjacent vertices Halfedge he = sourceP.edge.halfedge(); - VertexData logmapTail = computeLogMap(he.vertex()); - VertexData logmapTip = computeLogMap(he.twin().vertex()); + VertexData logmapTail = computeLogMap(he.vertex(), strategy); + VertexData logmapTip = computeLogMap(he.twin().vertex(), strategy); // Changes of basis Vector2 tailRot = geom.halfedgeVectorsInVertex[he].inv().normalize(); @@ -472,7 +662,7 @@ VertexData VectorHeatMethodSolver::computeLogMap(const SurfacePoint& so for (Halfedge he : sourceP.face.adjacentHalfedges()) { // Commpute logmap at vertex - VertexData logmapVert = computeLogMap(he.vertex()); + VertexData logmapVert = computeLogMap(he.vertex(), strategy); // Compute change of basis to bring it back to the face Vector2 rot = (geom.halfedgeVectorsInFace[he] / geom.halfedgeVectorsInVertex[he]).normalize();