Keenan Crane
Geodesics in Heat: A New Approach to Computing Distance Based on Heat Flow
ACM Transactions on Graphics 2013
(to appear in the Communications of the ACM)
Keenan Crane Clarisse Weischedel Max Wardetzky
Caltech University of Göttingen University of Göttingen
We introduce the heat method for computing the 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 resulting systems can be prefactored once and subsequently solved in near-linear time. In practice, distance is updated an order of magnitude faster than with state-of-the-art methods, while maintaining a comparable level of accuracy. The method requires only standard differential operators and can hence be applied on a wide variety of domains (grids, triangle meshes, point clouds, etc.). We provide numerical evidence that the method converges to the exact distance in the limit of refinement; we also explore smoothed approximations of distance suitable for applications where greater regularity is required.
Fast Forward
Additional Notes
None of the current implementations of the heat method (listed on this page) fully exploit its potential for speed or accuracy. To substantially improve performance on multi-core or GPU-based systems, one need only link against a parallel sparse linear solver that can handle symmetric positive-definite systems. (One might also parallelize matrix construction.) To further improve accuracy, one can incorporate the very nice iterative scheme of Belyaev & Fayolle from their paper On Variational and PDE-Based Distance Function Approximations. Like the heat method, this scheme lends itself nicely to prefactorization (hence low amortized cost relative to fast marching/fast sweeping or window-based methods), and would be easy to implement on top of the existing reference implementation.

For a very cool application of the heat method, see Floraform.

It is also possible to use the heat method to solve the single- or multiple-source shortest path problem on general graphs.

@article{Crane:2013:GH, author = {Keenan Crane and Clarisse Weischedel and Max Wardetzky}, title = {{Geodesics in Heat: A New Approach to Computing Distance Based on Heat Flow}}, journal = {ACM Trans. Graph.}, volume = {32}, issue = {5}, year = {2013}, publisher = {ACM}, address = {New York, NY, USA} }
ANSI Coptimized reference implementation written by the authors; depends on an external linear solver like CHOLMOD or HSL MA87.
Python from PyCortex package for cortical processing (Alex Huth/Gallant Lab)
MATLAB from Gabriel Peyré's Numerical Tours
Mathematica answer by John Dunlop on StackExchange
C++ from StarLab
C++ from our SIGGRAPH course)
Cool interactive demo by Ruoqi He & Chia-Man Hung from Ecole Polytechnique.
This work was funded by a Google PhD Fellowship and a grant from the Fraunhofer Gesellschaft. Thanks to Michael Herrmann for inspiring discussions, and Ulrich Pinkall for discussions about Lemma 1. Meshes are provided courtesy of the Stanford Computer Graphics Laboratory, the AIM@Shape Repository, Luxology LLC, and Jotero GbR.
The heat method is a general principle that can be applied to any geometric data structure, as long as one knows how to take the gradient of a scalar function. It has been implemented on a variety of data structures including subdivision surfaces [de Goes et al 2016], voxel grids [Coeurjolly et al 2015], spline surfaces [Nguyen et al 2015], point clouds [Crane et al 2013], tetrahedral meshes [Belyaev & Fayolle 2015], and regular grids (using standard finite differences).
Outline of the heat method. (I) Heat is allowed to diffuse for a brief period of time (left). (II) The temperature gradient (center left) is normalized and negated to get a unit vector field (center right) pointing along geodesics. (III) A function whose gradient follows recovers the final distance (right).
The heat method can be applied directly to point clouds that lack connectivity information connectivity information.
Since the heat method is based on well-established discrete operators like the Laplacian, it is easy to adapt to a variety of geometric domains. Above: distance on a hippo composed of high-degree nonplanar (and sometimes nonconvex) polygonal faces.
Distance on an extremely poor triangulation with significant noise – note that small holes are essentially ignored. Also note good approximation of distance even along thin slivers in the nose.
Medial axis of the hiragana letter “a” extracted by thresholding second derivatives of the distance to the boundary. Left: fast marching. Right: heat method.
Tests of robustness. Left: our smoothed distance appears similar on meshes of different resolution. Right: even for meshes with severe noise (top) we recover a good approximation of the distance function on the original surface (bottom, visualized on noise-free mesh).
Distance to the boundary on a region in the plane (left) or a surface in space (right) is achieved by simply placing heat along the boundary curve. Note good recovery of the cut locus, i.e., points with more than one closest point on the boundary.
Left: geodesic furthest-point sampling found by repeatedly extracting adding the point of maximum distance to a set of sites. Right: geodesic Voronoi diagram found by iteratively keeping track of the minimum distance to a list of sites.
Thin features do not adversely affect the accuracy of the solution: here each pair of spheres in a long strand of “pearls” is connected by cutting out small holes (just individual vertex 1-rings) and gluing them together. The green source point (left) produces the expected distance function, even very far from the source on the opposite end of the strand (right).
In any method based on a finite element approximation, mesh quality will affect the quality of the solution. However, because the heat method is based on solving low-order elliptic equations (rather than high-order or hyperbolic equations), it often produces fewer numerical artifacts. Here, for instance, we highlight spurious extrema in the distance function (i.e., local maxima and minima) produced by the fast marching method (left), biharmonic distance (middle) and the heat method (right) on an acute Delaunay mesh (top) and a badly degenerate mesh (bottom). Inset figures show closeup view of isolines for the bottom figure.