fast marching in spherical coordinates

In spherical, or polar, coordinates, we are effectively propagating a plane wave on a regular grid. For homogeneous media, this plane wave will have a front that is always parallel to the $\theta$-$\phi$ plane. As a result, traveltime calculation using equation (4) in spherical coordinates is always exact in homogeneous media.

 

march-pol Figure 4 The steps taken to implement the fast marching in polar coordinates. This implementation resembles the Cartesian coordinate one shown in Figure 2, but with a different grid orientation. Black circles imply computed traveltimes that are set because of their minimum traveltime value along the front. Gray circles constitute the front, and their values are ordered in the heap array from minimum to maximum. The minimum is always extracted first from the heap array at each step, its traveltime is set (given a black circle), and all surrounding grid points that are not set are computed and put into the heap array.

Figure 4 is similar to Figure 2, but using the polar coordinate system. The source is computed initially and set to zero for all angles $\theta$. When stretched on a regular grid, all points at the surface r=0 are set to zero. These points are inserted in the wavefront array and sorted from minimum traveltime to maximum. In the case of the source grid point, the sorting step is unnecessary because all traveltimes are equal to zero. The minimum is then extracted, and the traveltimes for neighboring grid points are computed using the following relation:

$$\begin{eqnarray} max(D_{ijk}^{-r}\,t,0)^2 & + & min(D_{ijk}^{+r}\,t,0)^2+ \nonum... ...ijk}^{-\phi}\,t,0)^2 & + & min(D_{ijk}^{+\phi}\,t,0)^2=s^2_{ijk}, \end{eqnarray}$$ (5)

where D_ijk^-r is the derivative of traveltime with respect to r at grid point i, j, k, given by

$$D_{ijk}^{-r}\,t = \frac{t_{i,j,k} - t_{i-1,j,k}}{\Delta r},$$

and

$$D_{ijk}^{+r}\,t = \frac{t_{i+1,j,k} - t_{i,j,k}}{\Delta r}.$$

The $D_{ijk}^{\theta}$ and $D_{ijk}^{\phi}$ derivatives are slightly different, given by

$$D_{ijk}^{-\theta}\,t = \frac{t_{i,j,k} - t_{i,j-1,k}}{r \Delta \theta},$$

$$D_{ijk}^{+\theta}\,t = \frac{t_{i,j+1,k} - t_{i,j,k}}{r \Delta \theta},$$

$$D_{ijk}^{-\phi}\,t = \frac{t_{i,j,k} - t_{i,j-1,k}}{r \sin \theta \Delta \phi},$$

and

$$D_{ijk}^{+\phi}\,t = \frac{t_{i,j+1,k} - t_{i,j,k}}{r \sin \theta \Delta \phi}.$$

Unlike the implementation in Cartesian coordinate, the heap array size in polar coordinates tends to be stable; as we extract a grid point, we usually insert another. In Cartesian coordinates, the heap array can become very large, especially in 3D media, as the wavefront expands. Figure 5 shows the size of the heap array for implementing the fast marching method on the model in Figure 7 using polar coordinates (solid curve) and Cartesian coordinates (dashed curve). Clearly, the heap array size in polar coordinate is practically the same through most of the computation, while the Cartesian version grows dramatically and then drops. The case becomes worse in 3-D media. Even with the polar coordinate system, the majority of the errors occur when the wavefront is at a 45 degree angle to the grid configuration. The time that the wavefront spends traveling diagonally along the polar coordinates is small, even in complex media, relative to that spent in the Cartesian coordinate system. Also, the wavefront curvature, another source of error in the scheme, though physically present because of a point source, does not computationally exist at the start of the wave along the spherical grid. As a result, the errors that accumulate at the start of wavefront propagation in the Cartesian coordinates are minimized in the spherical coordinate implementation.

 

count
Figure 5 The wavefront heap array size as a function of the number of iterations for traveltimes calculated in the model shown in Figure 7. The black curve corresponds to the fast marching method in polar coordinates, while the gray curve corresponds to the method in Cartesian coordinates.