A variational formulation of the fast marching eikonal solver

Variational principles on a grid

In this section, I derive a discrete traveltime computation procedure, based solely on Fermat's principle, and show that on a Cartesian rectangular grid it is precisely equivalent to the update formula (1) of the first-order eikonal solver.

triangle Figure 3. A geometrical scheme for the traveltime updating procedure in two dimensions.

For simplicity, let us focus on the two-dimensional case. Consider a line segment with the end points $A$ and $B$, as shown in Figure 3. Let $t_A$ and $t_B$ denote the traveltimes from a fixed distant source to points $A$ and $B$, respectively. Define a parameter $\xi$ such that $\xi=0$ at $A$, $\xi=1$ at $B$, and $\xi$ changes continuously on the line segment between $A$ and $B$. Then for each point of the segment, we can approximate the traveltime by the linear interpolation formula

$$t (\xi) = (1-\xi) t_A + \xi t_B\;.$$ (9)

Now let us consider an arbitrary point $C$ in the vicinity of $AB$. If we know that the ray from the source to $C$ passes through some point $\xi$ of the segment $AB$, then the total traveltime at $C$ is approximately

$$t_C = t (\xi) + s_C \sqrt{\vert AB\vert^2 (\xi-\xi_0)^2 + \rho_0^2}\;,$$ (10)

where $s_C$ is the local slowness, $\xi_0$ corresponds to the projection of $C$ to the line $AB$ (normalized by the length $\vert AB\vert$), and $\rho_0$ is the length of the normal from $C$ to $\xi_0$.

Fermat's principle states that the actual ray to $C$ corresponds to a local minimum of the traveltime with respect to raypath perturbations. According to our parameterization, it is sufficient to find a local extreme of $t_C$ with respect to the parameter $\xi$. Equating the $\xi$ derivative to zero, we arrive at the equation

$$t_B - t_A + \frac{s_C \vert AB\vert^2 (\xi-\xi_0)} {\sqrt{\vert AB\vert^2 (\xi-\xi_0)^2 + \rho_0^2}} = 0\;,$$ (11)

which has (as a quadratic equation) the explicit solution for $\xi$:

$$\xi = \xi_0 \pm \frac{\rho_0 (t_A - t_B)} {\vert AB\vert \sqrt{s_C^2 \vert AB\vert^2 - (t_A - t_B)^2}}\;.$$ (12)

Finally, substituting the value of $\xi$ from (12) into equation (10) and selecting the appropriate branch of the square root, we obtain the formula

$$c t_C = \rho_0 \sqrt{s_C^2 c^2 - (t_A - t_B)^2} + c t_A (1-\xi_0) + c t_B \xi_0 =$$

$$\rho_0 \sqrt{s_C^2 c^2 - (t_A - t_B)^2} + a t_A \cos{\beta} + b t_B \cos{\alpha}\;,$$ (13)

where $c = \vert AB\vert$, $a = \vert BC\vert$, $b = \vert AC\vert$, angle $\alpha$ corresponds to $\widehat{BAC}$, and angle $\beta$ corresponds to $\widehat{ABC}$ in the triangle $ABC$ (Figure 3).

square Figure 4. NR

A geometrical scheme for traveltime updating on a rectangular grid.

To see the connection of formula (13) with the eikonal difference equation (1), we need to consider the case of a rectangular computation cell with the edge $AB$ being a diagonal segment, as illustrated in Figure 4. In this case, $\cos{\alpha} = \frac{a}{c}$, $\cos{\beta} = \frac{b}{c}$, $\rho_0 = \frac{ab}{c}$, and formula (13) reduces to

$$t_C = \frac{ab \sqrt{s_C^2 (a^2 + b^2) - (t_A - t_B)^2} + a^2 t_A + b^2 t_B} {a^2+b^2}\;.$$ (14)

We can notice that (14) is precisely equivalent to the solution of the quadratic equation (13), which in our new notation takes the form

$$\left(\frac{t_C - t_A}{b}\right)^2 + \left(\frac{t_C - t_B}{a}\right)^2 = s_C^2\;.$$ (15)

What have we accomplished by this analysis? First, we have derived a local traveltime computation formula for an arbitrary grid. The derivation is based solely on Fermat's principle and a local linear interpolation, which provides the first-order accuracy. Combined with the fast marching evaluation order, which is also based on Fermat's principle, this procedure defines a complete algorithm of first-arrival traveltime calculation. On a rectangular grid, this algorithm is exactly equivalent to the fast marching method of Sethian (1996a) and Sethian and Popovici (1997). Second, the derivation provides a general principle, which can be applied to derive analogous algorithms for other eikonal-type (Hamilton-Jacobi) equations and their corresponding variational principles.

A variational formulation of the fast marching eikonal solver