The Heat Method for Distance Computation
Communications of the ACM (2017)
We introduce the heat method for solving the single- or multiple-source shortest path problem on both flat and curved domains. A key insight is that distance computation can be split into two stages: first find the direction along which distance is increasing, then compute the distance itself. The heat method is robust, efficient, and simple to implement since it is based on solving a pair of standard sparse linear systems. These systems can be factored once and subsequently solved in near-linear time, substantially reducing amortized cost. Real-world performance is an order of magnitude faster than state-of-the-art methods, while maintaining a comparable level of accuracy. The method can be applied in any dimension, and on any domain that admits a gradient and inner product—including regular grids, triangle meshes, and point clouds. Numerical evidence indicates 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 desired.
Fast Forward
Additional Notes
In the case of triangle meshes, any modern implementation of the heat method should really make use of the intrinsic Laplace operator, which adds negligible computational overhead, and has identical inputs/outputs, but is much (much) more robust on low-quality meshes. The GeometryCentral, CGAL, and libigl implementations listed below provide this option. The GeometryCentral version also works on general meshes, including non-manifold ones, by adopting the tufted Laplacian.
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.
Presentation
Geodesics on point clouds
(courtesy Oleg Memukhin)
@article{Crane:2017:HMD,
author = {Crane, Keenan and Weischedel, Clarisse and Wardetzky, Max},
title = {The Heat Method for Distance Computation},
journal = {Commun. ACM},
issue_date = {November 2017},
volume = {60},
number = {11},
month = oct,
year = {2017},
issn = {0001-0782},
pages = {90--99},
numpages = {10},
url = {http://doi.acm.org/10.1145/3131280},
doi = {10.1145/3131280},
acmid = {3131280},
publisher = {ACM},
address = {New York, NY, USA},
}
Source
| JavaScript— | interactive version running in-browser (requires WebGL support). |
| geometry-central— | robust C++ implementation (with optional Python bindings) that works on low-quality and non-manifold meshes. |
| point clouds— | implementation in geometry-central (C++ with optional Python bindings). Computes geodesic distance directly on point clouds, as well as other quantities like scalar extension, vector extension, and the log map via the vector heat method. |
| ANSI C— | optimized reference implementation written by the authors; depends on an external linear solver like CHOLMOD or HSL MA87. |
| CGAL— | C++ implementation in the Computational Geometry Algorithms Library (CGAL). This version uses the intrinsic Laplacian of Bobenko & Springborn to improve accuracy on low-quality meshes. |
| libigl— | header-only C++ implementation in the libigl geometry processing library. Like CGAL, this version (optionally) uses an intrinsic Delaunay triangulation to improve accuracy and robustness. |
| Houdini— | the heat method is built in to SideFX Software's Houdini package, as a “Heat Geodesic” geometry node. |
| 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 |
| Mathematica— | written by Henrik Schumacher. |
| C++— | from StarLab |
| C++— | from our SIGGRAPH course |
| Unity— [Mac] [Windows] [Linux] | Cool interactive demo by Ruoqi He & Chia-Man Hung from Ecole Polytechnique. |
| Grasshopper/Rhino— | Implementation by Daniel Piker, solved using a simple iterative scheme. |
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.
Figures
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.
An easy way to improve quality on low-quality meshes is to use the intrinsic Laplace-Beltrami operator of Bobenko & Springborn, then simply copy the results from the vertices of the intrinsic triangulation back to the vertices of the original triangulation. This strategy does not change the mesh, but more robustly deals with bad elements (as shown above). It is used by default in the CGAL version of the heat method.