Subject: WCS

From: mcalabre at grus.csiro.au (Mark Calabretta,123,412,6608870)
Date: Thu, 25 Feb 1993 06:28:57 GMT

The appended latex document clarifies the relationship between the
elements of the CD matrix and the CDELTn and CROTAn parameters.  It was
developed from a document by Bob Hanisch and Phil Hodge (STScI) in
consultation with Eric Greisen (NRAO) and will appear as an appendix of
the WCS document currently available for review from
fits.cv.nrao.edu:/FITS/Documents/FITS_wcs.

Mark Calabretta,
Australia Telescope

\documentstyle[11pt]{article}
\topmargin 0in
\oddsidemargin 0in
\evensidemargin 0in
\headheight .125in
\footheight .125in
\textheight 8.9in
\textwidth 6.5in
\parindent 0in
\parskip 1ex
\widowpenalty 100000
\clubpenalty 100000
\nofiles

\begin{document}
\pagestyle{plain}
\thispagestyle{empty}

\begin{center} \large \bf
The CD Matrix in FITS and its Relationship to CDELTn and CROTAn \\
\normalsize \bf
\medskip
Mark Calabretta (ATNF) \\
Developed from a document by Bob Hanisch and Phil Hodge (ST ScI) \\
in consultation with Eric Greisen (NRAO) \\
\end{center}

The FITS {\tt CD} matrix describes a general linear transformation
from pixel coordinates $(i,j)$ to physical coordinates $(x,y)$, where
$(i,j)$ and $(x,y)$ are offsets from a reference point.  In this
appendix we discuss the relationship of the {\tt CD} matrix to the
{\tt CDELT}{\it n} and {\tt CROTA}{\it n} parameters.  In particular,
we will derive expressions for the elements of the {\tt CD} matrix in
terms of {\tt CDELT}{\it n} and {\tt CROTA}{\it n}.

Being a general linear transformation, the {\tt CD} matrix may
describe the geometrical operations of scaling, rotation and skewness.
Although scaling and rotation can be handled via {\tt CDELT}{\it n}
and {\tt CROTA}{\it n}, skewness cannot.  Moreover, the {\tt
CROTA}{\it n} have generally been interpreted in a very limited way by
FITS readers.  Specifically, only one {\tt CROTA}{\it n} parameter is
usually used, and it applies only to the two spherical coordinate axes
within the {\tt NAXIS}-dimensional FITS data cube.

A sequence of scaling and rotation operations sufficiently general for
astronomical purposes may be cast into matrix form as
\begin{equation}
   \left[ \begin{array}{c}  x \\
                            y \end{array} \right] =
   \left[ \begin{array}{cc} S_x & 0 \\
                            0 & S_y \end{array} \right] 
   \left[ \begin{array}{cc} D_x & 0 \\
                            0 & D_y \end{array} \right] 
   \left[ \begin{array}{rc}  \cos \theta & \sin \theta \\
                            -\sin \theta & \cos \theta \end{array} \right] 
   \left[ \begin{array}{cc} S_i & 0 \\
                            0 & S_j \end{array} \right] 
   \left[ \begin{array}{cc} D_i & 0 \\
                            0 & D_j \end{array} \right] 
   \left[ \begin{array}{c}  i \\
                            j \end{array} \right] \,.
\end{equation}
This describes a scaling operation performed on the pixel coordinates,
followed by a rotation, followed by another scaling.  Each of the
scaling matrices have been separated into two so that $D_i$, $D_j$,
$D_x$, and $D_y$ are positive quantities which describe the magnitude
of the scaling, and the $S_i$, $S_j$, $S_x$, and $S_y$ assume values
of $\pm 1$ and denote axis inversions.  Note that, in general, the
operations of scaling and rotation do not commute.

