A fast algorithm for 3D azimuthally anisotropic velocity scan

Basic formulation

The right-hand side of equation 5 is a quotient of two (discrete) generalized Radon transforms (Beylkin, 1984). They can be expressed in a unified way as (to simplify the notation, we write $p=W_{\cos}$ , $q=W_{\sin}$ in this and next subsections)

$$(Rg)(\tau,p,q)=\sum_{x,y}g(\sqrt{\tau^2+p(x^2-y^2)+2qxy},x,y),$$ (7)

where $g$ is $d$ or some composite function of $d$ .

To construct the fast algorithm, we first rewrite equation 7 in the frequency domain as

$$(Rg)(\tau,p,q)=\sum_{f,x,y} e^{2\pi i f \sqrt{\tau^2+p(x^2-y^2)+2qxy}}\hat{g}(f,x,y),$$ (8)

where $f$ is frequency and $\hat{g}(f,x,y)$ is the Fourier transform of $g(t,x,y)$ in time. We next perform a linear transformation to map all discrete points in $(f,x,y)$ and $(\tau,p,q)$ domains to points in the unit cube $[0,1]^3$ ; i.e., a point $(f,x,y)\in[f_$min$,f_$max$]\times[x_$min$,x_$max$]\times [y_$min$,y_$max$]$ is mapped to $\mathbf{k}=(k_1,k_2,k_3)$ $\in[0,1]\times[0,1]\times[0,1]=K$ via

$$f=(f_ -f_ )k_1+f_ ,$$

$$x=(x_ -x_ )k_2+x_ ,$$

$$y=(y_ -y_ )k_3+y_ ;$$

a point $(\tau,p,q)\in[\tau_$min$,\tau_$max$]\times[p_$min$,p_$max$]\times[q_$min$,q_$max$]$ is mapped to $\mathbf{x}=(x_1,x_2,x_3)\in[0,1]\times[0,1]\times[0,1]=X$ via

$$\tau=(\tau_ -\tau_ )x_1+\tau_ ,$$

$$p=(p_ -p_ )x_2+p_ ,$$

$$q=(q_ -q_ )x_3+q_ .$$

If we define a phase function $\Phi(\mathbf{x},\mathbf{k})$ as

$$\Phi(\mathbf{x},\mathbf{k})= f \sqrt{\tau^2+p(x^2-y^2)+2qxy},$$ (9)

then equation 8 can be recast as

$$(Rg)(\mathbf{x})=\sum_{\mathbf{k}\in K} e^{ 2\pi i \Phi(\mathbf{x},\mathbf{k})}\hat{g}(\mathbf{k}), \quad \mathbf{x}\in X$$ (10)

(throughout the paper, $K$ and $X$ are used to denote either sets of discrete points or the cubic domains containing them).

A fast algorithm for 3D azimuthally anisotropic velocity scan