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. The formula for computing ia is

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

for a cube stored in its normal ``tightly packed'' form, or if the cube were stored in an array with the same dimensions as the data, 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.  

#$
#$=head1 NAME
#$
#$cartesian - convert to and from cartesian
#$
#$=head1 SYNOPSIS
#$
#$C<call line2cart(nn,i,ii)>
#$
#$C<call cart2line(nn,ii,i)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item nn - C<int(:)>
#$
#$      cartesian axes
#$
#$=item ii - C<int(:)>  
#$
#$      cartesian coordinates
#$
#$=item i  - int        
#$
#$      linear index
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$index transform (vector to matrix) and its inverse
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$

module cartesian {       # index transform (vector to matrix) and its inverse
contains
  subroutine line2cart( nn, i, ii) { 
    integer, dimension( :), intent( in) :: nn    # cartesian axes (n1,n2,n3,...)
    integer, dimension( :), intent(out) :: ii    # cartesn coords (i1,i2,i3,...)
    integer               , intent( in) :: i     # equivalent 1-D linear index
    integer                             :: axis, n123
    n123 = 1
    do axis = 1, size(  nn) { 
	ii( axis) = mod( ( i-1)/n123, nn( axis)) + 1
	n123 = n123 * nn( axis)
        }
    }
  subroutine cart2line( nn, ii, i) { 
    integer, dimension( :), intent( in) :: nn, ii
    integer                             :: i, axis, n123
    n123 = 1;    i = 1
    do axis = 1, size(  nn) { 
	i = i + ( ii( axis)-1)*n123
	n123 = n123 * nn( axis)
	}
    }
}

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 11) 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 11 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.  

#$=head1 NAME
#$ 
#$createhelix - create a helix filter
#$  
#$=head1 SYNOPSIS
#$
#$C<aa=createhelix(nd,center,gap,na)>
#$
#$=head1 INPUT PARAMETERS
#$
#$=over 4
#$
#$=item nd      - C<integer(:)>
#$
#$      size of data
#$
#$=item center  - C<integer(:)>
#$
#$      location of the 1 of the filter within na box
#$
#$=item gap     - C<integer(:)>
#$
#$      distance along each axis before filter coef
#$
#$=item na      - C<integer(:)>
#$
#$      size of box arround filter
#$
#$=back
#$
#$=head1 RETURN VALUE
#$
#$=over 4
#$
#$=item aa       - C<type(filter)>
#$
#$      Helix filter
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$Create helix filter
#$
#$=head1 SEE ALSO
#$
#$L<createnhelix>,L<helix>,L<createmshelix>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
module createhelixmod {                      # Create helix filter lags and mis
use helix
use cartesian
contains
  function createhelix( nd, center, gap, na)  result( aa) {
    type( filter)                     :: aa       # needed by helicon.
    integer, dimension(:), intent(in) :: nd, na   # data and filter axes
    integer, dimension(:), intent(in) :: center   # normally (na1/2,na2/2,...,1)
    integer, dimension(:), intent(in) :: gap      # normally ( 0,   0,  0,...,0)
    integer, dimension( size( nd))    :: ii       # cartesian indexes
    integer                           :: na123, ia, ndim, nh, lag0a,lag0d
    integer, dimension(:), allocatable:: lag
           nh= 0;   na123 = product( na);   ndim = size( nd)
    allocate(  lag( na123 ) )		  # filter cube size
    call cart2line ( na, center, lag0a)   # lag0a = index pointing to the "1.0"
    do ia = 1+lag0a, na123 {              #    ia  is fortran linear index.
       call line2cart( na, ia, ii)        # ii(ia) is fortran array  indices.
       if( any( ii <= gap))    next       # ignore some locations
       nh = nh + 1                        # got another live one.
       call cart2line( nd, ii, lag(nh))	  # get its fortran linear index
       }
    call cart2line(  nd, center, lag0d)   # lag0d is center shift for nd_cube
    call allocatehelix( aa, nh)		  # nh becomes size of filter on helix.
    aa%lag = lag(1:nh) - lag0d;		  # lag = fortran_linear_index - center
    aa%flt = 0.0;             deallocate( lag)
    }      
}

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 11 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:  

