IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10612


Ignore:
Timestamp:
Dec 10, 2006, 8:30:07 AM (19 years ago)
Author:
magnier
Message:

support for multi-level WCS, higher-order WCS

Location:
trunk/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/pmAstrometryWCS.c

    r10606 r10612  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2006-12-10 04:14:45 $
     9 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2006-12-10 18:30:07 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    2121#include "pmAstrometryWCS.h"
    2222
    23 // interpret header WCS (only handles traditional WCS for the moment)
    24 // plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
    25 bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic)
     23// the following functions support coordinate transformations direcly related to the FITS WCS
     24// keywords.  The FITS WCS allows for only a single level of transformation, thus it is not
     25// appropriate for mosaic astrometry consisting of telescope distortion plus chip terms.
     26// Below, we support the Elixir convention of using two connected FITS headers to define two
     27// levels of coordinate transformation.  In the pmFPA structure, the projection, distortion,
     28// and FPA-to-Chip transformations are carried independently.  NOTE: The FITS WCS keywords do
     29// not represent a simply polynomial.  Instead, they have no constant term, and the coordinates
     30// are corrected to a reference pixel before the polynomial transformation is applied.
     31
     32static void pmAstromWCSFree (pmAstromWCS *wcs)
     33{
     34
     35    if (!wcs)
     36        return;
     37    psFree (wcs->trans);
     38}
     39
     40pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder)
     41{
     42
     43    pmAstromWCS *wcs = (pmAstromWCS *) psAlloc(sizeof(pmAstromWCS));
     44    psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree);
     45
     46    wcs->trans = psPlaneTransformAlloc (nXorder, nYorder);
     47
     48    memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE);
     49    memset (wcs->ctype2, 0, PM_ASTROM_WCS_TYPE_SIZE);
     50    return wcs;
     51}
     52
     53bool pmAstromWCStoSky (psSphere *sky, pmAstromWCS *wcs, psPlane *chip)
     54{
     55
     56    if (chip == NULL)
     57        return false;
     58    if (sky == NULL)
     59        return false;
     60    if (wcs == NULL)
     61        return false;
     62
     63    psPlane *Chip = psPlaneAlloc();
     64    psPlane *FP = psPlaneAlloc();
     65
     66    Chip->x = chip->x - wcs->crpix1;
     67    Chip->y = chip->y - wcs->crpix2;
     68
     69    psPlaneTransformApply (FP, wcs->trans, Chip);
     70    psDeproject (sky, FP, wcs->toSky); // find the RA,DEC coord of the focal-plane coordinate
     71
     72    psFree (Chip);
     73    psFree (FP);
     74    return true;
     75}
     76
     77bool pmAstromWCStoChip (psPlane *chip, pmAstromWCS *wcs, psSphere *sky)
     78{
     79
     80    if (chip == NULL)
     81        return false;
     82    if (sky == NULL)
     83        return false;
     84    if (wcs == NULL)
     85        return false;
     86
     87    psError(PS_ERR_UNKNOWN, true, "not yet implemented: needs to invert the transformation");
     88    return false;
     89
     90    psPlane *Chip = psPlaneAlloc();
     91    psPlane *FP = psPlaneAlloc();
     92
     93    psProject (FP, sky, wcs->toSky); // find the RA,DEC coord of the focal-plane coordinate
     94
     95    // I need the inverse of wcs->transform at this point
     96    psPlaneTransformApply (Chip, wcs->trans, FP);
     97
     98    chip->x = Chip->x + wcs->crpix1;
     99    chip->y = Chip->y + wcs->crpix2;
     100
     101    psFree (Chip);
     102    psFree (FP);
     103    return true;
     104}
     105
     106// interpret header WCS keywords (only handles traditional WCS for the moment)
     107pmAstromWCS *pmAstromWCSfromHeader (psMetadata *header)
    26108{
    27109    psProjectionType type;
    28     psPlaneTransform *toFPA;
    29110    bool status, pcKeys, cdKeys, isPoly;
    30     float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
    31111    char name[16]; // used to store FITS keyword below (always < 8, so 16 should be safe!)
    32112
    … …  
    35115    if (!status) {
    36116        psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n");
    37         return false;
     117        return NULL;
    38118    }
    39119
    … …  
    51131    if (type == PS_PROJ_NTYPE) {
    52132        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);
    53         return false;
     133        return NULL;
    54134    }
    55135
    … …  
    67147        psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms");
    68148        // XXX we could default here to RA, DEC, ROTANGLE
    69         return false;
     149        return NULL;
    70150    }
    71151    if (isPoly && !pcKeys) {
    72152        psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms");
    73         return false;
     153        return NULL;
    74154        if (fitOrder == 0)
    75155            fitOrder = 1;
    76156        if ((fitOrder > 3) || (fitOrder < 1)) {
    77157            psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder);
    78             return false;
     158            return NULL;
    79159        }
    80160    } else {
    81161        fitOrder = 1;
    82162    }
     163
     164    pmAstromWCS *wcs = pmAstromWCSAlloc (fitOrder, fitOrder);
    83165
    84166    // construct a transformation from X,Y in pixels to L,M in pixels
    85167    // NOTE that the WCS keywords convert X,Y to degrees first (using cdelt1,2)
    86168    // and then define a transformation from degrees to degrees
    87     psPlaneTransform *wcsTrans = psPlaneTransformAlloc (fitOrder, fitOrder);
    88 
    89     crval1 = psMetadataLookupF32 (&status, header, "CRVAL1");
    90     crval2 = psMetadataLookupF32 (&status, header, "CRVAL2");
    91     crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1");
    92     crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2");
     169
     170    wcs->crval1 = psMetadataLookupF32 (&status, header, "CRVAL1");
     171    wcs->crval2 = psMetadataLookupF32 (&status, header, "CRVAL2");
     172    wcs->crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1");
     173    wcs->crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2");
     174    wcs->toSky = psProjectionAlloc (wcs->crval1*PM_RAD_DEG, wcs->crval2*PM_RAD_DEG, PM_RAD_DEG, PM_RAD_DEG, type);
    93175
    94176    // test the CDELTi varient
    95177    if (pcKeys) {
    96         cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1");
    97         cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2");
     178        wcs->cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1");
     179        wcs->cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2");
    98180
    99181        // test the CROTAi varient:
    … …  
    101183        double rotate = psMetadataLookupF32 (&status, header, "CROTA2");
    102184        if (status) {
    103             wcsTrans->x->coeff[1][0] = +cdelt1 * cos(rotate*PS_RAD_DEG); // == PC1_1
    104             wcsTrans->x->coeff[0][1] = -cdelt2 * sin(rotate*PS_RAD_DEG); // == PC1_2
    105             wcsTrans->y->coeff[1][0] = +cdelt1 * sin(rotate*PS_RAD_DEG); // == PC2_1
    106             wcsTrans->y->coeff[1][0] = +cdelt2 * cos(rotate*PS_RAD_DEG); // == PC2_2
    107             goto got_matrix;
     185            wcs->trans->x->coeff[1][0] = +wcs->cdelt1 * cos(rotate*PM_RAD_DEG); // == PC1_1
     186            wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2
     187            wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1
     188            wcs->trans->y->coeff[1][0] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
     189            return wcs;
    108190        }
    109191
    110192        // test the PC00i00j varient:
    111         wcsTrans->x->coeff[1][0] = cdelt1 * psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1
    112         wcsTrans->x->coeff[0][1] = cdelt2 * psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2
    113         wcsTrans->y->coeff[1][0] = cdelt1 * psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1
    114         wcsTrans->y->coeff[0][1] = cdelt2 * psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2
     193        wcs->trans->x->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1
     194        wcs->trans->x->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2
     195        wcs->trans->y->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1
     196        wcs->trans->y->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2
    115197
    116198        if (isPoly) {
    117199            // Elixir-style polynomial terms
    118             for (int i = 0; i < wcsTrans->x->nX; i++) {
    119                 for (int j = 0; j < wcsTrans->x->nX; j++) {
     200            // XXX currently, Elixir/DVO cannot accept mixed orders
     201            for (int i = 0; i < fitOrder; i++) {
     202                for (int j = 0; j < fitOrder; j++) {
    120203                    if (i + j < 2)
    121204                        continue;
    … …  
    123206                        continue;
    124207                    sprintf (name, "PCA1dX%1dY%1d", i, j);
    125                     wcsTrans->x->coeff[i][j] = pow(cdelt1, i) * pow(cdelt2, j) * psMetadataLookupF32 (&status, header, name);
     208                    wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name);
     209                    sprintf (name, "PCA2dX%1dY%1d", i, j);
     210                    wcs->trans->y->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name);
    126211                }
    127212            }
    128             // XXX currently, Elixir/DVO cannot accept mixed orders
    129             for (int i = 0; i < wcsTrans->y->nX; i++) {
    130                 for (int j = 0; j < wcsTrans->y->nX; j++) {
    131                     if (i + j < 2)
    132                         continue;
    133                     if (i + j > fitOrder)
    134                         continue;
    135                     sprintf (name, "PCA2dX%1dY%1d", i, j);
    136                     wcsTrans->y->coeff[i][j] = pow(cdelt1, i) * pow(cdelt2, j) * psMetadataLookupF32 (&status, header, name);
    137                 }
    138             }
    139         }
    140         goto got_matrix;
     213        }
     214        return wcs;
    141215    }
    142216
    143217    // test the CDi_j varient
    144218    if (cdKeys) {
    145         wcsTrans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD1_1"); // == PC1_1
    146         wcsTrans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD1_2"); // == PC1_2
    147         wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1
    148         wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2
    149         goto got_matrix;
     219        wcs->cdelt1 = 1.0;
     220        wcs->cdelt2 = 1.0;
     221        wcs->trans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD1_1"); // == PC1_1
     222        wcs->trans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD1_2"); // == PC1_2
     223        wcs->trans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1
     224        wcs->trans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2
     225        return wcs;
    150226    }
    151227    psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n");
    152     return false;
    153 
    154 got_matrix:
     228    return NULL;
     229}
     230
     231// interpret header WCS (only handles traditional WCS for the moment)
     232// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
     233bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic)
     234{
     235    psPlaneTransform *toFPA;
     236
     237    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     238    if (!wcs) {
     239        return false;
     240    }
    155241
    156242    /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
    157      * wcsTrans, which will convert X,Y in pixels to L,M in degrees.  we also have the following
     243     * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following
    158244     * elements defined:
    159245     * type (CTYPE)
    … …  
    163249     * plateScale (radians / physical TPA units)
    164250     *
    165      * now we convert wcsTrans to toFPA, which is different from wcsTrans in 3 important ways:
     251     * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
    166252     * 1) the output is in pixel (not degrees): divide by cdelt1,2 raised to an appropriate power
    167253     * 2) X,Y are applied directly, without an applied Xo,Yo offset
    … …  
    169255     */
    170256
    171     /*** XXXX need to extend these formulae to higher-order terms ***/
    172 
    173     // XXX free an existing toFPA?
     257    // convert wcs->trans to a matrix which yields L,M in pixels
     258    double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     259    double cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
     260    for (int i = 0; i <= wcs->trans->x->nX; i++) {
     261        for (int j = 0; j <= wcs->trans->x->nX; j++) {
     262            wcs->trans->x->coeff[i][j] /= cdelt1;
     263            wcs->trans->y->coeff[i][j] /= cdelt2;
     264        }
     265    }
     266
     267    // validate fit order
     268    int fitOrder = wcs->trans->x->nX;
    174269    toFPA = psPlaneTransformAlloc(fitOrder, fitOrder);
    175270
    176271    /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
    177      * we can transform L(x,y) into L'(x+xo,y+yo) by taking the derivatives of both sides and
     272     * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and
    178273     * noting that the constant term in each is the coefficient in the case of L(x,y) and is the
    179      * value of L'(xo,yo) in the second case.  in this case, xo,yo = crpix1,2
     274     * value of L'(-xo,-yo) in the second case.  in this case, xo,yo = crpix1,2
    180275     */
    181276
    182     psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcsTrans->x);
    183     psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcsTrans->y);
     277    psPolynomial2D *tmp;
     278
     279    psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcs->trans->x);
     280    psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcs->trans->y);
    184281
    185282    for (int i = 0; i <= fitOrder; i++) {
    … …  
    187284        psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
    188285        for (int j = 0; j <= fitOrder; j++) {
    189             toFPA->x->mask[i][j] = wcsTrans->x->mask[i][j];
    190             toFPA->y->mask[i][j] = wcsTrans->y->mask[i][j];
    191             toFPA->x->coeff[i][j] = (toFPA->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, crpix1, crpix2) / tgamma(i+1) / tgamma(j+1) / cdelt1;
    192             toFPA->y->coeff[i][j] = (toFPA->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, crpix1, crpix2) / tgamma(i+1) / tgamma(j+1) / cdelt2;
    193 
    194             // take the next derivative wrt y
    195             psPolynomial2D_dY(xPy, xPy);
    196             psPolynomial2D_dY(yPy, yPy);
    197         }
    198         psFree (xPy);
    199         psFree (yPy);
    200         // take the next derivative wrt x
    201         psPolynomial2D_dX(xPx, xPx);
    202         psPolynomial2D_dX(yPx, yPx);
    203     }
    204     psFree (xPx);
    205     psFree (yPx);
    206 
    207     // save until we verify the transformation
    208     # if (0)
    209         // basic transformation from chip to FPA (FPA in pixels)
    210         toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
    211     toFPA->x->coeff[1][0] = pc1_1;
    212     toFPA->x->coeff[0][1] = pc1_2;
    213     toFPA->x->mask[1][1]  = 1;
    214 
    215     toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2);
    216     toFPA->y->coeff[1][0] = pc2_1;
    217     toFPA->y->coeff[0][1] = pc2_2;
    218     toFPA->y->mask[1][1]  = 1;
    219     # endif
     286            toFPA->x->mask[i][j] = wcs->trans->x->mask[i][j];
     287            toFPA->y->mask[i][j] = wcs->trans->y->mask[i][j];
     288            toFPA->x->coeff[i][j] = (toFPA->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
     289            toFPA->y->coeff[i][j] = (toFPA->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
     290
     291            // take the next derivative wrt y, catch output (is NULL on last pass)
     292            tmp = psPolynomial2D_dY(NULL, xPy);
     293            psFree (xPy);
     294            xPy = tmp;
     295            tmp = psPolynomial2D_dY(NULL, yPy);
     296            psFree (yPy);
     297            yPy = tmp;
     298        }
     299        // take the next derivative wrt x, catch output (is NULL on last pass)
     300        tmp = psPolynomial2D_dX(NULL, xPx);
     301        psFree (xPx);
     302        xPx = tmp;
     303        tmp = psPolynomial2D_dX(NULL, yPx);
     304        psFree (yPx);
     305        yPx = tmp;
     306    }
    220307
    221308    // scale from FPA to TPA (microns / pixel)
    222     double pdelt1 = cdelt1*PS_RAD_DEG / plateScale;
    223     double pdelt2 = cdelt2*PS_RAD_DEG / plateScale;
     309    double pdelt1 = cdelt1*PM_RAD_DEG / plateScale;
     310    double pdelt2 = cdelt2*PM_RAD_DEG / plateScale;
    224311    float rX = 1.0;
    225312    float rY = 1.0;
    226313
    227314    // projection from TPA to SKY
    228     psProjection *toSky = psProjectionAlloc (crval1*PS_RAD_DEG, crval2*PS_RAD_DEG, plateScale, plateScale, type);
     315    psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, plateScale, plateScale, wcs->toSky->type);
    229316
    230317    if (fpa->toSky == NULL) {
    … …  
    263350        psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
    264351        psPlaneDistortApply (tp, fpa->toTPA, fp, 0.0, 0.0);
    265         p_psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
    266         p_psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     352        psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
     353        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
    267354        psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0);
    268355
    … …  
    309396             chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
    310397
    311     psLogMsg ("psastro", 3, "field center: %f,%f, plate scale: %f,%f (arcsec/pixel)\n",
    312               PS_DEG_RAD*fpa->toSky->R, PS_DEG_RAD*fpa->toSky->D,
    313               3600*PS_DEG_RAD*fpa->toSky->Xs, 3600*PS_DEG_RAD*fpa->toSky->Ys);
    314 
    315     psFree (wcsTrans);
     398    psFree (wcs->trans);
    316399
    317400    return true;
    … …  
    373456
    374457        // solve for CDELT1,2 (degrees / pixel)
    375         cdelt1 = PS_DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
    376     cdelt2 = PS_DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0];
     458        cdelt1 = PM_DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
     459    cdelt2 = PM_DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0];
    377460
    378461    // L,M = toFPA(X,Y)
    … …  
    445528        coords[0].crval1 -= 360.0;
    446529
    447     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
    448     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
     530    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PM_DEG_RAD);
     531    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PM_DEG_RAD);
    449532    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", Xo);
    450533    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", Yo);
    … …  
    477560    # else
    478561
    479         psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
    480     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
     562        psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PM_DEG_RAD);
     563    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PM_DEG_RAD);
    481564    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
    482565    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
    483     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale);
    484     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale);
     566    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PM_DEG_RAD*plateScale);
     567    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PM_DEG_RAD*plateScale);
    485568
    486569    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
    … …  
    566649
    567650    // XXX need to handle the plateScale??
    568     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
    569     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
     651    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PM_DEG_RAD);
     652    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PM_DEG_RAD);
    570653    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", 0.0);
    571654    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", 0.0);
    572     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale);
    573     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale);
     655    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PM_DEG_RAD*plateScale);
     656    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PM_DEG_RAD*plateScale);
    574657
    575658    if (toTP->x->nX != toTP->x->nY)
    … …  
    738821*****/
    739822
     823
     824// save until we verify the transformation
     825# if (0)
     826    // basic transformation from chip to FPA (FPA in pixels)
     827    toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
     828toFPA->x->coeff[1][0] = pc1_1;
     829toFPA->x->coeff[0][1] = pc1_2;
     830toFPA->x->mask[1][1]  = 1;
     831
     832toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2);
     833toFPA->y->coeff[1][0] = pc2_1;
     834toFPA->y->coeff[0][1] = pc2_2;
     835toFPA->y->mask[1][1]  = 1;
     836# endif
     837

  • trunk/psModules/src/astrom/pmAstrometryWCS.h

    r10603 r10612  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-12-10 02:06:47 $
     9*  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-12-10 18:30:07 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    1515#ifndef PM_ASTROMETRY_WCS_H
    1616#define PM_ASTROMETRY_WCS_H
     17
     18#define PM_ASTROM_WCS_TYPE_SIZE 80
     19typedef struct
     20{
     21    char ctype1[PM_ASTROM_WCS_TYPE_SIZE];
     22    char ctype2[PM_ASTROM_WCS_TYPE_SIZE];
     23    double crval1, crval2;
     24    double crpix1, crpix2;
     25    double cdelt1, cdelt2;
     26    psProjection *toSky;
     27    psPlaneTransform *trans;
     28}
     29pmAstromWCS;
     30
     31pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder);
     32bool pmAstromWCStoSky (psSphere *sky, pmAstromWCS *wcs, psPlane *chip);
     33bool pmAstromWCStoChip (psPlane *chip, pmAstromWCS *wcs, psSphere *sky);
     34pmAstromWCS *pmAstromWCSfromHeader (psMetadata *header);
    1735
    1836bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic);
    … …  
    2543bool psPlaneDistortIsIdentity (psPlaneDistort *distort);
    2644
    27 # define PS_DEG_RAD 57.295779513082322
    28 # define PS_RAD_DEG  0.017453292519943
     45# define PM_DEG_RAD 57.295779513082322
     46# define PM_RAD_DEG  0.017453292519943
    2947
    3048#endif // PM_ASTROMETRY_WCS_H
     49
     50/*
     51 * the wcs->trans component defines a polynomial which converts (x-crpix1),(y-crpix2) to
     52 * L,M in degrees
     53 */

  • trunk/psModules/test/astrom/tap_pmAstrometryWCS.c

    r10606 r10612  
    8383    }
    8484
     85    {
     86        diag("test pmAstromReadWCS");
     87        psMemId id = psMemGetId();
     88
     89        psPlaneTransform *wcsTrans = psPlaneTransformAlloc (1, 1);
     90        wcsTrans->x->coeff[0][0] = 0.0;
     91        wcsTrans->y->coeff[0][0] = 0.0;
     92        wcsTrans->x->coeff[1][0] = 1.0;
     93        wcsTrans->x->coeff[0][1] = 0.0;
     94        wcsTrans->y->coeff[1][0] = 0.0;
     95        wcsTrans->y->coeff[0][1] = 1.0;
     96
     97
     98
     99        // construct a header with a simple set of WCS values
     100        // convert to pmFPA components and check
     101
     102        psMetadata *header = psMetadataAlloc();
     103
     104        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN");
     105        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--TAN");
     106
     107        // center coords (R,D)
     108        psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", 0.0);
     109        psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", 0.0);
     110
     111        // center coords (X,Y)
     112        psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", 10.0);
     113        psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", 10.0);
     114
     115        // degrees per pixel
     116        psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", 1.0/3600.0);
     117        psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", 1.0/3600.0);
     118
     119        // rotation matrix
     120        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", 1.0);
     121        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", 0.0);
     122        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", 0.0);
     123        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", 1.0);
     124
     125        pmFPA *fpa = pmFPAAlloc (NULL);
     126        pmChip *chip = pmChipAlloc (NULL, "test");
     127
     128        // toFPA carries pixel scale (pixels per micron)
     129        // toTPA carries plate scale (microns per arcsecond)
     130        bool status = pmAstromReadWCS (fpa, chip, header, 25.0*PS_RAD_DEG/3600.0, false);
     131        ok (status, "converted WCS keywords to FPA astrometry");
     132        skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     133
     134        ok(fpa->toSky->type == PS_PROJ_TAN, "correct projection (TAN)");
     135
     136        // make these tests double
     137        ok_float(fpa->toSky->R*PS_DEG_RAD, 0.0, "projection center RA %f", fpa->toSky->R*PS_DEG_RAD);
     138        ok_float(fpa->toSky->R*PS_DEG_RAD, 0.0, "projection center DEC %f", fpa->toSky->R*PS_DEG_RAD);
     139
     140        ok_float(fpa->toSky->Xs, 25.0*PS_RAD_DEG/3600.0, "projection X scale %f", fpa->toSky->Xs);
     141        ok_float(fpa->toSky->Ys, 25.0*PS_RAD_DEG/3600.0, "projection X scale %f", fpa->toSky->Ys);
     142
     143        ok_float(fpa->toTPA->x->coeff[1][0][0][0], 0.04, "TP scale (mm per pixel): %f", fpa->toTPA->x->coeff[1][0][0][0]);
     144        ok_float(fpa->toTPA->y->coeff[0][1][0][0], 0.04, "TP scale (mm per pixel): %f", fpa->toTPA->x->coeff[1][0][0][0]);
     145
     146        ok_float(chip->toFPA->x->coeff[0][0], -10.0, "ref pixel X: %f", chip->toFPA->x->coeff[0][0]);
     147        ok_float(chip->toFPA->y->coeff[0][0], -10.0, "ref pixel Y: %f", chip->toFPA->y->coeff[0][0]);
     148
     149        ok_float(chip->toFPA->x->coeff[1][0], 1.0, "CD1_1: %f", chip->toFPA->x->coeff[1][0]);
     150        ok_float(chip->toFPA->x->coeff[0][1], 0.0, "CD1_2: %f", chip->toFPA->x->coeff[0][1]);
     151        ok_float(chip->toFPA->y->coeff[1][0], 0.0, "CD2_1: %f", chip->toFPA->y->coeff[1][0]);
     152        ok_float(chip->toFPA->y->coeff[0][1], 1.0, "CD2_2: %f", chip->toFPA->y->coeff[0][1]);
     153
     154        // apply both systems to real data:
     155
     156
     157        psFree (fpa);
     158        psFree (chip);
     159        psFree (header);
     160
     161        skip_end();
     162        ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     163    }
     164
    85165    return exit_status();
    86166}
Note: See TracChangeset for help on using the changeset viewer.