In the context of optical astronomy, equation (1) might
correspond to the following physical situation: one has a telescope
with a scale at the focal plane given by $D_x$ and $D_y$, forming an
image on a detector; the detector is rotated with respect to $(x,y)$
by the angle $\theta$ and it has pixels of size $D_i$ by $D_j$ (in
units such as mm/pixel).  $D_x$ and $D_y$ convert from mm on the
detector to the units of whatever is imaged onto the detector.  For an
image of the sky, the units of $D_x$ and $D_y$ would be, for example,
degrees/mm, while for a long-slit spectrograph the units could be
\AA/mm for one axis and degrees/mm for the other.  $S_i$ and $S_j$,
might assume values of $+1$ or $-1$ depending on whether the pixels
were read out of the detector in the forward or reverse direction, and
likewise the signs of $S_x$ and $S_y$ would depend on whether the
telescope optics had an even or odd number of reflections.  In the
radioastronomical context, only the first scaling operation is
required so that $D_x = D_y = 1$.

The {\tt CD} matrix formed by combining the above scaling and rotation 
matrices is
\begin{equation}
   \left[ \begin{array}{c}  x \\
                            y \end{array} \right] =
   \left[ \begin{array}{rc} S_x D_x S_i D_i \cos \theta 
                          & S_x D_x S_j D_j \sin \theta \\
                           -S_y D_y S_i D_i \sin \theta 
                          & S_y D_y S_j D_j \cos \theta \end{array} \right] 
   \left[ \begin{array}{c}  i \\
                            j \end{array} \right]
\end{equation}
whence
\begin{equation}
\begin{array}{lcr}
   {\tt CD}i\_i & = &   S_x D_x S_i D_i \cos \theta \,, \\
   {\tt CD}i\_j & = &   S_x D_x S_j D_j \sin \theta \,, \\
   {\tt CD}j\_i & = & - S_y D_y S_i D_i \sin \theta \,, \\
   {\tt CD}j\_j & = &   S_y D_y S_j D_j \cos \theta \,.
\end{array}
\end{equation}

In restricted cases, it is possible to associate the elements of the
{\tt CD} matrix with the {\tt CDELT}{\it n} and {\tt CROTA}{\it n}
parameters.  Consideration of the special case where the rotation
angle is zero allows us to identify the {\tt CDELT}{\it n} terms as 
\begin{equation}
\begin{array}{lcr}
   {\tt CDELT}i & = & S_x D_x S_i D_i \,, \\
   {\tt CDELT}j & = & S_y D_y S_j D_j \,.
\end{array}
\end{equation}
However, if the rotation angle is non-zero, the matrix elements of
equation (2) generally cannot be expressed as functions
of the {\tt CDELT}{\it n} and {\tt CROTA}{\it n} and so constraints
must be sought.  Possible constraints might be $D_i = D_j$ or $D_x =
D_y$.  Of these, the latter is appropriate in the radioastronomical
context and it is also the more plausible in optical astronomy.  (Note
that $D_x \neq D_y$ introduces a skew when the rotation angle is
not an integer multiple of $90^{\circ}$, a concept not supported in
the old notation.)

Now, since $S_x^2 = S_y^2 = +1$, we have 
\begin{equation}
\begin{array}{lcl}
   S_x S_j & = & S_x (S_y S_y) S_j   \\
           & = & (S_x S_y) \, (S_y S_j) \,.
\end{array}
\end{equation}
Likewise,
\begin{equation}
   S_y S_i = (S_x S_y) \, (S_x S_i) \,. \\
\end{equation}
Then, with $D_x = D_y$, and noting that $S_x S_y = \pm 1$ so that
$\cos \theta = \cos (S_x S_y \theta)$ and $S_x S_y \sin \theta = \sin
(S_x S_y \theta)$, equations (3) and (4) may
be combined to yield
\begin{equation}
\begin{array}{lcr}
   {\tt CD}i\_i & = &   {\tt CDELT}i\, \cos(S_x S_y \theta) \,, \\
   {\tt CD}i\_j & = &   {\tt CDELT}j\, \sin(S_x S_y \theta) \,, \\
   {\tt CD}j\_i & = & - {\tt CDELT}i\, \sin(S_x S_y \theta) \,, \\
   {\tt CD}j\_j & = &   {\tt CDELT}j\, \cos(S_x S_y \theta) \,.
\end{array}
\end{equation}

Likewise, equations (3) may be reduced to the form of
the Hanisch and Wells (1988) paper, by noting from
equations (4) that
\begin{equation}
\begin{array}{rcl}
   \left | {\tt CDELT}i \right | & = & D_x D_i  \\
                                 & = & D_y D_i \,, \\
   \left | {\tt CDELT}j \right | & = & D_y D_j  \\
                                 & = & D_x D_j \,, \\
        {\rm sign}({\tt CDELT}i) & = & S_x S_i \,, \\
        {\rm sign}({\tt CDELT}j) & = & S_y S_j \,,