#$
#$=head1 NAME
#$
#$box  - filter to hypercube
#$
#$=head1 SYNOPSIS
#$
#$C<call boxn (nd,center,na,aa,cube)>
#$
#$=head1 INPUT PARAMETERS
#$
#$=over 4
#$
#$=item nd  - C<int(:)>
#$
#$      Input data dimension
#$
#$=item center - C<int(:)>  
#$
#$      Location of the 1 coeficient in box
#$
#$=item na  -  C<int(:)>    
#$
#$      Dimensions of the filter
#$
#$=item aa  -
#$
#$      C<type(filter)> The filter
#$
#$=back
#$
#$=head1 OUTPUT PARAMETERS
#$
#$=over 4
#$
#$=item cube  -  real(:)      
#$
#$      Output filter values
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$ Convert helix filter to hypercube: cube(na(1),na(2),...)
#$
#$=head1 SEE ALSO
#$
#$B<helix>, B<cartesian>, B<helixcart>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$

module box {         # Convert helix filter to hypercube: cube(na(1),na(2),...)
use helix
use cartesian
contains
  subroutine boxn( nd, center, na, aa, cube) {
    integer, dimension (:), intent( in)  :: nd, center, na    # (ndim)
    type( filter),          intent( in)  :: aa
    real,    dimension( :), intent( out) :: cube 
    integer, dimension( size( nd))       :: ii
    integer                              :: j, lag0a, lag0d, id, ia
    cube = 0.;                          # cube=0
    call cart2line( na, center, lag0a)  # locate the 1.0 in the na_cube.
    cube( lag0a) = 1.                   # place it.
    call cart2line( nd, center, lag0d)  # locate the 1.0 in the nd_cube.
    do j = 1, size( aa%lag) {		# inspect the entire helix
       id = aa%lag( j) + lag0d		# index = helix_lag + center_d
       call line2cart( nd, id, ii)	# ii(id) = cartesian indices
       call cart2line( na,     ii, ia)	# ia(ii) = linear index in aa
       cube( ia) = aa%flt( j)		# copy the filter coefficient
       }
    }
}

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 .  

module unbox {				# helixfilter aa = cube(a1,a2,...)
use helix
use cartesian
contains
  function unboxn( nd, center, na, cube)  result( aa) {
    type( filter)                       :: aa
    integer, dimension( :), intent( in) :: nd, center, na      # (ndim)
    real,    dimension( :), intent( in) :: cube		       # cube(a1,a2,...)
    logical, dimension( size( cube))    :: keep		       # keep(a1*a2*...)
    integer, dimension( size( nd))      :: ii		       # (ndim)
    integer                             :: ic, lag0a, lag0d, i, h
    call cart2line(  na, center, lag0a)
    call cart2line(  nd, center, lag0d)
    keep = ( abs( cube) > epsilon( cube))     # epsilon is a Fortran intrinsic
    keep( lag0a) = .false.                    # throw away the 1.0.
    call allocatehelix( aa, count( keep));  h = 0
    do ic = 1, size( cube) {			       # sweep cube
           if( keep( ic) ) {                h = h + 1  # only the keepers
                call line2cart(  na, ic, ii)          # ii(ic)= indices on na
                call cart2line(  nd,     ii, i)       # i     = index   on nd
                aa%lag( h) = i - lag0d                 # lag = index - center
                aa%flt( h) = cube( ic)		       # copy coefs.
                }
           }
    }
}

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.

 

#$=head1 NAME
#$
#$regrid - convert a helix filter from one data space to another
#$
#$=head1 SYNOPSIS
#$
#$C<call regridn(nold,nnew,aa)>
#$  
#$=head1 PARAMETERS
#$  
#$=over 4
#$  
#$=item nold - C<integer(:)>
#$  
#$      Old data dimensions
#$
#$=item nnew - C<integer(:)>
#$  
#$      New data dimensions
#$
#$=item aa   -  type(filter)
#$  
#$      Helix filter 
#$  
#$=back 
#$
#$=head1 DESCRIPTION
#$  
#$Converts a helix filter from one ata size to another
#$
#$=head1 SEE ALSO
#$  
#$L<cartesian>,L<helix>
#$  
#$=head1 LIBRARY
#$  
#$B<geef90>
#$  
#$=cut

module regrid {         # convert a helix filter from one size data to another
use helix
use cartesian
contains
  subroutine regridn( nold, nnew, aa) {
    integer, dimension (:), intent (in) :: nold, nnew  # old and new helix grid
    type( filter)                       :: aa
    integer, dimension( size( nold))    :: ii
    integer                             :: i, h0, h1, h
    call cart2line( nold, nold/2, h0)	# lag of any near middle point on nold
    call cart2line( nnew, nold/2, h1)	# lag                          on nnew
    do i = 1, size( aa%lag) {		# forall given filter coefficients
       h = aa%lag( i) + h0		# what is this?
       call line2cart( nold, h, ii)	# 
       call cart2line( nnew,    ii, h)	# 
       aa%lag( i) = h - h1		# what is this
    }
  }
}