Source-linked AI summary
Geodesics in Heat
Keenan Crane, Clarisse Weischedel, Max Wardetzky
TL;DR
Computing geodesic distance requires avoiding inaccurate heat-kernel reconstruction and handling nonsmooth distance fields. The heat method separates gradient computation from distance recovery using sparse linear systems, achieving accuracy comparable to fast marching at lower computational cost while also supporting smoothed distances.
Problem
Inaccurate numerical reconstruction of the heat kernel can prevent Varadhan’s formula from producing correct distances, while geodesic distance is nonsmooth at cut-locus points.
Method
The heat method computes a distance gradient from a broad class of inputs, then recovers distance by solving sparse linear systems with a time step estimated as t = mh^2.
Results
The heat method exhibits errors of the same magnitude and rate of convergence as fast marching, while operating at lower computational cost.
Takeaways & Limitations
The method applies across regular and irregular grids, polygonal meshes, and point clouds, and can rapidly construct smoothed distance functions for applications requiring more regularity.
Takeaways & Limitations
The method’s accuracy could be further improved through alternative spatial discretizations, and the optimal time-step expression remains difficult to determine provably.
Abstract
from arXiv · showhide
We introduce the heat method for computing the shortest geodesic distance to a specified subset (e.g., point or curve) of a given domain. The heat method is robust, efficient, and simple to implement since it is based on solving a pair of standard linear elliptic problems. The method represents a significant breakthrough in the practical computation of distance on a wide variety of geometric domains, since the resulting linear systems can be prefactored once and subsequently solved in near-linear time. In practice, distance can be updated via the heat method an order of magnitude faster than with state-of-the-art methods while maintaining a comparable level of accuracy. We provide numerical evidence that the method converges to the exact geodesic distance in the limit of refinement; we also explore smoothed approximations of distance suitable for applications where more regularity is required.
1. INTRODUCTION
The heat method addresses the difficulty of recovering geodesic distance from numerically approximated heat kernels by separating gradient computation from distance recovery. It supports broad geometric discretizations and repeated distance queries through prefactorable sparse linear systems.
- Motivation: Varadhan’s formula recovers geodesic distance from the heat kernel but is sensitive to approximation and numerical error.The heat method is motivated by avoiding this sensitivity.
- Core idea: The heat method first computes the gradient of the distance field and then recovers the distance itself.
- Advantages: The method applies to regular and irregular grids, polygonal meshes, and unstructured point clouds.
- Advantages: Its sparse linear systems can be prefactored once and rapidly re-solved for repeated distance queries on a fixed domain.The paper identifies shape matching, path planning, and level set-based simulation as relevant applications.
2. RELATED WORK
Existing geodesic-distance methods commonly solve a nonlinear eikonal equation using propagation or relaxation schemes, but reuse, parallelization, implementation, and generalization remain challenging. The heat method instead offers a broadly applicable alternative based on heat flow and gradient information.
- Eikonal methods: The eikonal formulation is nonlinear and hyperbolic, so distance is typically computed with iterative schemes such as fast marching or fast sweeping.These methods are used on regular grids and triangulated surfaces, with indirect extensions to other domains.
- Eikonal methods: Fast marching and fast sweeping must recompute distances from scratch for different source subsets and present parallelization challenges.
- Polyhedral methods: Exact and approximate polyhedral-distance algorithms propagate information with priority queues, limiting reuse and complicating generalization beyond triangle meshes.
- Related heat-based methods: A related Schrödinger-equation method is limited to Euclidean domains and regular grids and requires arbitrary-precision arithmetic or multiple parameter choices for accurate results.
- Smooth distances: Existing smooth distances can be useful when differentiability is required but may substantially violate the unit-speed nature of geodesics.
3. THE HEAT METHOD
The heat method computes geodesic distance by diffusing heat, extracting a normalized gradient direction, and solving for a scalar potential whose gradient follows that direction. It supports point and piecewise-submanifold sources and approaches true distance as diffusion time vanishes.
- Algorithm: The method consists of heat diffusion, normalized negation of the temperature gradient, and a Poisson solve for the distance potential.The three steps are expressed as integrating heat flow, evaluating X = −∇u/|∇u|, and solving ∆φ = ∇·X.
- Algorithm: The resulting function φ approximates geodesic distance and approaches the true distance as t goes to zero.
- Sources: A Dirac delta source computes distance to a point, while a generalized Dirac over γ computes distance to a piecewise submanifold.
- Motivation: The method ignores heat-gradient magnitude and uses only its direction, reducing sensitivity to errors in the heat kernel’s decay rate.The distance potential is obtained by minimizing the discrepancy between ∇φ and the normalized field X.
3.1 Time Discretization
For time discretization, the heat method uses a single backward Euler step to solve a linear equation over the domain. Gradient normalization makes the resulting distance more accurate across diffusion times and supports smoothed distance construction.
- Time discretization: A single backward Euler step discretizes heat flow for fixed t and leads to a linear solve over the entire domain.Backward Euler also provides a maximum principle that prevents spurious oscillations.
- Robustness: Varadhan’s formula becomes inaccurate with approximate heat kernels, whereas gradient normalization produces more accurate solutions with evenly spaced isolines.
- Elliptic formulation: For a point source, the elliptic boundary-value solution is equal to the heat-flow solution up to a multiplicative constant.The relationship preserves the gradient direction needed by the subsequent steps.
3.2 Spatial Discretization
The heat method supports discretizations across meshes, polygonal domains, and point clouds using discrete differential operators. Its accuracy and cost depend on intrinsic dimension and the time-step choice.
- Spatial discretization: The heat method can be formulated on any domain with discrete gradient, divergence, and Laplace operators.This supports implementations across multiple geometric discretizations.
- Simplicial meshes: Standard Laplacian discretization on triangulated surfaces uses vertex areas and cotangent weights.The operator is represented as L = A^-1L_C, with A storing vertex areas and L_C the cotangent operator.
- Polygonal surfaces: Polygonal meshes use a polygonal Laplacian and approximate facewise gradient magnitudes to normalize the heat-flow vector field.The resulting operators are applied in the usual three heat-method steps.
- Point clouds: Point clouds without connectivity use a symmetric point-cloud Laplacian, least-squares gradients, and a final linear solve for the scalar potential.This produces distances even when only point samples are available.
- Point clouds: The computational cost depends primarily on intrinsic dimension n rather than ambient grid dimension m, which matters when m is significantly larger than n.The authors chose a discretization that was simple to implement in any dimension.
- Choice of time step: t = mh^2 provides a practical time-step estimate, with m = 1 yielding near-optimal accuracy across many irregular triangulations.The paper uses this time step uniformly except when explicitly seeking smoothed distance.
3.3 Smoothed Distance
Geodesic distance is nonsmooth at cut-locus points, motivating smooth approximations. The heat method obtains such approximations by using large integration times while retaining geodesic-like structure more effectively than several existing smooth distances.
- Motivation: Geodesic distance develops sharp cusps at cut-locus points where shortest paths are not unique.This nonsmoothness creates numerical difficulty for derivative-dependent applications and can be aesthetically undesirable.
- Existing smooth distances: Existing smooth distances provide regularity but can poorly approximate true geodesic distance, as indicated by uneven isoline spacing.They may also violate the unit-speed nature of geodesics.
- Comparison: The smoothed distance maintains isotropic circular profiles near the source, unlike biharmonic distance, whose level lines tend to be elliptical.Both smooth approximations mitigate the cusps visible in exact distance.
- Heat-based smoothing: Large values of t rapidly produce smoothed geodesic distances without changing computational cost.Normalization keeps isolines evenly spaced for any t, and the result is isometrically but not conformally invariant.
- Relation to other distances: Commute-time and biharmonic distances can be expressed using harmonic and biharmonic Green’s functions related to time-discrete heat flow.On constant-curvature manifolds, the relevant diagonal Green’s-function terms are constant.
- Limitation: For variable curvature, the Green’s functions provide only approximations of the corresponding distance functions.Finite t also regularizes the identity operator and prevents a logarithmic singularity.
3.4 Boundary Conditions
Boundary conditions influence smoothed geodesic distance and can be used to shape path-planning behavior. Averaged conditions aim to reduce boundary-shape effects, while Neumann and Dirichlet conditions favor different boundary interactions.
- Boundary effects: For exact distance, vanishing Neumann or Dirichlet conditions produce the same smooth-limit behavior.For smoothed geodesic distance at large t, however, boundary conditions alter the result.
- Boundary effects: Neumann, Dirichlet, and averaged conditions produce distinct smoothed-distance behaviors, with averaged conditions mimicking the corresponding boundary-free surface.Figure 11 compares these three conditions.
- Path planning: Neumann conditions encourage boundary adhesion, whereas Dirichlet conditions encourage boundary avoidance in path planning.The integration time also controls trajectories: small t yields straight-line paths, while large t yields more natural trajectories.
4. COMPARISON
The heat method matches fast marching in error magnitude and convergence while requiring lower computational cost, and remains effective across challenging meshes and discretizations.
- Accuracy and convergence: Both fast marching and the heat method exhibit apparent linear convergence, while exact polyhedral distance exhibits only quadratic approximation behavior.The comparison is made against mean edge length on triangulated surfaces.
- Computational comparison: Prefactorable sparse linear systems give the heat method lower computational cost than fast marching while retaining comparable error magnitude and convergence rate.Sparse systems can be prefactored and rapidly re-solved for repeated computations.
- Accuracy and convergence: Less than 1% error separates the heat method and fast marching from exact geodesic distance in the visual comparison using default parameters.The figure reports comparable accuracy for both methods.
- Applications: The method supports downstream geometric analyses including medial-axis extraction from second derivatives of distance.The comparison includes fast marching and the heat method on a hiragana letter.
- Robustness: The heat method recovers useful distances on poorly discretized, noisy, and extremely poor triangulations, including small holes and thin slivers.Smoothed distance remains similar across resolutions and approximates distance on noisy surfaces.
5. CONCLUSION
The paper concludes that the heat method is simple and general, while identifying discretization, parameter selection, and broader extensions as open directions.
- Conclusion: The heat method is simple, general, and easy to incorporate into a broad class of algorithms.The conclusion frames broad incorporability as a central practical property.
- Future work: Alternative discretizations and improved choices of the optimal t value remain areas for further investigation.The existing estimate typically already outperforms fast marching in mean error.
- Future work: Future extensions include transformations for broader Hamilton-Jacobi equations and weighted distance computation by rescaling source data.These are presented as open questions or possible directions.
A. A VARADHAN FORMULA FOR GRAPHS
The graph analogue shows that a small-time heat solution encodes combinatorial distance through the logarithmic asymptotics of its entries, with numerical experiments supporting the analysis.
- A. A VARADHAN FORMULA FOR GRAPHS: Lemma 1 considers a symmetric-matrix linear system on a graph with a source delta and positive time parameter t.The graph is induced by the matrix sparsity pattern.
- A. A VARADHAN FORMULA FOR GRAPHS: The first nonzero Neumann-series term corresponds to the number of edges between a vertex and the source.For a vertex n edges away, later terms become negligible as t approaches zero.
- A. A VARADHAN FORMULA FOR GRAPHS: As t goes to zero, log(ut)v / log t approaches graph distance, defined as the number of edges from source u to vertex v.The asymptotic is dominated by the first nonzero term.
- A. A VARADHAN FORMULA FOR GRAPHS: Numerical experiments on regular grids agree with the asymptotic analysis, showing the solution approaches combinatorial distance as t decreases.Figure 18 reports this behavior directly.
- Numerical validation: Distance-function approximations on spheres and tori are evaluated through convergence plots, while metric-property violations vanish under refinement.The tests examine absolute distance error, symmetry, and triangle inequality.