IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10825


Ignore:
Timestamp:
Dec 22, 2006, 11:23:06 AM (19 years ago)
Author:
magnier
Message:

cleaned up WCS conversions, added bilevel astrometry read/write functions

Location:
trunk/psModules/src/astrom
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/Makefile.am

    r10599 r10825  
    66        pmAstrometryObjects.c \
    77        pmAstrometryDistortion.c \
     8        pmAstrometryUtils.c \
    89        pmAstrometryWCS.c
    910
    … …  
    1112        pmAstrometryObjects.h \
    1213        pmAstrometryDistortion.h \
     14        pmAstrometryUtils.h \
    1315        pmAstrometryWCS.h
    1416

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

    r10603 r10825  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-12-10 02:06:47 $
     9*  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-12-22 21:23:06 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    187187    for (int i = 0; i < fpa->toTPA->x->nX; i++) {
    188188        for (int j = 0; j < fpa->toTPA->x->nY; j++) {
    189             if (fpa->toTPA->x->mask[i][j][0][0]) {
     189            if (fpa->toTPA->x->mask[i][j]) {
    190190                localX->mask[i-1][j] = 1;
    191191                localY->mask[i][j-1] = 1;
    … …  
    203203
    204204    // update fpa->toTP distortion terms
    205     fpa->toTPA->x->coeff[0][0][0][0] = 0;
     205    fpa->toTPA->x->coeff[0][0] = 0;
    206206    for (int i = 1; i < fpa->toTPA->x->nX; i++) {
    207         if (!fpa->toTPA->x->mask[i][0][0][0]) {
    208             fpa->toTPA->x->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
     207        if (!fpa->toTPA->x->mask[i][0]) {
     208            fpa->toTPA->x->coeff[i][0] = localX->coeff[i-1][0] / i;
    209209        }
    210210    }
    211211    for (int j = 1; j < fpa->toTPA->x->nY; j++) {
    212         if (!fpa->toTPA->x->mask[0][j][0][0]) {
    213             fpa->toTPA->x->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
     212        if (!fpa->toTPA->x->mask[0][j]) {
     213            fpa->toTPA->x->coeff[0][j] = localY->coeff[0][j-1] / j;
    214214        }
    215215    }
    216216    for (int i = 1; i < fpa->toTPA->x->nX; i++) {
    217217        for (int j = 1; j < fpa->toTPA->x->nY; j++) {
    218             if (!fpa->toTPA->x->mask[i][j][0][0]) {
    219                 fpa->toTPA->x->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
     218            if (!fpa->toTPA->x->mask[i][j]) {
     219                fpa->toTPA->x->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
    220220            }
    221221        }
    … …  
    232232    for (int i = 0; i < fpa->toTPA->y->nX; i++) {
    233233        for (int j = 0; j < fpa->toTPA->y->nY; j++) {
    234             if (fpa->toTPA->y->mask[i][j][0][0]) {
     234            if (fpa->toTPA->y->mask[i][j]) {
    235235                localX->mask[i-1][j] = 1;
    236236                localY->mask[i][j-1] = 1;
    … …  
    259259
    260260    // update fpa->toTP distortion terms
    261     fpa->toTPA->y->coeff[0][0][0][0] = 0;
     261    fpa->toTPA->y->coeff[0][0] = 0;
    262262    for (int i = 1; i < fpa->toTPA->y->nX; i++) {
    263         if (!fpa->toTPA->y->mask[i][0][0][0]) {
    264             fpa->toTPA->y->coeff[i][0][0][0] = localX->coeff[i-1][0] / i;
     263        if (!fpa->toTPA->y->mask[i][0]) {
     264            fpa->toTPA->y->coeff[i][0] = localX->coeff[i-1][0] / i;
    265265        }
    266266    }
    267267    for (int j = 1; j < fpa->toTPA->y->nY; j++) {
    268         if (!fpa->toTPA->y->mask[0][j][0][0]) {
    269             fpa->toTPA->y->coeff[0][j][0][0] = localY->coeff[0][j-1] / j;
     268        if (!fpa->toTPA->y->mask[0][j]) {
     269            fpa->toTPA->y->coeff[0][j] = localY->coeff[0][j-1] / j;
    270270        }
    271271    }
    272272    for (int i = 1; i < fpa->toTPA->y->nX; i++) {
    273273        for (int j = 1; j < fpa->toTPA->y->nY; j++) {
    274             if (!fpa->toTPA->y->mask[i][j][0][0]) {
    275                 fpa->toTPA->y->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
     274            if (!fpa->toTPA->y->mask[i][j]) {
     275                fpa->toTPA->y->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
    276276            }
    277277        }

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

    r10814 r10825  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2006-12-20 16:37:55 $
     9 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2006-12-22 21:23:06 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    2020#include "pmFPA.h"
    2121#include "pmAstrometryWCS.h"
     22#include "pmAstrometryUtils.h"
    2223
    2324// the following functions support coordinate transformations direcly related to the FITS WCS
    … …  
    2728// levels of coordinate transformation.  In the pmFPA structure, the projection, distortion,
    2829// 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// not represent a simple polynomial.  Instead, they have no constant term, and the coordinates
    3031// are corrected to a reference pixel before the polynomial transformation is applied.
    3132
    32 static void pmAstromWCSFree (pmAstromWCS *wcs)
    33 {
    34 
    35     if (!wcs)
    36         return;
    37     psFree (wcs->trans);
    38     psFree (wcs->toSky);
    39 }
    40 
    41 pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder)
    42 {
    43 
    44     pmAstromWCS *wcs = (pmAstromWCS *) psAlloc(sizeof(pmAstromWCS));
    45     psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree);
    46 
    47     wcs->trans = psPlaneTransformAlloc (nXorder, nYorder);
    48     wcs->toSky = NULL;
    49 
    50     memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE);
    51     memset (wcs->ctype2, 0, PM_ASTROM_WCS_TYPE_SIZE);
    52     return wcs;
    53 }
    54 
     33// interpret header WCS (only handles traditional WCS for the moment)
     34// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
     35bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale)
     36{
     37    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     38    if (!wcs) {
     39        return false;
     40    }
     41
     42    bool status = pmAstromWCStoFPA (fpa, chip, wcs, plateScale);
     43
     44    psFree (wcs);
     45    return status;
     46}
     47
     48// convert toFPA / toSky components to pmAstromWCS
     49// tolerance is convergence for inversion of non-linear terms in pixels
     50bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol)
     51{
     52    // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
     53    // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
     54    // XXX require chip->toFPA->x->nX == chip->toFPA->y->nX
     55    // XXX require chip->toFPA->nX == 1,2,3
     56
     57    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
     58    if (!psPlaneTransformIsDiagonal (fpa->toTPA))
     59        psAbort ("psastro", "invalid TPA transformation");
     60
     61    pmAstromWCS *wcs = pmAstromWCSfromFPA(fpa, chip, tol);
     62    pmAstromWCStoHeader (header, wcs);
     63
     64    psFree (wcs);
     65    return true;
     66}
     67
     68// interpret chip header WCS as bilevel chip components
     69// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
     70bool pmAstromReadBilevelChip (pmChip *chip, const psMetadata *header)
     71{
     72    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     73    if (!wcs) {
     74        return false;
     75    }
     76
     77    bool status = pmAstromWCSBileveltoChip (chip, wcs);
     78
     79    psFree (wcs);
     80    return status;
     81}
     82
     83// convert toFPA / toSky components to traditional WCS
     84// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     85bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header)
     86{
     87    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     88    if (!wcs) {
     89        return false;
     90    }
     91
     92    bool status = pmAstromWCSBileveltoFPA (fpa, wcs);
     93
     94    psFree (wcs);
     95    return status;
     96}
     97
     98// convert chip->toFPA components to bilevel WCS
     99bool pmAstromWriteBilevelChip (psMetadata *header, const pmChip *chip, float tol)
     100{
     101    pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (chip, tol);
     102
     103    pmAstromWCStoHeader (header, wcs);
     104
     105    psFree (wcs);
     106    return true;
     107}
     108
     109
     110// convert fpa->toTPA, fpa->toSky components to bilevel WCS
     111bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, float tol)
     112{
     113    pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol);
     114
     115    pmAstromWCStoHeader (header, wcs);
     116
     117    psFree (wcs);
     118    return true;
     119}
     120
     121// convert coordinates from chip to sky using a pmAstromWCS structure
    55122bool pmAstromWCStoSky (psSphere *sky, pmAstromWCS *wcs, psPlane *chip)
    56123{
    … …  
    77144}
    78145
     146// convert coordinates from sky to chip using a pmAstromWCS structure
    79147bool pmAstromWCStoChip (psPlane *chip, pmAstromWCS *wcs, psSphere *sky)
    80148{
    … …  
    106174}
    107175
    108 // interpret header WCS keywords (only handles traditional WCS for the moment)
     176// interpret header WCS keywords (valid for bilevel and traditional WCS)
    109177pmAstromWCS *pmAstromWCSfromHeader (const psMetadata *header)
    110178{
    … …  
    121189
    122190    // determine projection type
    123     // XXX add the Elixir DIS / WRP two-layer projection here
    124     type = PS_PROJ_NTYPE;
    125     if (!strcmp (&ctype[4], "-SIN"))
    126         type = PS_PROJ_SIN;
    127     if (!strcmp (&ctype[4], "-TAN"))
    128         type = PS_PROJ_TAN;
    129     if (!strcmp (&ctype[4], "-AIT"))
    130         type = PS_PROJ_AIT;
    131     if (!strcmp (&ctype[4], "-PAR"))
    132         type = PS_PROJ_PAR;
     191    // XXX there are two indications for higher-order terms: the type (DIS,WRP,PLY,ZPL) and
     192    // the value of NPLYTERM.
     193    type = psProjectTypeFromString (ctype);
    133194    if (type == PS_PROJ_NTYPE) {
    134195        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);
    … …  
    245306{
    246307    char name[16]; // used to store FITS keyword below (always < 8, so 16 should be safe!)
    247 
    248     switch (wcs->toSky->type) {
    249     case PS_PROJ_SIN:
    250         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---SIN");
    251         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--SIN");
    252         break;
    253     case PS_PROJ_TAN:
    254         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN");
    255         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--TAN");
    256         break;
    257     case PS_PROJ_AIT:
    258         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---AIT");
    259         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--AIT");
    260         break;
    261     case PS_PROJ_PAR:
    262         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---PAR");
    263         psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--PAR");
    264         break;
    265     default:
    266         psLogMsg ("psastro", 2, "warning: unknown projection type %d\n", wcs->toSky->type);
    267         return false;
    268     }
     308    char *type;
     309
     310    type = psProjectTypeToString (wcs->toSky->type, "RA--");
     311    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", type);
     312    psFree (type);
     313
     314    type = psProjectTypeToString (wcs->toSky->type, "DEC-");
     315    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", type);
     316    psFree (type);
    269317
    270318    psMetadataAddF64 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", wcs->toSky->R*PM_DEG_RAD);
    … …  
    312360// interpret header WCS (only handles traditional WCS for the moment)
    313361// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
    314 bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale, bool isMosaic)
     362bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale)
    315363{
    316364    psPlaneTransform *toFPA;
    … …  
    326374     *
    327375     * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
    328      * 1) the output is in pixel (not degrees): divide by cdelt1,2 raised to an appropriate power
     376     * 1) the output is in pixels (not degrees): divide by cdelt1,2 raised to an appropriate power
    329377     * 2) X,Y are applied directly, without an applied Xo,Yo offset
    330378     * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
    … …  
    341389    }
    342390
    343     // validate fit order
    344     int fitOrder = wcs->trans->x->nX;
    345     toFPA = psPlaneTransformAlloc(fitOrder, fitOrder);
    346 
    347     /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
    348      * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and
    349      * noting that the constant term in each is the coefficient in the case of L(x,y) and is the
    350      * value of L'(-xo,-yo) in the second case.  in this case, xo,yo = crpix1,2
    351      */
    352 
    353     psPolynomial2D *tmp;
    354 
    355     psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcs->trans->x);
    356     psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcs->trans->y);
    357 
    358     for (int i = 0; i <= fitOrder; i++) {
    359         psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
    360         psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
    361         for (int j = 0; j <= fitOrder; j++) {
    362             toFPA->x->mask[i][j] = wcs->trans->x->mask[i][j];
    363             toFPA->y->mask[i][j] = wcs->trans->y->mask[i][j];
    364             toFPA->x->coeff[i][j] = (toFPA->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
    365             toFPA->y->coeff[i][j] = (toFPA->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
    366 
    367             // take the next derivative wrt y, catch output (is NULL on last pass)
    368             tmp = psPolynomial2D_dY(NULL, xPy);
    369             psFree (xPy);
    370             xPy = tmp;
    371             tmp = psPolynomial2D_dY(NULL, yPy);
    372             psFree (yPy);
    373             yPy = tmp;
    374         }
    375         // take the next derivative wrt x, catch output (is NULL on last pass)
    376         tmp = psPolynomial2D_dX(NULL, xPx);
    377         psFree (xPx);
    378         xPx = tmp;
    379         tmp = psPolynomial2D_dX(NULL, yPx);
    380         psFree (yPx);
    381         yPx = tmp;
    382     }
     391    // create transformation with 0,0 reference pixel
     392    toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
    383393
    384394    // scale from FPA to TPA (microns / pixel)
    … …  
    392402
    393403    if (fpa->toSky == NULL) {
    394         fpa->toTPA = psPlaneDistortIdentity (1);
    395         fpa->fromTPA = psPlaneDistortIdentity (1);
    396         fpa->toTPA->x->coeff[1][0][0][0] = pdelt1;
    397         fpa->toTPA->y->coeff[0][1][0][0] = pdelt2;
    398         fpa->fromTPA->x->coeff[1][0][0][0] = 1.0 / pdelt1;
    399         fpa->fromTPA->y->coeff[0][1][0][0] = 1.0 / pdelt2;
     404        fpa->toTPA = psPlaneTransformIdentity (1);
     405        fpa->fromTPA = psPlaneTransformIdentity (1);
     406        fpa->toTPA->x->coeff[1][0] = pdelt1;
     407        fpa->toTPA->y->coeff[0][1] = pdelt2;
     408        fpa->fromTPA->x->coeff[1][0] = 1.0 / pdelt1;
     409        fpa->fromTPA->y->coeff[0][1] = 1.0 / pdelt2;
    400410        fpa->toSky = toSky;
    401411    } else {
     412
     413        // this section allows the loaded chip to be included in an fpa structure in which
     414        // other chips have already been loaded (ie, the fpa->toTPA, fpa->toSky components have
     415        // already been defined).  we have to adjust to match the existing transformation.
     416
    402417        if (fpa->toTPA == NULL)
    403418            psAbort ("wcs", "projection defined, tangent-plane not defined");
    … …  
    407422        // convert from pixels on this chip to pixels on reference chip
    408423        // rX has units of refpixels / pixel
    409         rX = pdelt1 / fpa->toTPA->x->coeff[1][0][0][0];
    410         rY = pdelt2 / fpa->toTPA->y->coeff[0][1][0][0];
    411         for (int i = 0; i <= fitOrder; i++) {
    412             for (int j = 0; j <= fitOrder; j++) {
     424        rX = pdelt1 / fpa->toTPA->x->coeff[1][0];
     425        rY = pdelt2 / fpa->toTPA->y->coeff[0][1];
     426        for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
     427            for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
    413428                toFPA->x->coeff[i][j] *= rX;
    414429                toFPA->y->coeff[i][j] *= rY;
    … …  
    425440
    426441        psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
    427         psPlaneDistortApply (tp, fpa->toTPA, fp, 0.0, 0.0);
     442        psPlaneTransformApply (tp, fpa->toTPA, fp);
    428443        psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
    429444        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
    430         psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0);
     445        psPlaneTransformApply (fp, fpa->fromTPA, tp);
    431446
    432447        toFPA->x->coeff[0][0] = fp->x;
    … …  
    454469        fpa->toSky->R -= 2.0*M_PI;
    455470
    456     // remove the correction to the common plate scale
    457     // NOTE: this assumes 1) we are reading in headers generated using per-chip astrometry
    458     // and 2) we are going to measure the mosaic distortion in the next step.
    459     // XXX perhaps make this its own function? (I'll need to store rX somewhere).
    460     if (isMosaic) {
    461         chip->toFPA->x->coeff[0][0] /= rX;
    462         chip->toFPA->x->coeff[1][0] /= rX;
    463         chip->toFPA->x->coeff[0][1] /= rX;
    464         chip->toFPA->y->coeff[0][0] /= rY;
    465         chip->toFPA->y->coeff[1][0] /= rY;
    466         chip->toFPA->y->coeff[0][1] /= rY;
    467     }
    468 
    469471    psTrace ("psastro", 5, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
    470472             chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
    … …  
    480482}
    481483
     484// convert a pmAstromWCS structure representing a bilevel chip into corresponding chip elements
     485bool pmAstromWCSBileveltoChip (pmChip *chip, const pmAstromWCS *wcs)
     486{
     487    /* we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
     488     * 1) the output is in pixel (not degrees): divide by cdelt1,2 raised to an appropriate power
     489     * 2) X,Y are applied directly, without an applied Xo,Yo offset
     490     * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
     491     */
     492
     493    // create transformation with 0,0 reference pixel
     494    chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
     495
     496    // XXX this needs to perform the full (non-linear) inversion
     497    // XXX we need to pull the region from the chip metadata
     498    // is chip trimsec defined?  do we need to ensure the
     499    // XXX should this function be defined for the CELL, not the CHIP?
     500    // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC");
     501    psRegion region = psRegionSet (0, 4000, 0, 4000);
     502    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, region, 50);
     503
     504    return true;
     505}
     506
     507// convert a pmAstromWCS structure representing a bilevel mosaic into corresponding fpa elements
     508bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs)
     509{
     510    // projection from TPA to SKY
     511    fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->toSky->Xs, wcs->toSky->Ys, wcs->toSky->type);
     512
     513    // create transformation with 0,0 reference pixel
     514    fpa->toTPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
     515
     516    // XXX this needs to perform the full (non-linear) inversion
     517    // XXX we need to pull the region from the chip metadata
     518    // is chip trimsec defined?  do we need to ensure the
     519    // XXX should this function be defined for the CELL, not the CHIP?
     520    // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC");
     521    psRegion region = psRegionSet (0, 4000, 0, 4000);
     522    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 50);
     523
     524    return true;
     525}
     526
    482527// convert toFPA / toSky components to pmAstromWCS
    483 // tolerance is in pixels
     528// tolerance is allowed error in center solution in pixels
    484529pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, float tol)
    485530{
    486531
    487532    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
    488     if (!psPlaneDistortIsDiagonal (fpa->toTPA))
     533    if (!psPlaneTransformIsDiagonal (fpa->toTPA))
    489534        psAbort ("psastro", "invalid TPA transformation");
    490535
    … …  
    494539    // XXX require chip->toFPA->nX == 1,2,3
    495540
    496     int fitOrder = chip->toFPA->x->nX;
    497     pmAstromWCS *wcs = pmAstromWCSAlloc(fitOrder, fitOrder);
     541    pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
    498542
    499543    // convert projection from TPA to SKY into wcs projection (degrees to radians)
    … …  
    502546    wcs->crval2 = fpa->toSky->D*PM_DEG_RAD;
    503547
    504     // crpix1,2 = X,Y(crval1,2)
    505     // start with linear solution for Xo,Yo:
    506     double R  = (chip->toFPA->x->coeff[1][0]*chip->toFPA->y->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->y->coeff[1][0]);
    507     double Xo = (chip->toFPA->y->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[0][1])/R;
    508     double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[1][0] - chip->toFPA->y->coeff[0][0]*chip->toFPA->x->coeff[1][0])/R;
    509 
    510     // iterate to actual solution: requires small non-linear terms
    511     if (fitOrder > 1) {
    512         psPolynomial2D *XdX = psPolynomial2D_dX(NULL, chip->toFPA->x);
    513         psPolynomial2D *XdY = psPolynomial2D_dY(NULL, chip->toFPA->x);
    514 
    515         psPolynomial2D *YdX = psPolynomial2D_dX(NULL, chip->toFPA->y);
    516         psPolynomial2D *YdY = psPolynomial2D_dY(NULL, chip->toFPA->y);
    517 
    518         psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32);
    519         psVector *Beta = psVectorAlloc (2, PS_DATA_F32);
    520 
    521         // XXX this loop is rather arbitrary in length...
    522         // XXX measure the error and use as criterion
    523         /* this is the Newton-Raphson method to solve for Xo,Yo - it needs the high order terms to be small */
    524         // Xo,Yo are in pixels;
    525         float dPos = tol + 1;
    526         for (int i = 0; (dPos > tol) && (i < 10); i++) {
    527             // NOTE: order is: [y][x]
    528             Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo);
    529             Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
    530             Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
    531             Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
    532 
    533             Beta->data.F32[0] = psPolynomial2DEval (chip->toFPA->x, Xo, Yo);
    534             Beta->data.F32[1] = psPolynomial2DEval (chip->toFPA->y, Xo, Yo);
    535 
    536             psMatrixGJSolveF32 (Alpha, Beta);
    537 
    538             Xo -= Beta->data.F32[0];
    539             Yo -= Beta->data.F32[1];
    540             dPos = hypot(Xo,Yo);
    541         }
    542         psFree (Alpha);
    543         psFree (Beta);
    544         psFree (XdX);
    545         psFree (XdY);
    546         psFree (YdX);
    547         psFree (YdY);
    548     }
    549     wcs->crpix1 = Xo;
    550     wcs->crpix2 = Yo;
    551 
    552     // convert the chip->toFPA polynomials (with 0,0 ref) into wcs polynomials, with Xo,Yo ref
    553     // chip->toFPA(x,y) = wcs->trans(x-xo,y-yo) -- see comment in pmAstromReadWCS
    554     psPolynomial2D *tmp;
    555     psPolynomial2D *xPx = psPolynomial2DCopy (NULL, chip->toFPA->x);
    556     psPolynomial2D *yPx = psPolynomial2DCopy (NULL, chip->toFPA->y);
    557 
    558     for (int i = 0; i <= fitOrder; i++) {
    559         psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
    560         psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
    561         for (int j = 0; j <= fitOrder; j++) {
    562             wcs->trans->x->mask[i][j] = chip->toFPA->x->mask[i][j];
    563             wcs->trans->y->mask[i][j] = chip->toFPA->y->mask[i][j];
    564             wcs->trans->x->coeff[i][j] = (wcs->trans->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, wcs->crpix1, wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
    565             wcs->trans->y->coeff[i][j] = (wcs->trans->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, wcs->crpix1, wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
    566 
    567             // take the next derivative wrt y, catch output (is NULL on last pass)
    568             tmp = psPolynomial2D_dY(NULL, xPy);
    569             psFree (xPy);
    570             xPy = tmp;
    571             tmp = psPolynomial2D_dY(NULL, yPy);
    572             psFree (yPy);
    573             yPy = tmp;
    574         }
    575         // take the next derivative wrt x, catch output (is NULL on last pass)
    576         tmp = psPolynomial2D_dX(NULL, xPx);
    577         psFree (xPx);
    578         xPx = tmp;
    579         tmp = psPolynomial2D_dX(NULL, yPx);
    580         psFree (yPx);
    581         yPx = tmp;
    582     }
     548    // given transformation, solve for coordinates which yields output coordinates of 0,0
     549    psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
     550
     551    // adjust wcs transform to use center as reference coordinate
     552    psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
     553
     554    // calculated center is crpix1,2
     555    wcs->crpix1 = center->x;
     556    wcs->crpix2 = center->y;
     557    psFree (center);
    583558
    584559    // cdelt1,2 convert from pixels->degrees
    585     double cdelt1 = fpa->toTPA->x->coeff[1][0][0][0]*fpa->toSky->Xs*PM_DEG_RAD;
    586     double cdelt2 = fpa->toTPA->y->coeff[0][1][0][0]*fpa->toSky->Ys*PM_DEG_RAD;
     560    double cdelt1 = fpa->toTPA->x->coeff[1][0]*fpa->toSky->Xs*PM_DEG_RAD;
     561    double cdelt2 = fpa->toTPA->y->coeff[0][1]*fpa->toSky->Ys*PM_DEG_RAD;
    587562    wcs->cdelt1 = cdelt1;
    588563    wcs->cdelt2 = cdelt2;
    … …  
    598573}
    599574
    600 // interpret header WCS (only handles traditional WCS for the moment)
    601 // plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
    602 bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale, bool isMosaic)
    603 {
    604     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
    605     if (!wcs) {
    606         return false;
    607     }
    608 
    609     bool status = pmAstromWCStoFPA (fpa, chip, wcs, plateScale, isMosaic);
    610 
    611     psFree (wcs);
    612     return status;
    613 }
    614 
    615 // convert toFPA / toSky components to pmAstromWCS, then write to the headers
    616 // tolerance is in pixels
    617 bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol)
     575/* the bilevel astrometry description consists of a polynomial warping from
     576   chip coordinates to FPA coordinates (coords->ctype = LIN---WRP), followed
     577   by a polynomial representation of the telescope distortion + the projection
     578   (coords->ctype = RA---DIS).
     579*/
     580
     581// convert the chip-level toFPA to a wcs polynomial transformation
     582pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, float tol)
    618583{
    619584    // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
    … …  
    622587    // XXX require chip->toFPA->nX == 1,2,3
    623588
    624     // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
    625     if (!psPlaneDistortIsDiagonal (fpa->toTPA))
    626         psAbort ("psastro", "invalid TPA transformation");
    627 
    628     pmAstromWCS *wcs = pmAstromWCSfromFPA(fpa, chip, tol);
    629     pmAstromWCStoHeader (header, wcs);
    630 
    631     psFree (wcs);
    632     return true;
    633 }
    634 
    635 // convert toFPA / toSky components to traditional WCS
    636 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
    637 bool pmAstromWriteBilevelChip (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol)
    638 {
    639     pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (fpa, chip, tol);
    640 
    641     pmAstromWCStoHeader (header, wcs);
    642 
    643     psFree (wcs);
    644     return true;
    645 }
    646 
    647 
    648 // convert toFPA / toSky components to traditional WCS
    649 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
    650 bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, float tol)
    651 {
    652     pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol);
    653 
    654     pmAstromWCStoHeader (header, wcs);
    655 
    656     psFree (wcs);
    657     return true;
    658 }
    659 
    660 pmAstromWCS *pmAstromWCSBilevelChipFromFPA (pmFPA *fpa, float tol)
    661 {
    662 
    663     // conv
    664 
    665     psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---WRP");
    666     psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--WRP");
    667 
    668     // XXX not really right: needs to deal with non-identity toTP coeffs & plateScale
    669     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
    670     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
    671     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", 0.0);
    672     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", 0.0);
    673     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", 1.0);
    674     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", 1.0);
    675 
    676     if (toFPA->x->nX != toFPA->x->nY)
    677         psAbort ("psastro", "mis-matched tangent plane orders (1)");
    678     if (toFPA->x->nX != toFPA->y->nX)
    679         psAbort ("psastro", "mis-matched tangent plane orders (2)");
    680     if (toFPA->x->nX != toFPA->y->nY)
    681         psAbort ("psastro", "mis-matched tangent plane orders (3)");
    682 
    683     switch (toFPA->x->nX) {
    684     case 3:
    685         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toFPA->x->coeff[3][0]);
    686         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toFPA->x->coeff[2][1]);
    687         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toFPA->x->coeff[1][2]);
    688         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toFPA->x->coeff[0][3]);
    689 
    690         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toFPA->y->coeff[3][0]);
    691         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toFPA->y->coeff[2][1]);
    692         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toFPA->y->coeff[1][2]);
    693         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toFPA->y->coeff[0][3]);
    694 
    695     case 2:
    696         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toFPA->x->coeff[2][0]);
    697         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toFPA->x->coeff[1][1]);
    698         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toFPA->x->coeff[0][2]);
    699 
    700         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toFPA->y->coeff[2][0]);
    701         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toFPA->y->coeff[1][1]);
    702         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toFPA->y->coeff[0][2]);
    703 
    704     case 1:
    705         psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
    706         psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
    707         psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
    708         psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
    709         break;
    710 
    711     case 0:
    712         psAbort ("psastro", "invalid tangent plane order");
    713     }
    714     return true;
    715 }
    716 
    717 // convert toFPA / toSky components to traditional WCS
    718 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
    719 psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale)
    720 {
    721 
    722     psMetadata *header = psMetadataAlloc ();
    723 
    724     psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---DIS");
    725     psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--DIS");
    726 
    727     // XXX need to handle the plateScale??
    728     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PM_DEG_RAD);
    729     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PM_DEG_RAD);
    730     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", 0.0);
    731     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", 0.0);
    732     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PM_DEG_RAD*plateScale);
    733     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PM_DEG_RAD*plateScale);
    734 
    735     if (toTP->x->nX != toTP->x->nY)
    736         psAbort ("psastro", "mis-matched tangent plane orders (1)");
    737     if (toTP->x->nX != toTP->y->nX)
    738         psAbort ("psastro", "mis-matched tangent plane orders (2)");
    739     if (toTP->x->nX != toTP->y->nY)
    740         psAbort ("psastro", "mis-matched tangent plane orders (3)");
    741 
    742     switch (toTP->x->nX) {
    743     case 3:
    744         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toTP->x->coeff[3][0][0][0]);
    745         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toTP->x->coeff[2][1][0][0]);
    746         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toTP->x->coeff[1][2][0][0]);
    747         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toTP->x->coeff[0][3][0][0]);
    748 
    749         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toTP->y->coeff[3][0][0][0]);
    750         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toTP->y->coeff[2][1][0][0]);
    751         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toTP->y->coeff[1][2][0][0]);
    752         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toTP->y->coeff[0][3][0][0]);
    753 
    754     case 2:
    755         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toTP->x->coeff[2][0][0][0]);
    756         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toTP->x->coeff[1][1][0][0]);
    757         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toTP->x->coeff[0][2][0][0]);
    758 
    759         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toTP->y->coeff[2][0][0][0]);
    760         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toTP->y->coeff[1][1][0][0]);
    761         psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toTP->y->coeff[0][2][0][0]);
    762 
    763     case 1:
    764         psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toTP->x->coeff[1][0][0][0]);
    765         psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toTP->x->coeff[0][1][0][0]);
    766         psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toTP->y->coeff[1][0][0][0]);
    767         psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toTP->y->coeff[0][1][0][0]);
    768         break;
    769 
    770     case 0:
    771         psAbort ("psastro", "invalid tangent plane order");
    772     }
    773     return header;
    774 }
    775 
    776 // construct a psPlaneDistort which is the identify transformation
    777 psPlaneDistort *psPlaneDistortIdentity (int order)
    778 {
    779 
    780     psPlaneDistort *distort;
    781 
    782     if (order < 1)
    783         psAbort ("psastro", "invalid order");
    784     if (order > 3)
    785         psAbort ("psastro", "invalid order");
    786 
    787     // all coeffs and masks initially set to 0
    788     distort = psPlaneDistortAlloc (order, order, 0, 0);
    789 
    790     for (int i = 0; i <= order; i++) {
    791         for (int j = 0; j <= order; j++) {
    792             if (i + j > order) {
    793                 distort->x->mask [i][j][0][0] = 1;
    794                 distort->y->mask [i][j][0][0] = 1;
    795             }
    796         }
    797     }
    798     distort->x->coeff[1][0][0][0] = 1;
    799     distort->y->coeff[0][1][0][0] = 1;
    800 
    801     return distort;
    802 }
    803 
    804 // check that the given psPlaneDistort is the identity * (Xs,Ys)
    805 bool psPlaneDistortIsDiagonal (psPlaneDistort *distort)
    806 {
    807 
    808     int order;
    809     bool status;
    810 
    811     // we currently only support up to 3rd order polynomials
    812     if (distort->x->nX < 1)
    813         return false;
    814     if (distort->x->nY < 1)
    815         return false;
    816     if (distort->y->nX < 1)
    817         return false;
    818     if (distort->y->nY < 1)
    819         return false;
    820 
    821     if (distort->x->nX > 3)
    822         return false;
    823     if (distort->x->nY > 3)
    824         return false;
    825     if (distort->y->nX > 3)
    826         return false;
    827     if (distort->y->nY > 3)
    828         return false;
    829 
    830     if (distort->x->nZ > 0)
    831         return false;
    832     if (distort->x->nT > 0)
    833         return false;
    834     if (distort->y->nZ > 0)
    835         return false;
    836     if (distort->y->nT > 0)
    837         return false;
    838 
    839     if (distort->x->nX != distort->x->nY)
    840         return false;
    841     if (distort->y->nX != distort->y->nY)
    842         return false;
    843 
    844     status = true;
    845     order = distort->x->nX;
    846     for (int i = 0; i <= order; i++) {
    847         for (int j = 0; j <= order; j++) {
    848             if (i + j > order) {
    849                 // high-order cross terms must be masked (eg, x^3 y^2)
    850                 status &= distort->x->mask[i][j][0][0];
    851             } else {
    852                 status &= !distort->x->mask[i][j][0][0];
    853                 if ((i == 1) && (i + j == 1)) {
    854                     // linear, diagonal terms must be 1.0
    855                     status &= (fabs(distort->x->coeff[i][j][0][0]) > FLT_EPSILON);
    856                 } else {
    857                     // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
    858                     status &= (fabs(distort->x->coeff[i][j][0][0]) < FLT_EPSILON);
    859                 }
    860             }
    861         }
    862     }
    863 
    864     order = distort->y->nX;
    865     for (int i = 0; i <= order; i++) {
    866         for (int j = 0; j <= order; j++) {
    867             if (i + j > order) {
    868                 // high-order cross terms must be masked (eg, x^3 y^2)
    869                 status &= distort->y->mask[i][j][0][0];
    870             } else {
    871                 status &= !distort->y->mask[i][j][0][0];
    872                 if ((j == 1) && (i + j == 1)) {
    873                     // linear, diagonal terms must be 1.0
    874                     status &= (fabs(distort->y->coeff[i][j][0][0]) > FLT_EPSILON);
    875                 } else {
    876                     // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
    877                     status &= (fabs(distort->y->coeff[i][j][0][0]) < FLT_EPSILON);
    878                 }
    879             }
    880         }
    881     }
    882     return status;
     589    // convert chip->toFPA to wcs format (WRP)
     590    pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
     591
     592    // Chip to FPA transformation is a Cartesian 'projection'
     593    wcs->toSky = psProjectionAlloc (0.0, 0.0, 1.0, 1.0, PS_PROJ_WRP);
     594
     595    // reference pixel for FPA is 0.0, 0.0
     596    wcs->crval1 = 0.0;
     597    wcs->crval2 = 0.0;
     598
     599    // given transformation, solve for coordinates which yields output coordinates of 0,0
     600    psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
     601
     602    // adjust wcs transform to use center as reference coordinate
     603    psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
     604
     605    // calculated center is crpix1,2
     606    wcs->crpix1 = center->x;
     607    wcs->crpix2 = center->y;
     608    psFree (center);
     609
     610    // output coordinates are in FPA pixels
     611    wcs->cdelt1 = 1.0;
     612    wcs->cdelt2 = 1.0;
     613
     614    return wcs;
     615}
     616
     617// convert the fpa-level toTPA, toSky to a wcs polynomial transformation
     618pmAstromWCS *pmAstromWCSBilevelMosaicFromFPA (const pmFPA *fpa, float tol)
     619{
     620    // XXX require fpa->toTPA->x->nX == fpa->toTPA->x->nY
     621    // XXX require fpa->toTPA->y->nX == fpa->toTPA->y->nY
     622    // XXX require fpa->toTPA->x->nX == fpa->toTPA->y->nX
     623    // XXX require fpa->toTPA->nX == 1,2,3
     624    // XXX require fpa->toSky->type == PS_PROJ_TAN
     625
     626    // convert fpa->toTPA + fpa->toSky to wcs format (DIS)
     627    pmAstromWCS *wcs = pmAstromWCSAlloc(fpa->toTPA->x->nX, fpa->toTPA->x->nY);
     628
     629    // convert projection from TPA to SKY into wcs projection (degrees to radians)
     630    wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, PS_PROJ_DIS);
     631    wcs->crval1 = fpa->toSky->R*PM_DEG_RAD;
     632    wcs->crval2 = fpa->toSky->D*PM_DEG_RAD;
     633
     634    // given transformation, solve for coordinates which yields output coordinates of 0,0
     635    psPlane *center = psPlaneTransformGetCenter (fpa->toTPA, tol);
     636
     637    // adjust wcs transform to use center as reference coordinate
     638    psPlaneTransformSetCenter (wcs->trans, fpa->toTPA, center->x, center->y);
     639
     640    // calculated center is crpix1,2
     641    wcs->crpix1 = center->x;
     642    wcs->crpix2 = center->y;
     643    psFree (center);
     644
     645    // output coordinates are in FPA pixels
     646    wcs->cdelt1 = 1.0;
     647    wcs->cdelt2 = 1.0;
     648
     649    return wcs;
     650}
     651
     652static void pmAstromWCSFree (pmAstromWCS *wcs)
     653{
     654
     655    if (!wcs)
     656        return;
     657    psFree (wcs->trans);
     658    psFree (wcs->toSky);
     659}
     660
     661pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder)
     662{
     663
     664    pmAstromWCS *wcs = (pmAstromWCS *) psAlloc(sizeof(pmAstromWCS));
     665    psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree);
     666
     667    wcs->trans = psPlaneTransformAlloc (nXorder, nYorder);
     668    wcs->toSky = NULL;
     669
     670    memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE);
     671    memset (wcs->ctype2, 0, PM_ASTROM_WCS_TYPE_SIZE);
     672    return wcs;
    883673}
    884674
    … …  
    900690*****/
    901691
    902 
    903 // save until we verify the transformation
    904 # if (0)
    905     // basic transformation from chip to FPA (FPA in pixels)
    906     toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
    907 toFPA->x->coeff[1][0] = pc1_1;
    908 toFPA->x->coeff[0][1] = pc1_2;
    909 toFPA->x->mask[1][1]  = 1;
    910 
    911 toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2);
    912 toFPA->y->coeff[1][0] = pc2_1;
    913 toFPA->y->coeff[0][1] = pc2_2;
    914 toFPA->y->mask[1][1]  = 1;
    915 # endif
    916 
    917692/* discussion of the coord transformations:
    918693   X,Y: coord on a chip in pixels

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

    r10783 r10825  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-12-17 09:46:56 $
     9*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-12-22 21:23:06 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    2929pmAstromWCS;
    3030
    31 // read wcs terms from the supplied header into the fpa hierarchy components
    32 bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale, bool isMosaic);
    33 
    34 // write the wcs terms from the fpa hierarchy components into the supplied header
    35 // tol is the convergence tolerance for the non-linear solution to the reference pixel
    36 bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol);
    37 
    3831// support function for the pmAstromWCS representation
    3932pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder);
    … …  
    4538pmAstromWCS *pmAstromWCSfromHeader (const psMetadata *header);
    4639
     40// convert from wcs terms to chip->toFPA, fpa->toSky,toTPA terms
     41bool pmAstromWCSBileveltoChip (pmChip *chip, const pmAstromWCS *wcs);
     42bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs);
     43
     44// convert from chip->toFPA, fpa->toSky,toTPA terms to wcs terms
     45pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, float tol);
     46pmAstromWCS *pmAstromWCSBilevelMosaicFromFPA (const pmFPA *fpa, float tol);
     47
    4748// convert the pmAstromWCS representation to the FPA representation
    48 bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale, bool isMosaic);
     49bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale);
    4950pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, float tol);
    5051
    51 bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale);
    52 psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale);
     52// read wcs terms from the supplied header into the fpa hierarchy components
     53bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale);
     54
     55// write the wcs terms from the fpa hierarchy components into the supplied header
     56// tol is the convergence tolerance for the non-linear solution to the reference pixel
     57bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol);
     58
     59bool pmAstromReadBilevelChip (pmChip *chip, const psMetadata *header);
     60bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header);
     61
     62bool pmAstromWriteBilevelChip (psMetadata *header, const pmChip *chip, float tol);
     63bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, float tol);
    5364
    5465// move to pslib
Note: See TracChangeset for help on using the changeset viewer.