[NOAO logo] [IRAF logo]

Functional Corrections and Axes Coupling in FITS WCS

Francisco Valdes (fvaldes@noao.edu)
Lindsey Davis (ldavis@noao.edu)
Doug Tody (dtody@noao.edu)
Draft: June 29, 2001


This document proposes that the model for transforming pixel coordinates to world coordinates in the FITS WCS consist of a linear transformation (the CD transformation), followed by a distortion/coupling function (DCF), followed by a world transformation (the various celestial and spectral transformations). This differs from the current model by explicitly including a common DCF applicable to all WCS transformations rather than incorporating it as specific cases in various world transformations. The proposal includes a syntax for describing this sequence in the CTYPE keywords.

This document also includes a sketch of some likely DCFs that can be defined for general use within the FITS WCS as well as a specific proposal for a polynomial DCF. The polynomial DCF is illustrated by several examples illustrating both its use for distortions and for coupling spatial and spectral axes. We emphasize that the ability to couple arbitrary axes through a DCF is a new aspect not provided by the current WCS proposals.

1. Introduction

The FITS world coordinate system (WCS) provides for a variety of coordinate system types. Under the current FITS WCS proposals (referred to as WCS-I, WCS-II, and WCS-III) each coordinate system type is entirely self-contained. The coordinate types provide two types of transformations, a specific geometric or mathematical transformation for an ideal instrument and a distortion transformation to apply corrections from a real instrument to the ideal instrument.

The problem with combining two types of transformations in a single self-contained definition is each and every coordinate system type needs to have its own distortion correction. Often the distortion correction will be the same for many coordinate system types leading to a great deal of repetitiveness in the definitions. For example, the current proposals provide a distortion correction for the tangent plane projection (TAN) but if one wants to uses the same distortion model with a cylindrical equal area projection (CEA) it would need to be included in the definition of that coordinate system type.

Another limitation of the current proposal is dealing with coupled axes. Currently to couple axes requires a set of coordinate type definitions for the set of axes to be coupled. For the ideal celestial transformations this makes sense and all the coordinate types are defined in pairs. But consider three-dimensional spectral data consisting of two celestial axes and a spectral axis. If distortions make the spectral coordinates dependent on spatial position then one cannot use the separate coordinate types for the spatial and spectral dimensions. Rather a set of three matched coordinate types would be needed, though in many aspects they would be the same as those for the separate spatial and spectral coordinate types.

Addressing the issue of coupled axes, especially dissimilar axes such as spatial and spectral axes, is a major argument for this proposal. In the current set of WCS proposals the issue of coupling spectral and spatial axes, as noted explicitly in WCS-III, is not addressed. This proposal has the generality to cover all types of coupling but was specifically developed to address the case of spectral data. Examples of applying this proposal to such data are included here as well as in Spectral WCS Conventions.

The proposal presented here separates the ideal two-dimensional celestial and one-dimensional spectral coordinate types from the distortion and coupling function (DCF). This allows using the existing ideal coordinate system types in a variety of situations without change and defines a general way to deal with distortions and axes coupling.

The sequence of transformations is shown in the following figure.

Figure 1: Sequence of Transformations

linear              distortion/coupling    ideal world         
transformation  ->  transformation      -> transformation      
(CRPIX & CD)        (DV/DS)                (projection & CRVAL)
Image pixel coordinates enter on the left and are transformed to intermediate coordinates. The intermediate coordinates are then optionally transformed to ideal intermediate coordinates. Finally the world coordinates are produced on the right.

This sequence is the same as described in the current WCS proposal. What this proposal does is separate the last two transformations rather than combine them under one coordinate type definition. The way axes are coupled is by assigning the same multivariate distortion correction to the intermediate coordinates of all the coupled axes while allowing the final ideal intermediate coordinates to be transformed to the final world coordinates through the independent or ideally coupled (e.g. RA/DEC) transformations.

Before going into the details of this proposal we need to say a few words about the pixel regularization image defined in the current WCS-I proposal. The purpose of this is also to correct data from an instrument with distortions to that of an ideal instrument. Clearly it can address the same concerns though it is located in a different place in the chain of transformations. However, it has several drawbacks. A separate image extension must be used and interpolated. Also many, if not most, types of distortions are physically describable by a function with a small set of parameters. Therefore, using a separate image to apply such a function is not an elegant, compact, and self-contained approach.

