[fitsbits] top or buttom

[Mark Calabretta]

On Tue 2000/04/25 11:12:58 MST, Stephen Walton wrote
in a message to: fitsbits at cv.nrao.edu

>J=1 to N.  Thus, C will be dimensioned (L,N).  This is standard
>mathematical matrix multiplication carried over into Fortran, where the
>first index of a matrix is the row index and the second is the column
>index.  Is there any equivalent convention in C/C++?

The first index in FORTRAN is whatever you want it to be, likewise in C,
it's just a storage notation.  For programming purposes you only need to
know the storage order, i.e. "which index varies fastest".

However, when dealing with matrices it is usual (but not necessary) to
consider the first index to be the row number of the matrix in both
FORTRAN and C.  That's why FORTRAN is said to have column-major storage
(first index varies fastest) while C has row-major storage (first index
varies slowest) .

Neither is better than the other when performing a matrix multiply since
in C(I,J) = A(I,K)*B(K,J), the summation index, K, runs over the rows of
one matrix and the columns of the other.

When you read a FITS image into ARRAY(NAXIS1,NAXIS2) in FORTRAN (or MATLAB)
it is not useful to think of ARRAY as being a matrix in the mathematical
sense, at least not for visualization purposes - the rows of the image
become the columns of such a matrix (your MATLAB problem).  You just need
to know that NAXIS1 varies fastest for FITS and for ARRAY so they're a good
match for FITS I/O.

>When we display FITS images with the first index varying horizontally
>and the second vertically, the first index is, in effect, the column
>index if we think (erroneously, I suppose) of the FITS data as a matrix
>in the mathematical sense.  MATLAB thinks of an array representing image
>data as a mathematical matrix and its Image Toolbox displays the array
>with its first index increasing vertically downward, consistent with the
>way a mathematics text displays a matrix.

While MATLAB uses matrix-oriented display conventions, we may think of the
FITS convention as being Cartesian-oriented; the simple FITS linear
coordinate mapping when reduced to simplest terms, i.e.

   (i,j,k, ...) -> (x,y,z, ...)

and compared to the usual Cartesian convention suggests that i, like x,
runs horizontally from left to right.  For a right-handed system the
remaining axes are then determined - j runs vertically from bottom to
top, and so on.

>How do you all normally treat the CD matrix?  This is supposed to be a
>true mathematical matrix with CD1_2 being the element in the first row,
>second column.  Does this get stored in your programs as element [1,2]
>or as element [2,1]?  It seems to me that either choice gives a
>problem.  Suppose I write a standard matrix-vector multiplication
>routine which takes an array A(M,N) and a vector V(N) and computes as
>output D(I) = A(I,K)*V(K) for I=1 to M.  I can compute the (x,y) WCS
>coordinates from pixel coordinates (I,J) if I set A(I,J) = CD_IJ, V(1) =
>I, V(2) = J.  However, now the corresponding (x,y) are for element
>ARRAY[J,I] if you dimension ARRAY[NAXIS2,NAXIS1], aren't they?

This is an implementation issue - the CD matrix and image array are
stored in the most efficient and convenient way for the programmer
implementing the standard.

For the record, the main difficulty in handling the CD array is that of
not knowing its size in advance.  In the FORTRAN implementation of WCSLIB,
for example, since FORTRAN doesn't have structs or pointers, I decided to
create a pseudo-data structure with everything packed into one array.  The
portion containing the CD matrix looks like:

*      LIN(4+MAXIS)
*         The first element of an array of size MAXIS*MAXIS containing
*         the NAXIS*NAXIS elements of the PC (pixel coordinate)
*         transformation matrix stored in column-major order with padded
*         columns, that is PC(1,1), PC(2,1), ..., PC(MAXIS,1), PC(1,2),
*         PC(2,2) ..., where PC(1,2) = PC1_2, etc.

Other elements of the array contain control data, MAXIS, NAXIS, the
CRPIXn, and the inverse CD matrix.

It was not sufficient (and wouldn't have helped much) to use COMMON
alignment for this "data structure" since it needs to exist as a separate
entity which can be passed as a subroutine argument so that an indefinite
number of coordinate systems can be maintained.

The C implementation of WCSLIB uses structs containing pointers to arrays
and is rather cleaner.  However, not knowing the size of the CD matrix in
advance still forces us to handle it as a 1-D array according to the
following prescription (from an unreleased version of WCSLIB):

*      double *cd
*         Pointer to the first element of the CD (pixel coordinate)
*         transformation matrix.  The expected order is
*
*            lin.cd = {CD1_1, CD1_2, CD2_1, CD2_2};
*
*         This may be conveniently constructed from a two-dimensional array
*         via
*
*            double m[2][2] = {{CD1_1, CD1_2},
*                              {CD2_1, CD2_2}};
*
*         which is equivalent to,
*
*            double m[2][2];
*            m[0][0] = CD1_1;
*            m[0][1] = CD1_2;
*            m[1][0] = CD2_1;
*            m[1][1] = CD2_2;
*
*         for which the storage order is
*
*            CD1_1, CD1_2, CD2_1, CD2_2
*
*         so it would be legitimate to set lin.cd = *m.


As for the image array, the only reason for storing it as
ARRAY[NAXIS2][NAXIS1] in C would be for I/O efficiency.  The consequent
inconvenient reversal of indices is the price paid for the mismatch
between the FITS and C storage conventions.  However, FORTRAN is a good
match and they can't both be good!

You could store the image array in ARRAY[NAXIS1][NAXIS2] in C for
bookkeeping convenience but would then have to pay the price of one
transpose on read and another on write.

Mark Calabretta
ATNF