Variational formulation of fast-marching algorithm

The fast-marching algorithm consists of two parts: minimum traveltime selection and a local traveltime updating scheme. The selection scheme is essentially based on Fermat's principle. The local traveltime updating scheme can also be derived from the same principle using a local linear interpolation, which provides the first-order accuracy.

 

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

For simplicity, let us focus on the 2-D case Consider a line segment with the end points A and B, as shown in Figure 6. 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\;.$$ (6)

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}\;,$$ (7)

where s_C is the local slowness, $\xi_0$ corresponds to the projection of C to the line AB (normalized by the length |AB|), 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\;,$$ (8)

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}}\;.$$ (9)

Finally, substituting the value of $\xi$ from (9) into equation (7) 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} + a\,t_A\,\cos{\beta} + b\,t_B\,\cos{\alpha}\;,$$ (10)

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

The above general procedure can be greatly simplified when applied to some regular grids, such as the rectangular grid or the tetragonal grid. This expression is even valid for unstructured grids. As pointed out by Fomel 1997, unstructured grids have some attractive computational advantages over the regular ones. Moreover, 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.