In section 3 there is a proposal for a type of DCF that provides a sampled correction. The pixel regularization image can be redefined to be a DCF of this type with a string parameter pointing to the image extension. As noted in the previous paragraph, this is not identical in that the correction occurs at a different point in the transformation chain, but otherwise the function and effect are equivalent. This is not to say we encourage use of a PRI. But by using the DCF mechanism it is possible to include it and remove it from the explicit description in WCS-I and the DCF might be defined in such a way that the sample points are contained in the same extension as the image.

2. General Form of the CTYPE Keywords

The CTYPE keywords contains all the information about the WCS transformation for each of the world coordinate axes. Encapsulating this in a single keyword, as is done in the current proposals, is reasonable. This proposal also uses the CTYPE keyword to completely identify the transformation from pixel to world coordinates for a world coordinate axes. But instead of treating the keyword value as the name of a single transformation it may identify two transformations.

In this proposal, any of the general CTYPE labels for the ideal transformations may be modified by appending a string, beginning with a hyphen, which identifies a type of functional correction and an instance to allow for multiple cases. So the general suffix format is -[type][instance]. For example, the string "-PLYa" defines instance or group a of a polynomial function "PLY". Note that a distortion/coupling function need not be present in which case the ideal coordinate transformation is applied with no intermediate transformation.

One interpretation of this format is that the transformations, excluding the standard linear CD transformation, may be read from right to left. So "DEC--TAN-PLY" says to transform with the PLY function, then transform with the TAN projection to end up with declination.

The dimensionality of the multivariate function and the axes coupled are determined by the number of axes which have the same function extension. For example, if three axes are coupled by a polynomial transformation, the three axes would all have the modifier "-PLY1". If there are additional axes which also require a polynomial correction they would use "-PLY2".

Figure 2 illustrates the form of the CTYPE keyword. This shows the keywords with a multiple WCS letter for generality though for the primary WCS there would be no letter. The first two world coordinate axes correspond to celestial RA and DEC and the third corresponds to wavelength. The celestial axes are projected from the ideal intermediate coordinates using the tangent plane model given in WCS-II. The wavelength intermediate coordinates are expressed in velocity coordinates using the wavelength to velocity methods of WCS-III. The fact that all three axes have "-PLY" indicates that the intermediate coordinates are transformed by a polynomial function of the three coordinates before applying the ideal projection or transformation.

Figure 2: CTYPE keyword example.


In figure 2 there is only one instance of the polynomial function which couples three axes. This is expected to be the most common situation and the optional instance character need not be used. It may not be clear why one needs multiple instances. Consider the above case of three axes where the velocity axis is not coupled to the celestial coordinates but does require an intermediate polynomial transformation to deal with a non-linear wavelength sampling. In the case the keywords might be

Figure 3: CTYPE keyword example with multiple instances.


In this case the first function, PLYa, would have two independent variables from the spatial axes (see section 4 to see that there are actually three variables where the third variable is the cartisian radius of the other variables) and the second function would be a function of only the wavelength axis. Note that if there are four axes one might have two pairs of coupled axes represented identified by PLYa and PLYb.

To illustrate the concept of different DCF types consider the example figure 3 where the spectral DCF is a sampled interpolation.

Figure 4: CTYPE keyword example with multiple DCF types.

CTYPE1s = 'RA---TAN-PLY' / Polynomial DCF
CTYPE2s = 'DEC--TAN-PLY' / Polynomial DCF
CTYPE3s = 'VELO-WAV-SMP' / Sampled DCF   

3. Distortion Coupling Transformation Types

This section identifies three possible distortion coupling transformations. An important point of the proposal presented here is that it defines a distortion coupling transformation stage in the WCS model and not just a specific transformation. This is a general mechanism by which a variety of DCFs may be defined by agreement and leaves a path open to defining new DCFs in the future.

