Arrival time and fronts

A fire starts at a single point in a landscape and spreads outwards. If the question is when it reaches one particular farmhouse, the cheapest path from the ignition point to that farmhouse answers it and the previous page’s machinery is enough. If the question is when it reaches every farmhouse, and where the fire line will be at four in the afternoon, then solving for one path at a time is the wrong shape of computation altogether, and enumerating all of them is hopeless.

The object that answers the second question is the arrival time itself. Fix the source and let $T(x)$ be the least cost of reaching $x$, for every point of the domain at once. This field can be computed directly rather than assembled out of individual paths. Away from the source and from shocks, its gradient lies on the unit sphere of the norm dual to $F$, the one that measures covectors as $F$ measures vectors, $F^*_x(\xi) = \sup\{\, \xi(v) : F_x(v) = 1 \,\}$:

The sublevel set $\{T \leq t\}$ is then everything reachable on a budget of $t$, which for a spreading process is the region burnt, infected or reacted by time $t$, so the fire line is a level set of $T$ and the entire history of the front is already contained in the one field. The optimal route to a destination can be recovered by following the characteristic backwards. Its forward tangent is the vector $v$ that attains the supremum in $F^*(x, \nabla T)$; backward tracing uses its negative. Except in the isotropic case, this is not simply the negative Euclidean gradient. At a shock several minimising routes may meet, so the characteristic is not unique there.

When $F$ is symmetric, the dual unit sphere containing $\nabla T$ is centred on the origin — an ellipsoid in the Riemannian case, and in the isotropic case a sphere of radius $1/c$, which gives the familiar equation $\lvert \nabla T \rvert = 1/c$. For a Randers metric in affine form $F = \sqrt{v^\top G v} + B^\top v$ the dual unit sphere is again an ellipsoid, displaced by the drift, and the equation becomes

For fixed $G$, the one-form $B$ translates the dual ellipsoid and breaks reversibility. It makes the arrival time grow differently with and against the drift, leaving the field asymmetric about its source even in a uniform medium. Setting $B = 0$ leaves $\nabla T^\top G^{-1} \nabla T = 1$, the symmetric anisotropic equation, and taking $G$ proportional to the identity on top of that gives the isotropic one. In the Zermelo parametrization of the previous page $(G, B)$ are algebraic functions of $(h, W, \lambda)$, so both fields are derived from the same navigation data used by the geodesic solver.

The numerics follow from where $B$ points. An update at a grid point has to draw its information from the directions the front actually arrives from, and once there is a drift, those are no longer the directions a symmetric stencil would consult. A stencil chosen without regard to the anisotropy produces a field that is smooth, plausible and wrong, which is a harder failure to notice than a solver that visibly diverges. Constructing stencils that respect the anisotropy has a substantial literature behind it, including Mirebeau’s work on adaptive grid discretisations. The forward sweep and implicit adjoint implemented here follow Gahtan, Shpund and Bronstein (2026), with the Randers conversion performed by finslerax.

EikonalSolver applies a fully anisotropic Godunov update on a dense two-dimensional Cartesian grid, alternating four directional Gauss–Seidel sweeps until convergence. MeshEikonalSolver uses local tangent-plane stencils on a triangulation, while VolumetricEikonalSolver uses a three-dimensional Cartesian stencil. All three consume Zermelo data and admit drift.

Each solve is differentiable with respect to its metric data at a memory cost independent of the number of sweeps. The converged field is a fixed point of the update operator, and the number of sweeps depends on the geometry rather than being known in advance, so differentiating by unrolling would store an unpredictable stack of intermediate fields. Instead a custom vector–Jacobian product runs an adjoint fixed-point iteration at the converged solution. The backward pass requires an additional fixed-point solve and $O(1)$ memory in the sweep count, which is what permits an entire front computation to sit inside a training loss and a metric to be fitted to observed arrival times directly.