Heating Up

When experimenting with geometry and graphics algorithms, I rely pretty heavily on visual feedback so I often feel the need to write little interactive prototypes. I’ve been looking for a low-effort way to share some of these on the web over the past couple years and I think I’ve finally landed on something manageable so it’s time to put this blog to use.

To get started, here’s a simple demo of the heat method being used to approximate geodesic distance on triangle meshes. It’s written in C++ and compiled to WebAssembly via Emscripten.

Launch in new tab

I’m using Sokol libraries for context creation, graphics, and input which made building for the web significantly less painful than expected. Other dependencies are listed below, all of which worked with Emscripten out of the box.

Eigen for sparse linear algebra
Dear ImGui for UI
stb and hapPLY for asset loading

The Heat Method

I’ll give a brief overview of the heat method here though it’s mostly for my own benefit. The original publication and accompanying presentation offer excellent explanations that are chock-full of helpful diagrams and visualizations so most readers would likely be better off going straight to the source.

Summary

In general, geodesic distance can be expressed as the solution to the Eikonal equation i.e. a scalar function $\phi$ over some domain $\Omega$ that’s zero within a given subset $\partial \Omega$ and whose gradient is unit length everywhere else.

$$ \begin{aligned} |\nabla \phi(x)| &= 1 \quad \forall x \in \Omega \\ \phi(x) &= 0 \quad \forall x \in \partial \Omega \end{aligned} $$

Most algorithms for solving the Eikonal equation resemble something like a breadth-first search i.e. they model an advancing front that incrementally works its way out from $\partial \Omega$. In contrast, the heat method frames the problem as a pair of elliptic PDEs based on a clever observation about the relationship between distance and short-time heat flow.

If we treat $\partial \Omega$ as a heat source and allow temperature to diffuse over $\Omega$ for some small duration $t$, we end up with a temperature distribution $u_t$ whose gradient is almost parallel to the gradient of geodesic distance $\nabla \phi$. We can, therefore, approximate $\nabla \phi$ by simply normalizing the negated temperature gradient since we know $|\nabla \phi| = 1$ by definition.

$$ \frac{-\nabla u_t}{|\nabla u_t|} \approx \nabla \phi $$

From this, we can recover an approximation of geodesic distance $\phi^\ast$ by solving Poisson’s equation where the right hand side comes from the standard definition of the Laplacian i.e. the divergence of the gradient.

$$ \Delta \phi^\ast = \nabla \cdot \frac{-\nabla u_t}{|\nabla u_t|} $$

Implementation

From an implementation perspective, the heat method essentially boils down to solving a pair of sparse linear systems. The first is a single backward Euler step of heat flow to get temperature on $\Omega$ at time $t$.

$$ u_t - t \Delta u_t = u_0 $$

$$ u_0(x) = \begin{cases} 1, & \text{if $x \in \partial \Omega$} \\ 0, & \text{otherwise} \end{cases} $$

For triangle meshes, the Laplace operator $\Delta$ is commonly discretized as the product of an inverse diagonal “lumped mass” matrix $M^{-1}$ and a sparse semidefinite “stiffness” matrix $C$ (see paper section 3.2.1 for details).

$$ \Delta = M^{-1} C $$

Using this discrete Laplace operator, we get the following linear system for the backward Euler step which can be rearranged to get a semidefinite matrix on the left, allowing us to use a faster solver.

$$ \tag{1} \begin{aligned} u_t - t M^{-1} C u_t &= u_0 \\ (I - t M^{-1} C) u_t &= u_0 \\ M^{-1}(M - t C) u_t &= u_0 \\ (M - t C) u_t &= M u_0 \end{aligned} $$

Solving (1) gives us the temperature at time $t$. From here, we can assemble the second sparse linear system representing Poisson’s equation. This makes use of the same discrete Laplace operator and, again, can be rearranged to get a semidefinite matrix on the left.

$$ \tag{2} \begin{aligned} M^{-1}C \phi^\ast &= \nabla \cdot \frac{-\nabla u_t}{|\nabla u_t|} \\ C \phi^\ast &= M \left(\nabla \cdot \frac{-\nabla u_t}{|\nabla u_t|} \right) \end{aligned} $$

Evaluating the right hand side of (2) turns out to be fairly straightforward in the case of triangle meshes. This term can be interpreted as the divergence of our approximated distance gradient integrated over the dual cell of each vertex. Taking advantage of Stokes theorem, we can get the same result by integrating flux over the boundary of each vertex dual cell which ends up being quite a simple thing to do since the gradient is constant per face (see paper section 3.2.1 for details).

The solution to (2) is only unique up to translation since $C$ is semidefinite. To recover the final approximation of geodesic distance, we can simply subtract off the average of the solution in $\partial \Omega$ i.e. where we expect geodesic distance to be zero.

So, to recap, implementing the heat method can be broken down into the following steps:

Solve a backward Euler step of heat flow (1) to get temperature at time $t$
Normalize the negated gradient of temperature to get a vector field that closely resembles the gradient of geodesic distance
Solve the Poisson equation (2) to recover geodesic distance

Notice how none of this assumes that $\Omega$ has been discretized in a particular way. The only assumption is that we have working definitions of the Laplacian, gradient, and divergence operators for the chosen discretization. We can, therefore, follow this exact same procedure to compute geodesic distance on any spatial discretization (e.g. triangle meshes, polygon meshes, tetrahedral meshes, point clouds) so long as these differential operators are available. Now we’re cooking.

2021-06-08