# Heating Up

When experimenting with geometry and graphics algorithms, I rely pretty heavily on visual feedback so I often feel compelled 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++17 and compiled to WebAssembly via Emscripten.

It uses 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 out of the box with Emscripten.

- 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 unfamiliar 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.

Most algorithms for solving the Eikonal equation fall into the “advancing front” category, incrementally working their way out from $\partial \Omega$ in a concentric manner. The heat method, on the other hand, takes a fundamentally different approach 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 $$

For triangle meshes, the Laplace operator $\Delta$ is commonly represented 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 more efficient solver.

The second sparse linear system is the aforementioned Poisson problem. For this, we can use the same discrete Laplace operator and, again, rearrange to get a semidefinite matrix on the left.

Evaluating the right hand side of (2) (i.e. the integrated divergence of our approximated distance gradient) is fairly straightforward for triangle meshes - it just amounts to computing a (constant) gradient vector per face, scaling it, and then integrating the result over the boundary of each vertex dual cell (see paper section 3.2.1 for details).

It’s also worth noting that the solution to (2) is only unique up to translation since $C$ is *semi*definite. To recover the final approximation of geodesic distance, then, we 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.