Teaching a Mesh Which Way to Swim: Intrinsic Triangulations and the Anisotropic Laplacian

Introduction

Swimming with a current feels effortless; swimming across it is exhausting — even over the same distance in meters. That asymmetry isn’t a fact about water, it’s a fact about which notion of “distance” you’re using. This week’s project builds exactly that idea into a triangle mesh: a metric that makes some directions cheaper to move through than others. Along the way we needed two core tools from discrete differential geometry — intrinsic triangulations and the anisotropic Laplacian — and had to make both of them survive a real 3D shark scan, not just a hand-drawn toy mesh.

Background

A quick glossary before the derivations:

Mesh (V, F): vertex positions and the triangles connecting them.
Halfedge: a directed half of an edge, belonging to one triangle; its twin is the matching half from the triangle on the other side. This lets code ask “what’s across this edge?” in constant time.
Tensor field / metric S: a small symmetric, positive-definite matrix per triangle (2×2 in 2D, 3×3 on a 3D surface) describing how expensive movement is in each direction. S = I recovers ordinary Euclidean distance.
Intrinsic triangulation: a different triangulation of the same surface, changing only connectivity — never vertex positions.
Edge flip: the only move allowed — swap the diagonal of the quad formed by two triangles sharing an edge.
Delaunay triangulation: every edge satisfies the empty-circumcircle property; avoids thin, needle-like triangles.
Cotangent weight (wij): built from the two angles opposite an edge; decides both whether to flip it and what enters the Laplacian.
Harmonic field: a function solving $Lf=0$ with fixed boundary values — “ink diffusing smoothly” between two fixed points.

Method: Intrinsic Triangulations

An intrinsic triangulation keeps a mesh’s lengths and angles — what an ant walking on the surface could measure — while reconnecting vertices via edge flips:

  i                 i
 /|\               / \
/ | \             /   \
m |  k   ====>   m --- k
\ | /             \   /
 \|/               \ /
  j                 j

This matters because it can fix badly-shaped triangles without adding a single vertex — useful whenever an algorithm only needs lengths and angles, not raw coordinates.

The special case is the Delaunay triangulation, equivalent to:
$w_{ij} = \tfrac12(\cot\alpha+\cot\beta) \geq 0 \iff \alpha+\beta \leq 180°$

The proof is short: $\cot\alpha+\cot\beta = \frac{\sin(\alpha+\beta)}{\sin\alpha\sin\beta}$ , and since $\alpha,\beta\in(0°,180°)$ the denominator is always positive — so the sign of wij depends only on sin(α+β). Empty-circumcircle and angle-sum are the same condition. In code, this is a queue: flip any edge with negative weight, re-check its neighbors, repeat until none remain.

Method: The Anisotropic Laplacian

To make distance direction-dependent, attach a metric $M$ M to each triangle: $\|v\|_M=\sqrt{v^\top Mv}$ . M must be symmetric (for real eigenvalues/perpendicular axes) and positive-definite (so no length comes out negative or imaginary).

The pleasant surprise: the Delaunay proof above never used Euclidean angles specifically — only that they’re triangle angles in (0°,180°). Measuring α,β with M instead, the same algebra holds, so the criterion is unchanged in form. On a fixed test quad, switching S from identity to diag(1,16) flips the verdict entirely (w: +0.72 → -1.37) — same vertices, different verdict on whether the edge should exist.

The same $w_{ij}$ assembles the Laplacian ( $L_{ij}=-w_{ij}$ , $L_{ii}=\sum_j w_{ij}$ ), which must satisfy: rows sum to zero, symmetric, constants in the kernel. When all $w_{ij}\geq0$ (i.e., Delaunay), $L$ also satisfies the maximum principle — solutions to $Lf=0$ extremize only on the boundary. That guarantee is the entire reason the flip exists before building anything on top of $L$ .

Implementation

The build order, including what broke:

Synthetic 2D test bed: HalfedgeMesh with flip_edge/flip_to_delaunay, validated on a hand-built non-Delaunay quad and a randomly-scrambled disc.

Toy 2D fish: solved a harmonic field (nose=0, tail=1) with the ordinary Laplacian, took its per-triangle gradient as the “body direction” $d$ , and built $S=\alpha\,dd^\top+\beta(I-dd^\top)$ , $\alpha<\beta$ — cheap along the body, expensive across it.

Real 3D shark scan — the actual porting work:

Cleanup: weld coincident vertices, drop degenerate (zero-area) faces, since these send cotangent weights to nan and corrupt connectivity.

Automatic nose/tail: SVD on centered vertices finds the principal axis; nose and tail are its two extremes.

A genuinely 3D gradient: no single “perpendicular” exists on a surface, so the normal $n$ and cross product replace 2D rotation: $\nabla f = \frac{1}{2A}\sum_i f_i(n\times e_i)$ , and the tensor grows to 3×3 with a transverse direction $t=n\times d$ .

Hardening the flip queue: index by (face, corner) instead of hashing mutable Halfedge objects; return None instead of nan for degenerate triangles; reject flips that would duplicate an existing edge.

Results

The Polyscope UI toggles between the original mesh and the re-triangulated one, with a live counter of non-Delaunay edges before/after and total flips performed. Recomputing the direction field after flipping shows the payoff: geometry never moved, but connectivity reorganized to better respect the metric — the field still tracks the shark’s body, now on better-conditioned triangles.

Author

0262205