The first possible DCF is a multivariate polynomial transformation. Polynomial functions have many advantages such as being separable in each coordinate and fully general for situations with low frequencies; that is smooth functions without high frequency wiggles. This type of DCF can replace the polynomial function of the WCS-II proposal for TAN projections as well as the non-linear part of the WCS-III proposal for 1D spectral coordinates. A specific proposal is given in section 4 .

The second possible DCF or class of DCF is a sampled transformation. This type of transformation interpolates between points where the function values are given explicitly. Examples of this type of transformation are splines, lookup tables, and the pixel regularization image. Note that this type of DCF can replace the PRI and can incorporate the practice of representing spectral coordinates by columns in spectra stored in tables with the WCS model.

A third possible DCF is one based on trigonometric functions. We suggest this as a possiblity because it has a physical instrumental connection. Grating equations describing the mapping between spatial incidence on a grating and the spectral coordinate in terms of trigonometric functions would naturally fit this type of DCF.

There is the potential for confusion between the multiple WCS identification character and the instance character defined here. The instance character only appears in the CTYPE value string. It is used to identify the set of axes which are coupled by a particular DCF. The same DCF type may be used to couple independent sets of axes thus requiring the instance character. However, the instance character is not used in the keywords for the parameters of the DCF. The parameter keywords will use the multiple WCS character.

4. PLYs: Polynomial Corrections including Radial Terms

Polynomial corrections of any number of variables are identified by the string -PLYi, where i is an optional instance character, at the end of the CTYPE value string. The coefficients for this function are given by keywords DVl_ns, where l is the axis number, n is a coefficient number, and s is a letter for multiple WCS (which is distinct from the character used for multiple instances of the PLY function). Currently no string valued parameters are defined for the PLY function but the keywords for string valued parameters are DSl_ns. The axis number and the coefficient number are written with no leading zeros. The number of digits for n is the maximum allowed after the axis number and the optional multiple WCS letter are counted. Thus, for a single digit axis number with no multiple WCS letter the coefficient number may run from 0 to 9999. If a coefficient keyword is absent the default value is zero.

Let x_l and x'_l represent the intermediate coordinate for WCS axis l before and after the correction is applied. Since the set of N axes coupled by having the same PLYs identification may be any subset of the WCS axes, we let m index just the axes in the polynomial in increasing order of the WCS axes. For example, if the WCS axes covered by the polynomial are 1, 3 and 6 then m=1 for l=1, m=2 for l=3, and m=3 for l=6.

The polynomial correction is defined by a set of coefficients DVl_ns. The association between the coefficient index n and the elements of the transformation is given following the function description. The letter s is optional and is used when multiple WCS are defined.

The first step is to normalize the coordinates.

    (1) n_m = x_l / (1 + DVl_ns)
This is defined in this form so that when the DV coefficients are zero or absent there is no normalization. The normalization is optional but it is useful to avoid taking large powers of very large or small numbers. It is also useful in reducing dissimilar axes to similar magnitudes.

The next step is to compute an additional variable r given by

    (2) r = sqrt (sum_m DVl_ns n_m^2)
Whether or not this extra variable is used and which axes contribute to it is determined by the DV coefficients. Normally the DV coefficients in (2) will be either zero or one. With the addition of the r variable the dimensionality of the polynomial, the number of variables, is D=N+1.

The multivariate polynomial of the D variables has the form

    (3) x'_l = x_l + sum_ij...z {DVl_ns n_1^i n_2^j ... r^z}
This form is chosen so that values of zero for the DV terms will reduce to no correction.

The DVl_ns coefficients are ordered as follows. The first D coefficients, starting with 0, are the maximum powers for each of the variables in the order indicated in (3). The next N coefficients are the normalization values for (1). This is followed by the N coefficients in (2). Finally the coefficients in (3) are specified in the order with the first variable varying fastest and the r variable varying slowest. The coefficient numbers for (3) are then related to the polynomial powers by

    (4) n = 3D - 2 + i + j (1+DVl_0s) + k (1+DVl_0s)(1+DVl_1s) + ...

The coefficients have been carefully defined so that a default of zero is sensible. Therefore, any DV keywords not specified take the value zero.

The definition given here is quite general. In practice, only one, two, or three axes are likely to be coupled. Also in many cases one or more of the variables will not be used. When a variable is not used the corresponding maximum power coefficient will be zero or absent.