\end{array}
\end{equation}
whence
\begin{equation}
\begin{array}{lcrcll}
   {\tt CD}i\_i & = &&        {\tt CDELT}i
                    & \cos(S_i S_j \theta) \,, \\
   {\tt CD}i\_j & = &   {\rm sign} ( {\tt CDELT}i )
                    & \left | {\tt CDELT}j \right |
                    & \sin(S_i S_j \theta) \,, \\
   {\tt CD}j\_i & = & - {\rm sign} ( {\tt CDELT}j )
                    & \left | {\tt CDELT}i \right |
                    & \sin(S_i S_j \theta) \,, \\
   {\tt CD}j\_j & = &&        {\tt CDELT}j
                    & \cos(S_i S_j \theta) \,.
\end{array}
\end{equation}
Note that the arguments of the trigonometric functions in
equations (7) and (9) may differ in sign.

At this point we have yet to associate the {\tt CROTA}{\it n} with
$\theta$.  The original FITS specification given by Wells, Greisen,
and Harten (1981) did not define how the {\tt
CROTA}{\it n} were to be interpreted but instead left each software
package to provide its own interpretation:
\begin{quote}
``The rotation parameters {\tt CROTA}{\it n} describe a coordinate
system which is rotated with respect to the normal (specified)
coordinate system.  Users of this option should provide extensive
explanatory comments.''
\end{quote}
Consequently, different interpretations are possible for different
packages.  It was stated above that, conventionally, only the {\tt
CROTA}{\it n} applying to the pair of spherical coordinate axes was
used.  The scheme adopted by AIPS and other packages has
\begin{equation}
\begin{array}{lcl}
   x & = & {\tt CDELT}i \,\cos ({\tt CROTA}j) \,i -
           {\tt CDELT}j \,\sin ({\tt CROTA}j) \,j \,, \\
   y & = & {\tt CDELT}i \,\sin ({\tt CROTA}j) \,i + 
           {\tt CDELT}j \,\cos ({\tt CROTA}j) \,j \,,
\end{array}
\end{equation}
where $i$ is the longitude-like axis and $j$ is the latitude-like
axis.  Comparing these with equations (7) the
association
\begin{equation}
   S_x S_y \theta = - {\tt CROTA}j
\end{equation}
may therefore be made, and the same conclusion can be drawn by
comparison with equations (9).

   The final form of the equations for the elements of the {\tt CD}
matrix in terms of the {\tt CDELT}{\it n} and {\tt CROTA}{\it n} is
\begin{equation}
\begin{array}{lcr}
   {\tt CD}i\_i & = &   {\tt CDELT}i \,\cos({\tt CROTA}j) \,, \\
   {\tt CD}i\_j & = & - {\tt CDELT}j \,\sin({\tt CROTA}j) \,, \\
   {\tt CD}j\_i & = &   {\tt CDELT}i \,\sin({\tt CROTA}j) \,, \\
   {\tt CD}j\_j & = &   {\tt CDELT}j \,\cos({\tt CROTA}j) \,.
\end{array}
\end{equation}

These equations may also be derived directly from
equations (9) by noting from equations (8)
and (11) that 
\begin{equation}
\begin{array}{lcl}
   S_i S_j \theta & = & S_i \,( S_x S_x )\, ( S_y S_y )\, S_j \,\theta \\
                  & = & ( S_i S_x )\, ( S_x S_y )\, ( S_y S_j )\, \theta \\
                  & = & - {\rm sign} ( {\tt CDELT}i ) \,
                          {\rm sign} ( {\tt CDELT}j ) \,{\tt CROTA}j \,.
\end{array}
\end{equation}
Although there was general agreement to support the conventions used
by AIPS in the NASA-sponsored conference of 1988 (Hanisch and Wells,
1988), the ambiguities in the original FITS papers mean
that there is no absolute guarantee that all packages have interpreted
the {\tt CROTA}{\it n} in accordance with equations (10).

\end{document}