[fitswcs] Polynomial Correction Function & Alternative-WCS
Don Wells
dwells at NRAO.EDU
Wed May 6 14:08:04 EDT 1998
Polynomial Correction Function, Alternative-WCS Descriptions
Don Wells
DRAFT-for-comments 1998-05-06
-=- Introduction -=-
This memo proposes an n-D polynomial correction function specified by
keywords in the FITS header to replace the IMAGE extension numerical
map of the correction which was proposed in Appendix A of Greisen and
Calbretta (http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps).
This new concept supports coding correction values for IMAGE
extensions as well as for primary matricies of FITS files. The CPn
keywords which define this function are capable of replacing the
functionality of CRVALi, CDi_j, CDELTi and PCiiijjj, if the community
wishes to make this change. This memo also proposes a keyword
convention for alternative-WCS descriptions.
-=- Polynomial Correction Function -=-
The correction function specifies a perturbed position in each axis as
a function of the positions in all axes. The polynomial coefficients
will be specified by keywords of the form:
CPn = <E-format-value> / coeff of polynomial term n
In the NAXIS=1 (spectroscopic) case, the ordinal 'n' will be computed as
n = 1000 + i
where 'i' is the power 0,1,2,.. of the pixel coordinate. In the
NAXIS=2 (image) case, the ordinal 'n' will be computed as
n = axis*1000 + ((i + j + 1) * (i + j)) / 2 + j
where 'axis' is the ordinal 1..NAXIS of the polynomial, 'i' is the
power 0,1,2,.. of the first pixel axis and 'j' is the power of the
second pixel axis. For the NAXIS=2 case, the tableau of terms for Y'
(perturbed-Y) is
CP2000*1 CP2001*X CP2003*X^2 CP2006*X^3
CP2002*Y CP2004*X*Y CP2007*X^2*Y CP2011*X^3*Y
CP2005*Y^2 CP2008*X*Y^2 CP2012*X^2*Y^2 CP2017*X^3*Y^2
CP2009*Y^3 CP2013*X*Y^3 CP2018*X^2*Y^3 CP2024*X^3*Y^3
where X = (i - CRPIX1) and Y = (j - CRPIX2). Note the diagonal
pattern of the increasing ordinal. The capability for 999 terms in
the polynomial will support terms up to q=(i+j)=21 (X^21, X^20Y,
X^19Y^2, .., XY^20, Y^21) for the NAXIS=2 case. If five characters are
available for 'n' (see the alternative-WCS suffix proposal below) the
keyword names will support up to NAXIS=99. The inversion from n back
to i, j and axis for NAXIS=2 is (expressed in AWK):
axis = int(n / 1000);
np = n - axis * 1000;
ip = int((sqrt(8 * np + 1) - 1) / 2);
j = int(np - (ip *(ip + 1) / 2));
i = ip - j;
The value of n which includes all terms with (i+j)<=q is obtained by
setting i=0, j=q in the formula for n given above:
n_q = axis*1000 + ((q + 1) * q) / 2 + q)
For NAXIS=2 and q=3 (cubic), the following keywords are used:
CP1000 = <E-format-value> / coeff for 1 for axis 1
CP1001 = <E-format-value> / coeff for X for axis 1
CP1002 = <E-format-value> / coeff for Y for axis 1
CP1003 = <E-format-value> / coeff for X^2 for axis 1
CP1004 = <E-format-value> / coeff for X*Y for axis 1
CP1005 = <E-format-value> / coeff for Y^2 for axis 1
CP1006 = <E-format-value> / coeff for X^3 for axis 1
CP1007 = <E-format-value> / coeff for X^2*Y for axis 1
CP1008 = <E-format-value> / coeff for X*Y^2 for axis 1
CP1009 = <E-format-value> / coeff for Y^3 for axis 1
CP2000 = <E-format-value> / coeff for 1 for axis 2
CP2001 = <E-format-value> / coeff for X for axis 2
CP2002 = <E-format-value> / coeff for Y for axis 2
CP2003 = <E-format-value> / coeff for X^2 for axis 2
CP2004 = <E-format-value> / coeff for X*Y for axis 2
CP2005 = <E-format-value> / coeff for Y^2 for axis 2
CP2006 = <E-format-value> / coeff for X^3 for axis 2
CP2007 = <E-format-value> / coeff for X^2*Y for axis 2
CP2008 = <E-format-value> / coeff for X*Y^2 for axis 2
CP2009 = <E-format-value> / coeff for Y^3 for axis 2
The author has not yet derived the general formula for n and its
inversion formulae for NAXIS>2. Somebody should do this, and should
determine the q=(i+j) order limit for the 999-term capability as a
function of NAXIS=1..7 (the largest NAXIS known to be potentially
useful).
Users will be free to use any analytical form (splines, Chebyshevs,
Zernikes, etc) for fitting or computing pixel corrections to data; a
sufficient number of CPn keywords can be supplied to represent any
such function with negligible error. Also, fitting/computing
implementations may use normalized coordinates, and may be referred to
origins different from the FITS reference pixel specified by the
CRPIXi. This is because the conventional (non-orthogonal) polynomial
correction function specified in this proposal is simply the universal
interchange convention for correction information, and users can
transform to and/or from local conventions which may be based on other
functional forms.
If a submatrix or a decimation is extracted from a matrix which has an
attached correction function, the function must be transformed to suit
the new matrix. The required transformations for the most common
operations will need to be derived.
The question of conventions for coding corrections to matricies in
fields of BINTABLE extensions is TBD.
It appears that this proposed function could replace the functionality
of CRVALi, CDi_j, CDELTi and PCiiijjj. The CDi_j (partial derivatives
at the reference point) notation is equivalent to the vector-matrix
product of the CDELTi and PCiiijjj. For the NAXIS=2 case, the CRVALi
and CDi_j can be expressed as:
CRVAL1 --> CP1000
CD1_1 --> CP1001
CD1_2 --> CP1002
CRVAL2 --> CP2000
CD2_1 --> CP2001
CD2_2 --> CP2002
This concept of replacing all of the prior notations with a single
generalized functional form is certainly attractive. The notation is
equivalent to the 'plate-constants' of classical optical astrometry.
The FITS community will need to decide whether they wish to make this
change or whether they would prefer to retain CRVALi and CDi_j and
apply the CPn correction function to transform the raw pixel
coordinates to perturbed pixel coordinates. An argument for doing the
latter is that when a submatrix, decimation and/or reflection of a
matrix has been produced, the original WCS description could be
retained if a 4-term linear transformation of the pixel coordinates
were specified with the CPn (using CP1000, CP1001, CP2000 and CP2002
for NAXIS=2 case). The pros and cons of these possibilities are
matters on which reasonable people can differ; the author currently
takes no position on the question.
An API for a polynomial correction function in C could be:
int wcsPolyCorr( /* return != 0 indicates error */
int direction, /* 0: x[]-->xp[], !0: xp[]-->x[] */
int naxes, /* NAXIS value */
int ncoeff, /* number of cp[] coefficients */
int cpa[] /* axis number for coefficients */
int cpi[][], /* powers of pixel coords */
double cp[], /* coefficient values */
double x[], /* original coordinates */
double xp[]) /* perturbed coordinates */
The formal parameters x[] and xp[] are either input or output
arguments depending on the value of the 'direction' argument. The
cpi[] and cp[] arrays are filled from the FITS keywords in the header;
the inversion formulae should be used to set cpa[ncoeff] with the axis
number and cpi[ncoeff][naxis] with the powers of the pixel coordinates
for the cp[] values. This initialization will avoid recomputation of
the index inversion formulae for every pixel in an n-D data matrix.
When 'direction' is non-zero the values in xp[] will be copied to x[]
and the forward calculation performed to get the perturbed position in
all axes. The coordinate change (the perturbation) will then be
subtracted from x[] and the calculation repeated until convergence.
The mapping defined by the cp[] coefficients should be 1-to-1 (e.g.,
monotonic increasing). Implementors should note that in the case
where only CPi000 through CPi002 are present, the mapping is linear
and is invertible by analytic means.
+------------------------------+
| Alternative WCS descriptions |
+------------------------------+
A matrix of pixels can usually be mapped to more than one World
Coordinate System. For example, in imaging systems pixels can be
mapped to linear coordinate systems associated with the detector
and/or telescope as well as to spherical coordinates in the sky; both
types of WCS mappings are valid, both are useful and we should be able
to carry both in the FITS headers. Another example is radio
spectroscopy, where the pixels containng emission from a particular
spectral line can be mapped either to observed frequency or to Doppler
velocity (using the rest frequency of the spectral line); again, both
WCS notations are useful and we should be able to carry both in the
FITS headers.
This memo proposes that an optional suffix character be used to
distinguish keywords for such alternative-WCS descriptions. The
character can take values of blank or [A-Z], which will support up to
26 alternative-WCS descriptions of a matrix.
For example, suppose that we wish to specify an alternative quadratic
version of the cubic correction polynomial whose terms were listed
above. We would append a suffix ('Q' in this example) to the CRPIXi
and to the 12 CPn keywords:
CRPIX1Q = 327.76 / X-axis
CRPIX2Q = 225.43 / Y-axis
CP1001Q = <E-format-value> / coeff for 1 for axis 1
CP1002Q = <E-format-value> / coeff for X for axis 1
CP1003Q = <E-format-value> / coeff for Y for axis 1
CP1004Q = <E-format-value> / coeff for X^2 for axis 1
CP1005Q = <E-format-value> / coeff for X*Y for axis 1
CP1006Q = <E-format-value> / coeff for Y^2 for axis 1
CP2001Q = <E-format-value> / coeff for 1 for axis 2
CP2002Q = <E-format-value> / coeff for X for axis 2
CP2003Q = <E-format-value> / coeff for Y for axis 2
CP2004Q = <E-format-value> / coeff for X^2 for axis 2
CP2005Q = <E-format-value> / coeff for X*Y for axis 2
CP2006Q = <E-format-value> / coeff for Y^2 for axis 2
We would append 'Q' to all WCS keywords, such as CTYPEiQ.
It may be appropriate to provide a label string for alternative-WCS
mappings. The author suggests WCSNAME and WCSNAMEv (WCSNAMEQ in the
case above):
WCSNAME = 'Cubic Polynomial Case' /
WCSNAMEQ= 'Quadratic Polynomial Case' /
-=- wcs_pcf.awk: the program which computed the terms above -=-
#!/opt/local/gnu/bin/awk -f
# wcs_pcf.awk --- proof-of-concept test for Polynomical Correction Function
# D.Wells, NRAO-CV, 1998-05-06.
#
BEGIN {
# compute a matrix of polynomial terms up to i-limit and j-limit:
ij_max = 3; axis = 2;
printf("\n");
for (j = 0; j <= ij_max; j++) {
for (i = 0; i <= ij_max; i++) {
n = axis*1000 + ((i + j + 1) * (i + j)) / 2 + j;
printf(" CP%d*%-8s", n, term(i, j));
}
printf("\n\n");
}
# compute a list of keywords up to q=(i+j) exponent limit for NAXIS=2:
q = 3;
for (axis = 1; axis <= 2; axis++) {
for (n = nq(axis,0); n <= nq(axis,q); n++) {
axis = int(n / 1000);
np = n - axis * 1000;
ip = int((sqrt(8 * np + 1) - 1) / 2);
j = int(np - (ip *(ip + 1) / 2));
i = ip - j;
keyword = sprintf("CP%d", n);
printf(" %-8s= <E-format-value> / coeff for %-5s for axis %d\n",
keyword, term(i, j), axis);
}
}
exit;
}
function nq(axis, q) {
return(axis*1000 + ((q + 1) * q) / 2 + q);
}
function term(i, j, termq) {
termq = term_ij("X", i) "*" term_ij("Y", j);
if (termq ~ /1\*.+/) termq = substr(termq, 3);
if (termq ~ /.+\*1/) termq = substr(termq, 1, length(termq)-2);
return(termq);
}
function term_ij(a, i, term_q) {
if (i == 0) {
term_q = "1";
} else {
if (i == 1) {
term_q = a;
} else {
term_q = sprintf("%s^%d", a, i);
}
}
return(term_q);
}
More information about the fitswcs
mailing list