4.1 A non-linear one dimensional spectrum

Suppose one wants to described the dispersion relation of a non-linear, 1000 point, one-dimensional spectrum (though see Spectral WCS Conventions for a proposal that all spectra be described by a three-dimensional WCS even if the image dimensionality is one or two). The transformations are given by

x_1  = CD1_1 * (p_1 - CRPIX1)
n_1  = x_1 / (1+DV1_2)       
x'_1 = x_1 + DV1_5 n_1^2     
w_1  = CRVAL1 + x'_1         
where p_1 is the pixel coordinate and w_1 is the wavelength and CRVAL1, CRPIX1, and CD1_1 are WCS values given by those keywords. The non-linear polynomial in the example is quadratic. The keywords that describe this WCS are

Figure 5: WCS keywords for a non-linear 1D spectrum.

NAXIS  =                  1 / Number of image raster axes       
NAXIS1 =               1000 / Number of pixels                  
CTYPE1 = 'WAVE-WAV-PLY'     / Wavelength axis                   
CUNIT1 = 'Angstrom'         / Wavelength units                  
CRVAL1 =              4500. / Central wavelength (Angstrom)     
CD1_1  =                2.1 / Linear dispersion (Angstrom/pixel)
CRPIX1 =               500. / Reference pixel                   
DV1_0  =                  2 / Order of non-linear function      
DV1_2  =               500. / Normalization                     
DV1_5  =                 5. / Maximum quadratic deviation       

4.2 A two-dimensional long slit spectrum.

A slit is centered at some celestial coordinate and oriented with some position angle. The slit image is dispersed and imaged on a 2D detector. For simplicity in this example the slit is placed at (12h, 32d) and the length is aligned with right ascension. The pixel scales are approximately 1 arcsecond per pixel along the spatial axis and 10 Angstroms per pixel along the dispersion axis. The detected image is 100x470 with the dispersion along the longer second axis.

We define a three dimensional WCS, one dispersion and two spatial axes, in which one spatial axis is reduced to a single pixel corresponding to the missing third image dimension (see Proposal on FITS WCS Dimensionality for how the missing dimension may be treated). The position of the slit on the sky defines the coordinate reference value and the orientation of the slit on the sky and the pixel scales define the CD terms.

The keywords for this WCS with coupled polynomial corrections is given below.

Figure 6: WCS keywords for a long slit spectrum.

NAXIS   =                     2 / Number of image raster axes       
NAXIS1  =                   100 / Number of pixels                  
NAXIS2  =                   470 / Number of pixels                  
CTYPE1  = 'WAVE-WAV-PLY'       / Wavelength axis with distortion    
CTYPE2  = 'RA---TAN-PLY'       / RA axis with distortion            
CTYPE3  = 'DEC--TAN-PLY'       / DEC axis with distortion           
CRVAL1  =                5560. / Reference wavelength (Angstrom)    
CUNIT1  = 'Angstrom'           / Wavelength unit                    
CRVAL2  =                 180. / Reference RA (deg)                 
CRVAL3  =                  32. / Reference DEC (deg)                
CD1_1   =                  10. / Wavelength scale (Angstrom/pixel)  
CD2_2   =            2.7778E-4 / RA scale (deg)                     
CD3_3   =            2.7778E-4 / DEC scale (deg)                    
CRPIX1  =                 256. / Reference pixel                    
CRPIX2  =                  50. / Reference pixel                    
CRPIX3  =                   1. / Reference pixel                    

DV1_0   =                   2. / 2nd order in wavelength            
DV1_1   =                   2. / 2nd order in RA                    
DV1_4   =                2569. / Wavelength normalization           
DV1_5   =              -0.9858 / RA normalization                   
DV1_12  =                 500. / Wavelength=wavelength^2 coefficient
DV1_16  =                  50. / Wavelength=RA^2 coefficient        

DV2_0   =                   3. / 3rd order in wavelength            
DV2_4   =                2569. / Wavelength normalization           
DV2_12  =            0.0013889 / RA=wavelength^2 coefficient        
DV2_13  =            0.0027778 / RA=wavelength^3 coefficient        

