A variational formulation of the fast marching eikonal solver |
For a detailed description of level set methods, the reader is referred to Sethian's recently published book (Sethian, 1996b). More details on the fast marching method appear in articles by Sethian (1996a) and Sethian and Popovici (1997). This section serves as a brief introduction to the main bulk of the algorithm.
The key feature of the algorithm is a carefully selected order of traveltime evaluation. At each step of the algorithm, every grid point is marked as either Alive (already computed), NarrowBand (at the wavefront, pending evaluation), or FarAway (not touched yet). Initially, the source points are marked as Alive, and the traveltime at these points is set to zero. A continuous band of points around the source are marked as NarrowBand, and their traveltime values are computed analytically. All other points in the grid are marked as FarAway and have an ``infinitely large'' traveltime value.
An elementary step of the algorithm consists of the following moves:
An update procedure is based on an upwind first-order approximation to
the eikonal equation. In simple terms, the procedure starts with
selecting one or more (up to three) neighboring points around the
updated point. The traveltime values at the selected neighboring
points need to be smaller than the current value. After the selection,
one solves the quadratic equation
Except for the updating scheme (1), the fast marching algorithm bears a very close resemblance to the famous shortest path algorithm of Dijkstra (1959). It is important to point out that unlike Moser's method, which uses Dijkstra's algorithm directly (Moser, 1991), the fast marching approach does not construct the ray paths from predefined pieces, but dynamically updates traveltimes according to the first-order difference operator (1). As a result, the computational error of this method goes to zero with the decrease in the grid size in a linear fashion. The proof of validity of the method (omitted here) is also analogous to that of Dijkstra's algorithm (Sethian, 1996a,b). As in most of the shortest-path implementations, the computational cost of extracting the minimum point at each step of the algorithm is greatly reduced [from to operations] by maintaining a priority-queue structure (heap) for the NarrowBand points (Cormen et al., 1990).
Figure 1 shows an example application of the fast marching eikonal solver on the three-dimensional SEG/EAGE salt model. The computation is stable despite the large velocity contrasts in the model. The current implementation takes about 10 seconds for computing a 100x100x100 grid on one node of SGI Origin 200. Alkhalifah and Fomel (1997) discuss the differences between Cartesian and polar coordinate implementations.
salt
Figure 1. Constant-traveltime contours of the first-arrival traveltime, computed in the SEG/EAGE salt model. A point source is positioned inside the salt body. The top plot is a diagonal slice; the bottom plot, a depth slice. | |
---|---|
The difference equation (1) is a finite-difference approximation
to the continuous eikonal equation
A variational formulation of the fast marching eikonal solver |