next up previous [pdf]

Next: INVERSE FILTERS AND OTHER Up: The helical coordinate Previous: Filtering mammograms

SUBSCRIPTING A MULTIDIMENSIONAL HELIX

Basic utilities transform back and forth between multidimensional matrix coordinates and helix coordinates. The essential module used repeatedly in applications later in this book is createhelix. We begin here from its intricate underpinnings.

Fortran77 has a concept of a multidimensional array being equivalent to a 1-dimensional array. Given that the hypercube specification nd=(n1,n2,n3,...) defines the storage dimension of a data array, we can refer to a data element as either dd(i1,i2,i3,...) or dd( i1 +n1*(i2-1) +n1*n2*(i3-1) +...). The helix says to refer to the multidimensional data by its equivalent 1-dimensional index (sometimes called its vector subscript or linear subscript).

The filter, however, is a much more complicated story than the data: First, we require all filters to be causal. In other words, the Laplacian does not fit very well, because it is intrinsically noncausal. If you really want noncausal filters, you need to provide your own time shifts outside the tools supplied here. Second, a filter is usually a small hypercube, say aa(a1,a2,a3,...) and would often be stored as such. For the helix we must store it in a special 1-dimensional form. Either way, the numbers na= (a1,a2,a3,...) specify the dimension of the hypercube. In cube form, the entire cube could be indexed multidimensionally as aa(i1,i2,...) or it could be indexed 1-dimensionally as aa(ia,1,1,...) or sometimes aa[ia] by letting ia cover a large range. When a filter cube is stored in its normal ``tightly packed'' form, the formula for computing its 1-dimensional index ia is:

  ia = i1 +a1*i2 +a1*a2*i3 + ...
When the filter cube is stored in an array with the same dimensions as the data, data[n3][n2][n1], the formula for ia is:
  ia = i1 +n1*i2 +n1*n2*i3 + ...

The following module decart contains two subroutines that explicitly provide us the transformations between the linear index i and the multidimensional indices ii= (i1,i2,...). The two subroutines have the logical names cart2line and line2cart.

api/c/decart.c
void sf_line2cart(int dim         /* number of dimensions */, 
		  const int* nn /* box size [dim] */, 
		  int i         /* line coordinate */, 
		  int* ii       /* cartesian coordinates [dim] */)
/*< Convert line to Cartesian >*/
{
    int axis;
 
    for (axis = 0; axis < dim; axis++) {
	ii[axis] = i%nn[axis];
	i /= nn[axis];
    }
}

int sf_cart2line(int dim         /* number of dimensions */, 
		 const int* nn /* box size [dim] */, 
		 const int* ii /* cartesian coordinates [dim] */) 
/*< Convert Cartesian to line >*/
{
    int i, axis;

    if (dim < 1) return 0;

    i = ii[dim-1];
    for (axis = dim-2; axis >= 0; axis-) {
	i = i*nn[axis] + ii[axis];
    }
    return i;
}

The Fortran linear index is closely related to the helix. There is one major difference, however, and that is the origin of the coordinates. To convert from the linear index to the helix lag coordinate, we need to subtract the Fortran linear index of the ``1.0'' usually taken at center= (1+a1/2, 1+a2/2, ..., 1). (On the last dimension, there is no shift, because nobody stores the volume of zero values that would occur before the 1.0.) The decart module fails for negative subscripts. Thus, we need to be careful to avoid thinking of the filter's 1.0 (shown in Figure 8) as the origin of the multidimensional coordinate system although the 1.0 is the origin in the 1-dimensional coordinate system.

Even in 1-D (see the matrix in equation ([*])), to define a filter operator, we need to know not only filter coefficients and a filter length, but we also need to know the data length. To define a multidimensional filter using the helix idea, besides the properties intrinsic to the filter, also the circumference of the helix, i.e., the length on the 1-axis of the data's hypercube as well as the other dimensions nd=(n1,n2,...) of the data's hypercube.

Thinking about convolution on the helix, it is natural to think about the filter and data being stored in the same way, that is, by reference to the data size. Such storeage would waste so much space, however, that our helix filter module helix instead stores the filter coefficients in one vector and the lags in another. The i-th coefficient value of the filter goes in aa->flt[i] and the i-th lag ia[i] goes in aa->lag[i]. The lags are the same as the Fortran linear index except for the overall shift of the 1.0 of a cube of data dimension nd. Our module for convolution on a helix, helicon. has already an implicit ``1.0'' at the filter's zero lag, so we do not store it. (It is an error to do so.)

Module createhelix allocates memory for a helix filter and builds filter lags along the helix from the hypercube description. The hypercube description is not the literal cube seen in Figure 8 but some integers specifying that cube: the data cube dimensions nd, likewise the filter cube dimensions na, the parameter center identifying the location of the filter's ``1.0'', and a gap parameter used in a later chapter. To find the lag table, module createhelix first finds the Fortran linear index of the center point on the filter hypercube. Everything before that has negative lag on the helix and can be ignored. (Likewise, in a later chapter, we see a gap parameter that effectively sets even more filter coefficients to zero so those extra lags can also be ignored.) Then, it sweeps from the center point over the rest of the filter hypercube calculating for a data-sized cube nd, the Fortran linear index of each filter element.

user/gee/createhelix.c
sf_filter createhelix(int ndim    /* number of dimensions */, 
		      int* nd     /* data size [ndim] */, 
		      int* center /* filter center [ndim] */, 
		      int* gap    /* filter gap [ndim] */, 
		      int* na     /* filter size [ndim] */)