Figure 7 shows a display of the 2D image with lines of constant wavelength and right ascension overlayed. The grid is also in constant intervals of the world coordinates. Because of the non-linear dispersion the vertical (wavelength) grid intervals are non-uniform. Because of distortions and a curved slit the lines of constant spatial position and wavelength are not straight. This example has been exaggerated over a typical long slit spectrum to illustrate the form.

Figure 7: Long slit spectrum with constant coordinate lines.

4.3 Distorted tangent plane image.

One of the uses for the additional r variable is to provide a replacement for the distortion part of the distorted tangent plane system given in section 4.1.2 of the Calabretta and Greisen proposal for celestial coordinates. In other words, the distortion part of the TAN celestial coordinate type is eliminated leaving the standard TAN projection definition. Distortions of the type described there, along with the basis and justifications given for the choice of polynomial, apply to the PLYs function defined here.

An exaggerated example of a "pincushion" distortion is shown below. The form of this simplified distortion is

x'_1 = x_1 + DV1_10 x_1 r = x_1 (1 + DV1_10 r)
x'_2 = x_2 + DV2_10 x_2 r = x_2 (1 + DV2_10 r)
Figure 8 shows the WCS keywords for an exaggerated pincushion distortion and the corresponding lines of constant right ascension and declination in uniform increments are shown overlayed on an image in figure 9.

Figure 8: WCS keywords for an image with pincushion distortion.

NAXIS   =                    2 / Number of image raster axes
NAXIS1  =                  512 / Number of pixels           
NAXIS2  =                  512 / Number of pixels           
CTYPE1  = 'RA---TAN-PLY'       / RA---TAN with distortion   
CTYPE2  = 'DEC--TAN-PLY'       / DEC--TAN with distortion   
CRVAL1  =      201.94541667302 / RA reference (deg)         
CRVAL2  =             47.45444 / DEC reference (deg)        
CD1_1   =        -2.1277777E-4 / RA axis scale (deg/pixel)  
CD2_2   =         2.1277777E-4 / DEC axis scale (deg/pixel) 
CRPIX1  =               257.75 / Reference pixel            
CRPIX2  =               258.93 / Reference pixel            

DV1_0   =                    1 / Up to 1st order in x       
DV1_2   =                    1 / Up to 1st order in r       
DV1_5   =                    1 / Include x in r             
DV1_6   =                    1 / Include y in r             
DV1_10  =                  -3. / xr coefficient             

DV2_1   =                    1 / Up to 1st order in y       
DV2_2   =                    1 / Up to 1st order in r       
DV2_5   =                    1 / Include x in r             
DV2_6   =                    1 / Include y in r             
DV2_10  =                  -3. / yr coefficient             

Figure 9: Image with pincushion distortion.

4.4. Interference filters and Fabry-Perot data.

Imaging through an interference filter, such as occurs when using a scanning Fabry-Perot instrument, produces a spatially dependent wavelength. This is a three-dimensional function in which two of the axes are spatial and the third is spectral.

The basic distortion in wavelength as a function of spatial position is

    (5) w = w(0) cos (arctan (r/C)) = w(0) / sqrt (1 + (r/C)^2)
    (6) w ~ w(0) (1 - 1/2 (r/C)^2 + 3/8 (r/C)^4 ...)
where w(0) is the wavelength at the optical axis. Equation (6) is the power series expansion of (5). Generally r/C<<1 and the power series can be truncated after just a few terms.

The quantity w(0) is related to the intermediate dispersion coordinate x_l by the final ideal projection

    (7) w(0) = CRVALl + x'_l = CRVALl + f(x_l)
So then
    (8) w = CRVALl + f(x_l) - 1/2 CRVALl (r/C)^2 - 1/2 f(x_l) (r/C)^2 ...
    (9) x'_l = f(x_l) - 1/2 CRVALl (r/C)^2 - 1/2 f(x_l) (r/C)^2 ...
For the most common case where f(x_l) = x_l (a linear dispersion) then
    (10) x'_l = x_l - 1/2 CRVALl (r/C)^2 - 1/2 x_l (r/C)^2 ...
which is the PLYs representation with powers of (r^2) and powers of 0 and 1 in x_l. Typically x_l, the deviation from the reference wavelength, is much smaller than CRVALl (for instance it is zero at the reference image plane), and so one could ignore the xr^2 term depending on the accuracy needed.

