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 createhelixmod . We begin here from its intricate underpinnings.

Fortran77 has a concept of a multidimensional array being equivalent to a one-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 one-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 doesn't fit very well, since it is intrinsically noncausal. If you really want noncausal filters, you will 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 one-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 one-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 one-dimensional index ia is

        ia = i1 +a1*(i2-1) +a1*a2*(i3-1) + ...

When the filter cube is stored in an array with the same dimensions as the data, data(n1,n2,n3,...), the formula for ia is

        ia = i1 +n1*(i2-1) +n1*n2*(i3-1) + ...

The fortran compiler knows how to convert from the multidimensional cartesian indices to the linear index. We will need to do that, as well as the converse. Module cartesian below 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. cartesianhelical-cartesian coordinate conversion

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'' which is 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 cartesian module fails for negative subscripts. Thus we need to be careful to avoid thinking of the filter's 1.0 (shown in Figure ) as the origin of the multidimensional coordinate system although the 1.0 is the origin in the one-dimensional coordinate system.

Even in one dimension (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, we also need to know 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 hypecube.

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. This would waste so much space, however, that our helix filter module helix instead stores the filter coefficients in one vector and their 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 createhelixmod 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  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 createhelixmod 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 their lags can be ignored too.) 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. createhelixmodconstructing helix filter in N-D Near the end of the code you see the calculation of a parameter lag0d. This is the count of the number of zeros that a data-sized fortran array would store in a filter cube before 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 (7) shows us in two dimensions or like Figure  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 box: boxConvert helix filter to (n1,n2,...) The box module is normally used to display or manipulate a filter that was estimated in helical form (usually estimated by the least-squares method).

The inverse process to box is to convert a fortran hypercube to a helix filter. For this we have module unbox. It abandons all zero-valued coefficients such as those that should be zero before the box's 1.0. It abandons the ``1.0'' as well, because it is implicitly present in the helix convolution module helicon . unboxConvert hypercube filter to helix An example of using unbox would be copying some numbers such as the factored laplacian in equation (11) into a cube and then converting it to a helix.

A reasonable arrangement for a small 3-D filter is na=(5,3,2) and center=(3,2,1). Using these arguments, I used createhelixmod to create a filter. I set all the helix filter coefficients to 2. Then I used module box to put it in a convenient form for display. After this conversion, the coefficient aa(3,2,1) is 1, not 2. 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's 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.

regridConvert filter to different data size