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,

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.regridConvert filter to different data sizeWe 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.randshould be larger than (or equal to)centerand smaller thanmin (nold, nnew), otherwise the filter might stick outside the grid (our piece of paper.)rand=nold/2will do (assuming the filter is small), although nothing should change if you replacenold/2with a random integer array betweencenterandnold - na.The linear coordinate of

randish0on the old helix andh1on the new helix. Recall that the helix lagsaa%lagare relative to the center. Therefore, we need to addh0to get the absolute helix coordinate (h). Likewise, we need to subtracth1to return to a relative coordinate system.

4/27/2004