The Heat Method

2 min read•1 page

Heat Equation Formulation:

1. Backward Euler Integration: Solves (M - t L_C) u = δ_x where L_C is cotangent Laplacian and M is lumped area. 2. Time Step t = h^2: Scaled to mean edge length h, balancing numerical stability with Varadhan's asymptotic limit. 3. Linear System: Solvable via a single sparse Cholesky factorization for millisecond queries.

python

1import numpy as np
2
3def integrate_heat_diffusion(L_cot, M_lumped, source_idx, t_step=0.05):
4    """
5    Step 1 of the Heat Method (Crane et al. 2013).
6    Integrate backward Euler heat diffusion: (M - t*L) * u = delta_source.
7    t_step is chosen as h^2 (mean edge length squared).
8    """
9    n = L_cot.shape[0]
10    rhs = np.zeros(n)
11    rhs[source_idx] = 1.0  # Dirac delta heat source
12    
13    # System matrix A = M - t * L_cot (positive definite)
14    A = np.diag(M_lumped) - t_step * L_cot
15    u = np.linalg.solve(A, rhs)
16    return u

The Heat Method: Heat Diffusion Step: Rather than running slow Dijkstra or Fast Marching on discrete graph edges, the Heat Method computes exact polyhedral geodesic distances by simulating short-time heat diffusion from a source vertex du/dt = Δu. Heat flows smoothly across triangle interiors independent of mesh tessellation.

Diffusion Time (t)

0.12

5 min read•2 pages

Vector Field Integration:

1. Face Gradients: Evaluate linear piecewise gradient ∇u across each triangle. 2. Unit Normalization: X = -∇u / ||∇u|| eliminates exponential decay and preserves direction. 3. Poisson Reconstruction: Solve symmetric positive semi-definite system L_C φ = div(X) with Dirichlet shift φ(x_0)=0.

python

1import numpy as np
2from scipy.sparse.linalg import spsolve
3
4def compute_heat_geodesics(vertices, faces, M, L_cot, source_idx, t):
5    """
6    The 3-step Heat Method for Geodesic Distance (Crane et al. 2013).
7    """
8    n_verts = len(vertices)
9    u0 = np.zeros(n_verts)
10    u0[source_idx] = 1.0
11
12    # Step 1: Integrate heat flow (M - t * L_cot) u = u0
13    heat_u = spsolve(M - t * L_cot, u0)
14
15    # Step 2: Evaluate face gradients and normalize X = -grad_u / ||grad_u||
16    div_X = np.zeros(n_verts)
17    for i, j, k in faces:
18        vi, vj, vk = vertices[i], vertices[j], vertices[k]
19        e_ij, e_jk, e_ki = vj - vi, vk - vj, vi - vk
20        normal = np.cross(e_ij, -e_ki)
21        dbl_area = np.linalg.norm(normal)
22        n_hat = normal / dbl_area
23
24        grad_u = (
25            heat_u[i] * np.cross(n_hat, e_jk) +
26            heat_u[j] * np.cross(n_hat, e_ki) +
27            heat_u[k] * np.cross(n_hat, e_ij)
28        ) / dbl_area
29        X = -grad_u / (np.linalg.norm(grad_u) + 1e-12)
30
31        # Integrated divergence per vertex using cotangent weights
32        cot_i = -np.dot(e_ij, e_ki) / dbl_area
33        cot_j = -np.dot(e_jk, e_ij) / dbl_area
34        cot_k = -np.dot(e_ki, e_jk) / dbl_area
35        div_X[i] += 0.5 * (cot_j * np.dot(e_ij, X) + cot_k * np.dot(-e_ki, X))
36        div_X[j] += 0.5 * (cot_k * np.dot(e_jk, X) + cot_i * np.dot(-e_ij, X))
37        div_X[k] += 0.5 * (cot_i * np.dot(e_ki, X) + cot_j * np.dot(-e_jk, X))
38
39    # Step 3: Solve Poisson equation L_cot * phi = div_X
40    phi = spsolve(L_cot, div_X)
41    return phi - phi.min()

The Heat Method: Vector Field & Poisson Solve: Geodesics satisfy the Eikonal equation‖∇ϕ‖ = 1. Varadhan's formula guarantees that as t → 0, -∇u / ‖∇u‖ aligns precisely with the geodesic gradients ∇ϕ. Solving the Poisson problem Δϕ = ∇ · X integrates this vector field back into exact geodesic distances ϕ.

Geodesic Isolines

6.00

Related Concepts

Computational Metrics

Volume, mass properties, and surface curvature.

Explore concept →

Area And Volume

Exact polyhedral surface area via cross-product magnitudes and enclosed volume via signed tetrahedral summation.

Explore concept →

Centroid And Moments

Computing shell and solid centroids, 3x3 inertia tensors, principal axes of inertia, and static balance.

Explore concept →