add the affine heat method (#219)

* initial work on affine

* add center option to logmap strategy

* use source initial guess for surface centers

* Revert "use source initial guess for surface centers"

This reverts commit cb3dea1.

* add adaptive affine heat method

* docs for affine heat method

* add caveat about intrinsic tangent space alignment

Author: nmwsharp

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<double> scalarExtension = vhmSolver->extendScalar(points);
 ```
 
-??? func "`#!cpp VertexData<double> VectorHeatSolver::extendScalar( const std::vector<std::tuple<Vertex, double>>& sources)`"
+??? func "`#!cpp VertexData<double> VectorHeatSolver::extendScalar(const std::vector<std::tuple<Vertex, double>>& 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<double> VectorHeatSolver::extendScalar( const std::vector<std::tuple<SurfacePoint, double>>& sources)`"
+??? func "`#!cpp VertexData<double> VectorHeatSolver::extendScalar(const std::vector<std::tuple<SurfacePoint, double>>& 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<Vector2> 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<Vector2> 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<Vector2> VectorHeatSolver::computeLogMap(const SurfacePoint& sourceP)`"
+??? func "`#!cpp VertexData<Vector2> 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<Vector2> logmapIntrinsic = solver->computeLogMap(sourceVert, logMapStrategy);
+
+    // Copy the logmap back to your original mesh
+    VertexData<Vector2> 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}
+}
 ```

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<Vertex>& vertexPts, int p = 2);
+SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom,
+                        const std::vector<Vertex>& vertexPts, int p = 2,
+                        LogMapStrategy strategy = LogMapStrategy::VectorHeat);
 SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, VectorHeatMethodSolver& solver,
-                        const std::vector<Vertex>& vertexPts, int p = 2);
+                        const std::vector<Vertex>& vertexPts, int p = 2,
+                        LogMapStrategy strategy = LogMapStrategy::VectorHeat);
 
 // Find a center of distribution
-SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const VertexData<double>& distribution, int p = 2);
+SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom,
+                        const VertexData<double>& distribution, int p = 2,
+                        LogMapStrategy strategy = LogMapStrategy::VectorHeat);
 SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, VectorHeatMethodSolver& solver,
-                        const VertexData<double>& distribution, int p = 2);
+                        const VertexData<double>& distribution, int p = 2,
+                        LogMapStrategy strategy = LogMapStrategy::VectorHeat);
 
 } // namespace surface
 } // namespace geometrycentral

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<Vector2> transportTangentVectors(const std::vector<std::tuple<Vertex, Vector2>>& sources);
   VertexData<Vector2> transportTangentVectors(const std::vector<std::tuple<SurfacePoint, Vector2>>& sources);
 
-
   // === The Logarithmic map
-  VertexData<Vector2> computeLogMap(const Vertex& sourceVert, double vertexDistanceShift = 0.);
-  VertexData<Vector2> computeLogMap(const SurfacePoint& sourceP);
+  VertexData<Vector2> computeLogMap(const Vertex& sourceVert, LogMapStrategy strategy = LogMapStrategy::VectorHeat);
+  VertexData<Vector2> 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<PositiveDefiniteSolver<double>> scalarHeatSolver;
   std::unique_ptr<LinearSolver<std::complex<double>>> vectorHeatSolver;
+  std::unique_ptr<LinearSolver<double>> affineHeatSolver;
   std::unique_ptr<PositiveDefiniteSolver<double>> poissonSolver;
   SparseMatrix<double> massMat;
 
   // Helpers
   void ensureHaveScalarHeatSolver();
   void ensureHaveVectorHeatSolver();
+  void ensureHaveAffineHeatSolver();
   void ensureHavePoissonSolver();
 
+  // Log map approaches, see documentation on the strategy enum
+  VertexData<Vector2> computeLogMap_VectorHeat(const Vertex& sourceVert, double vertexDistanceShift = 0.);
+  VertexData<Vector2> computeLogMap_AffineLocal(const Vertex& sourceVert);
+  VertexData<Vector2> computeLogMap_AffineAdaptive(const Vertex& sourceVert);
+
   void addVertexOutwardBall(Vertex v, Vector<std::complex<double>>& distGradRHS);
 };
 

src/surface/surface_centers.cpp

@@ -10,37 +10,39 @@ namespace geometrycentral {
 namespace surface {
 
 
-SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const std::vector<Vertex>& vertexPts, int p) {
+SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom,
+                        const std::vector<Vertex>& vertexPts, int p, LogMapStrategy strategy) {
   VertexData<double> 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<Vertex>& vertexPts, int p) {
+                        const std::vector<Vertex>& vertexPts, int p, LogMapStrategy strategy) {
   VertexData<double> 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<double>& distribution, int p) {
+SurfacePoint findCenter(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom,
+                        const VertexData<double>& 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<double>& distribution, int p) {
+                        const VertexData<double>& 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<double, Vector2> {
     // Compute the current log map
     // Solve at the face point
-    VertexData<Vector2> logmap = solver.computeLogMap(aboutPoint);
+    VertexData<Vector2> logmap = solver.computeLogMap(aboutPoint, strategy);
 
     // Evaluate energy and update step
     double thisEnergy = 0.;