/*< allocate and output a helix filter >*/ 
{
    sf_filter aa;
    int ii[SF_MAX_DIM], na123, ia, nh, lag0a,lag0d, *lag, i;
    bool skip;

    for (na123 = 1, i=0; i < ndim; i++) na123 *= na[i];
    lag = (int*) alloca(na123*sizeof(int));

    /* index pointing to the "1.0" */
    lag0a = sf_cart2line (ndim, na, center); 

    nh=0;
    /* loop over linear index. */
    for (ia = 1+lag0a; ia < na123; ia++) { 
		sf_line2cart(ndim, na, ia, ii);
	
		skip = false;
		for (i=0; i < ndim; i++) {
			if (ii[i] < gap[i]) {
				skip = true;
				break;
			}
		}
		if (skip) continue;
	
		lag[nh] = sf_cart2line(ndim, nd, ii);
		nh++;                        /* got another live one */
    }
    /* center shift for nd cube */
    lag0d = sf_cart2line(ndim,  nd, center); 
    aa = sf_allocatehelix(nh); /* nh becomes size of filter */

    for (ia=0; ia < nh; ia++) {
		aa->lag[ia] = lag[ia] - lag0d; 
		aa->flt[ia] = 0.;
    }

    return aa;
}
Near the end of the code you see the calculation of a parameter lag0d, which is the count of the number of zeros that a data-sized Fortran array would store in a filter cube preceding the filter's 1.0. We need to subtract this shift from the filter's Fortran linear index to get the lag on the helix.

A filter can be represented literally as a multidimensional cube like equation (9) shows us in two dimensions or like Figure 8 shows us in three dimensions. Unlike the helical form, in literal cube form, the zeros preceding the ``1.0'' are explicitly present, so lag0 needs to be added back in to get the Fortran subscript. To convert a helix filter aa to Fortran's multidimensional hypercube cube(n1,n2,...) is module boxfilter:

user/gee/boxfilter.c
void box (int dim            /* number of dimaneions */, 
	  const int *nd      /* data size [dim] */, 
	  const int *center  /* filter center [dim] */, 
	  const int *na      /* filter size [dim] */, 
	  const sf_filter aa /* input filter */, 
	  int nc             /* box size */, 
	  float* cube        /* output box [nc] */)
/*< box it >*/ 
{
    int ii[SF_MAX_DIM];
    int j, lag0a, lag0d, id, ia;

    for (ia=0; ia < nc; ia++) {
	cube[ia] = 0.;
    }
    lag0a = sf_cart2line(dim, na, center);  /* 1.0 in na. */
    cube[lag0a] = 1.;                       /* place it. */
    lag0d = sf_cart2line(dim, nd, center);  /* 1.0 in nd. */
    for (j=0; j < aa->nh; j++) { /* inspect the entire helix */
	id = aa->lag[j] + lag0d;
	sf_line2cart(dim, nd, id, ii);	/* ii = cartesian  */
	ia = sf_cart2line(dim, na, ii);	/* ia = linear in aa */
	cube[ia] = aa->flt[j];		/* copy the filter */
    }
}
The boxfilter module is normally used to display or manipulate a filter that was estimated in helical form (usually estimated by the least-squares method).

A reasonable arrangement for a small 3-D filter is na={5,3,2} and center={3,2,1}. Using these arguments, I used createhelix to create a filter. I set all the helix filter coefficients to 2. Then I used module boxfilter to put it in a convenient form for display. Finally, I printed it:

          0.000  0.000  0.000  0.000  0.000
          0.000  0.000  1.000  2.000  2.000
          2.000  2.000  2.000  2.000  2.000
          ---------------------------------
          2.000  2.000  2.000  2.000  2.000
          2.000  2.000  2.000  2.000  2.000
          2.000  2.000  2.000  2.000  2.000

Different data sets have different sizes. To convert a helix filter from one data size to another, we could drop the filter into a cube with module cube. Then, we could extract it with module unbox specifying any data set size we wish. Instead, we use module regrid prepared by Sergey Fomel which does the job without reference to an underlying filter cube. He explains his regrid module thus:

Imagine a filter being cut out of a piece of paper and glued on another paper, which is then rolled to form a helix.

We start by picking a random point (let us call it rand) in the cartesian grid and placing the filter so that its center (the leading 1.0) is on top of that point. rand should be larger than (or equal to) center and smaller than min (nold, nnew), otherwise the filter might stick outside the grid (our piece of paper.) rand=nold/2 will do (assuming the filter is small), although nothing should change if you replace nold/2 with a random integer array between center and nold - na.

The linear coordinate of rand is h0 on the old helix and h1 on the new helix. Recall that the helix lags aa->lag are relative to the center. Therefore, we need to add h0 to get the absolute helix coordinate (h). Likewise, we need to subtract h1 to return to a relative coordinate system.

user/gee/regrid.c
void regrid( int dim         /* number of dimensions */, 
	     const int* nold /* old data size [dim] */, 
	     const int* nnew /* new data size [dim] */, 
	     sf_filter aa    /* modified filter */) 
/*< change data size >*/
{
    int i, h0, h1, h, ii[SF_MAX_DIM];

    for (i=0; i < dim; i++) {
	ii[i] = nold[i]/2-1;
    }
  
    h0 = sf_cart2line( dim, nold, ii); /* midpoint lag on nold */
    h1 = sf_cart2line( dim, nnew, ii); /*              on nnew */
    for (i=0; i < aa->nh; i++) { /* for all filter coefficients */
	h = aa->lag[i] + h0;
	sf_line2cart( dim, nold, h, ii);
	aa->lag[i] = sf_cart2line( dim, nnew, ii) - h1;
    }
}


next up previous [pdf]

Next: INVERSE FILTERS AND OTHER Up: The helical coordinate Previous: Filtering mammograms

2015-03-25