SUMMATION OVER THE DIFFRACTION SURFACE

In order to sum along the diffraction surface of a scattering point, the energy can be summed along circles in the time-slice representation of a 3D data volume. This process can be approximated by multiplying the 2D Fourier transform of each time slice by a cosine filter, corresponding to the radius of the intersection of the time slice with the diffraction surface.

The spatial Fourier domain is computationally attractive, since convolution in space reduces to a simple multiplication in the spatial frequency domain. The first task is then to perform a 2D Fourier transform of all the time slices. For that I used the CMSSL library on the Connection Machine, which includes an optimized multidimensional FFT subroutine. Figure shows the layout of the transformed time slices and their filters. The 3D data cube has now one serial (temporal) dimension and two parallel axes, namely k_x and k_y (spatial frequencies). Each spatial frequency in a time slice is handled by a separate processor, and all the spatial frequencies are calculated independently from each other. For each spatial frequency component the values for all times and three filter coefficients are stored in local memory. No interprocessor communication is necessary to migrate each frequency component. All filtered time-slice components that are summed, come directly from each processor's local memory. Thus the spatial frequency domain maps perfectly onto the Connection Machine layout. The summation over the diffraction surface proceedes along the serial time axis, an is carried out for all spatial frequency components simultaneously.

 

cmlayout
Figure 1 The spatial frequencies of individual time slices and filters are laid out on the Connection Machine in parallel. The time axis is kept as a serial dimension. No interprocessor communication is necessary to filter and sum individual spatial frequency components.

 

filter
Figure 2 Time-slice filters CSF( (n+1)*dr ) are generated by recursion as needed. This can be done for each spatial frequency independently. No interprocessor communication is necessary.

The pseudocode for migrating time slices is

$\bullet$

2D FFT of all time slices to be summed

$\bullet$

For each increment in diffraction radius

$\lbrace$

Generate time-slice filter CSF

Interpolate time slice

Filter time slice

Sum into migrated slice

$\rbrace$

$\bullet$

Inverse 2D FFT of migrated time slices

The most efficient way to obtain the filter values for each slice that is summed, is to apply the trigonometric recursion

$$CSF_{n+1} = 2~CSF_n ~CSF_1 - CSF_{n-1},\qquad\qquad\qquad n\gt 1 ,$$ (3)

where the subscript is a multiplication factor in the argument of CSF, such that $CSF_i~=~CSF(~ i~DR~2\pi\sqrt{k_x^2 + k_y^2}~)$, with the initialization CSF₀=1 and $CSF_1=CSF(~DR~2\pi\sqrt{k_x^2 + k_y^2}~)$.Thus each processor can calculate the new filter coefficient for each slice which is summed, with just three-floating point operations. Starting the summation at the vertex of the diffraction surface is advantageous in that the migration aperture does not have to be fixed at the beginning of the process. One can use this for an interactive session where the migration aperture can be chosen by inspection. The algorithm, as implemented now, is running typically at a speed of 190 Mflops on a 4k processor machine. It seems to me fast enough to be used as a tool for a very simple interactive 3D migration velocity analysis.