Note that there are several ways to handle the constant C and the factors of CRVALl. One is simply to fold them into the polynomial coefficients for each power. However the most straightforward way is to normalize the intermediate coordinates by (CRVALl-1) in equation (1) and use (CRVALl/C) as the the radius coefficients in (2). In this case the polynomial coefficients would be the power series coefficients multiplied by CRVALl.

Figures 10 and 11 show an artificial example of a Fabry-Perot data cube. In the example the wavelength steps along the third image dimension are linear but have the basic interference filter behavior as a function of position in the spatial image plane. The spatial coordinates have the same pincushion distortion as shown in the previous example. Figure 10 presents the WCS keywords. This includes 4 terms for illustrations though only the first r^2 term is really significant.

Figure 11 displays two planes of the data cube. The left display is a spatial plane. The lines of constant RA and DEC, in equal increments, are overlayed in yellow. Lines of constant wavelength, also in equal intervals, appear as red circles. The circles are not equally spaced in radius because of the basic r^2 behavior.

The display on the right is a plane with wavelength along the horizontal dimension and RA along the vertical dimension. Lines of constant wavelength and RA, in equal increments, are overlayed. Note the curvature of the wavelength grid which match the night sky spectral lines.

Figure 10: WCS keywords for a Fabry-Perot data cube.

NAXIS   =                     3 / Number of image raster axes      
NAXIS1  =                   512 / Number of pixels                 
NAXIS2  =                   512 / Number of pixels                 
NAXIS3  =                   512 / Number of pixels                 
CTYPE1  = 'RA---TAN-PLY'        / RA---TAN with distortion         
CTYPE2  = 'DEC--TAN-PLY'        / DEC--TAN with distortion         
CTYPE3  = 'WAVE-WAV-PLY'        / WAVE-WAV with interference filter
CRVAL1  =       201.94541667302 / RA reference (deg)               
CRVAL2  =              47.45444 / DEC reference (deg)              
CRVAL3  =                 5560. / Wavelength (Angstrom)            
CUNIT3  = 'Angstrom'            / Units                            
CD1_1   =         -2.1277777E-4 / RA axis scale (deg/pixel)        
CD2_2   =          2.1277777E-4 / DEC axis scale (deg/pixel)       
CD3_3   =                   -1. / Wavelength scale (Angstrom)      
CRPIX1  =                257.75 / Reference pixel                  
CRPIX2  =                258.93 / Reference pixel                  
CRPIX3  =                  256. / Reference pixel                  

DV1_0   =                     1 / Up to 1st order in x             
DV1_3   =                     1 / Up to 1st order in r             
DV1_7   =                     1 / Include x in r                   
DV1_8   =                     1 / Include y in r                   
DV1_13  =                   -3. / xr coefficient                   

DV2_1   =                     1 / Up to 1st order in y             
DV2_3   =                     1 / Up to 1st order in r             
DV2_7   =                     1 / Include x in r                   
DV2_8   =                     1 / Include y in r                   
DV2_13  =                   -3. / yr coefficient                   

DV3_2   =                     1 / Up to 1st order in z             
DV3_3   =                     4 / Up to 4th order in r             
DV3_6   =                 5559. / Normalization (CRVAL-1)          
DV3_7   =                  1.38 / (1/C)                            
DV3_8   =                  1.38 / (1/C)                            
DV3_14  =                -2780. / r^2 coefficient (-0.5*CRVAL)     
DV3_15  =                -2780. / zr^2 coefficient (-0.5*CRVAL)    
DV1_18  =                 2085. / r^4 coefficient (3/8*CRVAL)      
DV1_19  =                 2085. / zr^4 coefficient (3/8*CRVAL)     

Figure 11: Planes of a Fabry-Perot data cube with constant coordinate lines.


This section provides a record of comments and replies received since the public release of this proposal. Comments on this proposal may lead to modifying this document, to clarifications, or may simply be included in this document for the record. Such mail may be edited to include only material relevant to this document and the discussion concerning it. A more complete record of discussions concerning all FITS WCS issues is given elsewhere.