A fast butterfly algorithm for generalized Radon transforms

Algorithm

Assume $d(t,h)$ is a function in the data space. The hyperbolic Radon transform $R$ maps $d$ to a function $(Rd)(\tau,p)$ in the model space (Thorson and Claerbout, 1985) through

$$(Rd)(\tau,p)=\int d(\sqrt{\tau^2+p^2h^2},h)\,dh,$$ (1)

where $t$ is time, $h$ is offset, $\tau$ is intercept time, and $p$ is slowness. Fixing $(\tau,p)$ , hyperbola $t=\sqrt{\tau^2+p^2h^2}$ describes the traveltime for the event; hence integration along these curves can be used to identify different reflections.

Instead of approximating the integral in equation 1 directly, we reformulate it as a double integral,

$$(Rd)(\tau,p)=\iint\hat{d}(f,h) e^{ 2\pi i f\sqrt{\tau^2+p^2h^2} } \,df\,dh.$$ (2)

Here $f$ is the frequency, $\hat{d}(f,h)$ is the Fourier transform of $d(t,h)$ in $t$ . A simple discretization of equation 2 yields

$$(Rd)(\tau,p)=\sum_{f,h} e^{ 2\pi i f\sqrt{\tau^2+p^2h^2} }\hat{d}(f,h)$$ (3)

(the area element is omitted; the same symbols $f$ , $h$ , $\tau$ , and $p$ are used for both continuous and discrete variables). The reason that hyperbolic RT is harder to compute than linear RT ($t=\tau+ph$ ) or parabolic RT ( $t=\tau+ph^2$ ) should be clear from equation 3: product $f\tau$ in the phase cannot be decoupled from other terms.

To construct the fast algorithm, we first perform a linear transformation to map all discrete points in $(f,h)$ and $(\tau,p)$ domains to points in the unit square $[0,1]\times[0,1]$ ( $[a,b]\times[c,d]$ represents a 2D rectangular domain in the $xy$ -plane with $x\in[a,b]$ and $y\in[c,d]$ ), i.e., a point $(f,h)\in[f_$min$,f_$max$]\times[h_$min$,h_$max$]$ is mapped to $\mathbf{k}=(k_1,k_2)\in[0,1]\times[0,1]=K$ via

$$f=(f_ -f_ )k_1+f_ , \quad h=(h_ -h_ )k_2+h_ ;$$ (4)

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

$$\tau=(\tau_ -\tau_ )x_1+\tau_ , \quad p=(p_ -p_ )x_2+p_ .%\footnote{The notations $\mathbf{x}$ and $\mathbf{k}$ are adopted to be consistent with the Cand\\lq {e}s et al. paper referenced below.}$$ (5)

If we define input $g(\mathbf{k})=\hat{d}(f(k_1),h(k_2))$ , output $u(\mathbf{x})=(Rd)(\tau(x_1),p(x_2))$ , and phase function $\Phi(\mathbf{x},\mathbf{k})=f(k_1)\sqrt{\tau(x_1)^2+p(x_2)^2h(k_2)^2}$ , then equation 3 can be written as

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

(throughout the paper, $K$ and $X$ will either be used for sets of discrete points or square domains containing them; the meaning should be clear from the context). This form is the discrete version of an oscillatory integral of the type

$$u(\mathbf{x})=\int_K e^{ 2\pi i \Phi(\mathbf{x},\mathbf{k})}g(\mathbf{k}) \,d\mathbf{k},\quad \mathbf{x}\in X,$$ (7)

whose fast evaluation has been considered in Candès et al. (2009). Our method for computing the summation in equation 6 follows the FIO butterfly algorithm introduced there.

A fast butterfly algorithm for generalized Radon transforms