next up previous [pdf]

Next: KOLMOGOROFF SPECTRAL FACTORIZATION Up: FILTERING ON A HELIX Previous: Examples of simple 2-D

Coding multidimensional convolution and deconvolution

Let us unroll the filter helix seen previously in Figure 2, and see what we have. Start from the idea that a 2-D filter is generally made from a cluster of values near one another in two dimensions similar to the Laplacian operator in the figure. We see that in the helical approach, a 2-D filter is a 1-D filter containing some long intervals of zeros. These intervals complete the length of a single 1-D seismogram.

Our program for 2-D convolution with a 1-D convolution program, could convolve with the somewhat long 1-D strip, but it is much more cost effective to ignore the many zeros, which is what we do. We do not multiply by the backside zeros, nor do we even store those zeros in memory. Whereas, an ordinary convolution program would do time shifting by a code line like iy=ix+lag, Module helicon ignores the many zero filter values on backside of the tube by using the code iy=ix+lag[ia] where a counter ia ranges over the nonzero filter coefficients. Before operator helicon is invoked, we need to prepare two lists, one list containing nonzero filter coefficients flt[ia], and the other list containing the corresponding lags lag[ia] measured to include multiple wraps around the helix. For example, the 2-D Laplace operator can be thought of as the 1-D filter:

\begin{displaymath}\begin{array}{\vert r\vert r\vert r\vert r\vert r\vert r\vert...
...& 1 & \\ \hline 1 & -4 & 1\\ \hline & 1 & \\ \hline \end{array}\end{displaymath} (6)

The first filter coefficient in equation (6) is $ +1$ as implicit to module helicon. To apply the Laplacian on a 1,000 $ \times$ 1,000 mesh requires the filter inputs:

                i    lag[i]   flt[i]
               ---   ------   -----
                0      999      1
                1     1000     -4
                2     1001      1
                3     2000      1

Here, we choose to use ``declaration of a type'', a modern computer language feature that is absent from Fortran 77. Fortran 77 has the built in complex arithmetic type. In module helix, we define a type filter, actually, a helix filter. After making this definition, it is used by many programs. The helix filter consists of three vectors, a real valued vector of filter coefficients, an integer valued vector of filter lags, and an optional vector that has logical values ``true'' for output locations that are not computed (either because of boundary conditions or because of missing inputs). The filter vectors are the size of the nonzero filter coefficients (excluding the leading 1.), while the logical vector is long and relates to the data size. The helix module allocates and frees memory for a helix filter.

api/c/helix.c
typedef struct sf_helixfilter {
    int     nh;
    float* flt;
    int*   lag;
    bool*  mis;
    float   h0;
} *sf_filter;

For those of you with no C experience, the ``->'' appearing in the helix module denotes a pointer. Fortran 77 has no pointers (or everything is a pointer). The behavior of pointers is somewhat different in each language. In C, pointer behavior is straightforward. In module helicon you see the expression aa->flt[ia]. It refers to the filter named aa. Any filter defined by the helix module contains three vectors, one of which is named flt. The second component of the flt vector in the aa filter is referred to as aa->flt[1] which in the foregoing example refers to the value 4.0 in the center of the Laplacian operator. For data sets like above with 1,000 points on the 1-axis, this value 4.0 occurs after 1,000 lags, thus aa->lag[1]=1000.

Our first convolution operator tcai1 was limited to one dimension and a particular choice of end conditions. With the helix and C pointers, the operator helicon is a multidimensional filter with considerable flexibility (because of the mis vector) to work around boundaries and missing data.

api/c/helicon.c
void sf_helicon_lop( bool adj, bool add, 
		     int nx, int ny, float* xx, float*yy) 
/*< linear operator >*/
{
    int ia, iy, ix;
    
    sf_copy_lop(adj, add, nx, nx, xx, yy);

    for (ia = 0; ia < aa->nh; ia++) {
	for (iy = aa->lag[ia]; iy < nx; iy++) {
	    if( aa->mis != NULL && aa->mis[iy]) continue;
	    ix = iy - aa->lag[ia];
	    if(adj) {
		xx[ix] += yy[iy] * aa->flt[ia];
	    } else {
		yy[iy] += xx[ix] * aa->flt[ia];
	    }
	}
    }
}
The code fragment aa->lag[ia] corresponds to b-1 in tcai1.

Operator helicon did the convolution job for Figure 1. As with tcai1, the adjoint of helix filtering is reversing the screw--filtering backwards.

The companion to convolution is deconvolution. The module polydiv is essentially the same as polydiv1 but here it was coded using our new filter type in module helix which simplifies our many future uses of convolution and deconvolution. Although convolution allows us to work around missing input values, deconvolution does not (any input affects all subsequent outputs), so polydiv never references aa->mis[ia].

api/c/polydiv.c
void sf_polydiv_lop( bool adj, bool add, 
		     int nx, int ny, float* xx, float*yy) 
/*< linear operator >*/
{
    int ia, iy, ix;
    
    sf_adjnull( adj, add, nx, ny, xx, yy);
    
    for (ix=0; ix < nx; ix++) tt[ix] = 0.;
    
    if (adj) {
	for (ix = nx-1; ix >= 0; ix-) {  
	    tt[ix] = yy[ix];
	    for (ia = 0; ia < aa->nh; ia++) {
		iy = ix + aa->lag[ia];
		if( iy >= ny) continue;
		tt[ix] -= aa->flt[ia] * tt[iy];
	    }
	}
	for (ix=0; ix < nx; ix++) xx[ix] += tt[ix];
    } else {
	for (iy = 0; iy < ny; iy++) { 
	    tt[iy] = xx[iy];
	    for (ia = 0; ia < aa->nh; ia++) {
		ix = iy - aa->lag[ia];
		if( ix < 0) continue;
		tt[iy] -= aa->flt[ia] * tt[ix];
	    }
	}
	for (iy=0; iy < ny; iy++) yy[iy] += tt[iy];
    }
}


next up previous [pdf]

Next: KOLMOGOROFF SPECTRAL FACTORIZATION Up: FILTERING ON A HELIX Previous: Examples of simple 2-D

2015-03-25