IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10712


Ignore:
Timestamp:
Dec 13, 2006, 8:26:15 PM (20 years ago)
Author:
magnier
Message:

fixing up the 2nd & 3rd-order astrometry

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

Legend:

Unmodified
Added
Removed

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

    r10643 r10712  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2006-12-12 08:00:44 $
     9 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2006-12-14 06:26:15 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    107107
    108108// interpret header WCS keywords (only handles traditional WCS for the moment)
    109 pmAstromWCS *pmAstromWCSfromHeader (psMetadata *header)
     109pmAstromWCS *pmAstromWCSfromHeader (const psMetadata *header)
    110110{
    111111    psProjectionType type;
    … …  
    234234}
    235235
    236 // interpret header WCS (only handles traditional WCS for the moment)
    237 // plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
    238 bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic)
    239 {
    240     psPlaneTransform *toFPA;
    241 
    242     pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
    243     if (!wcs) {
    244         return false;
    245     }
    246 
    247     /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
    248      * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following
    249      * elements defined:
    250      * type (CTYPE)
    251      * crval1,2 (in RA,DEC degrees)
    252      * crpix1,2
    253      * cdelt1,2 (in degrees / pixel)
    254      * plateScale (radians / physical TPA units)
    255      *
    256      * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
    257      * 1) the output is in pixel (not degrees): divide by cdelt1,2 raised to an appropriate power
    258      * 2) X,Y are applied directly, without an applied Xo,Yo offset
    259      * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
    260      */
    261 
    262     // convert wcs->trans to a matrix which yields L,M in pixels
    263     double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
    264     double cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
    265     for (int i = 0; i <= wcs->trans->x->nX; i++) {
    266         for (int j = 0; j <= wcs->trans->x->nX; j++) {
    267             wcs->trans->x->coeff[i][j] /= cdelt1;
    268             wcs->trans->y->coeff[i][j] /= cdelt2;
    269         }
    270     }
    271 
    272     // validate fit order
    273     int fitOrder = wcs->trans->x->nX;
    274     toFPA = psPlaneTransformAlloc(fitOrder, fitOrder);
    275 
    276     /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
    277      * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and
    278      * noting that the constant term in each is the coefficient in the case of L(x,y) and is the
    279      * value of L'(-xo,-yo) in the second case.  in this case, xo,yo = crpix1,2
    280      */
    281 
    282     psPolynomial2D *tmp;
    283 
    284     psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcs->trans->x);
    285     psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcs->trans->y);
    286 
    287     for (int i = 0; i <= fitOrder; i++) {
    288         psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
    289         psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
    290         for (int j = 0; j <= fitOrder; j++) {
    291             toFPA->x->mask[i][j] = wcs->trans->x->mask[i][j];
    292             toFPA->y->mask[i][j] = wcs->trans->y->mask[i][j];
    293             toFPA->x->coeff[i][j] = (toFPA->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
    294             toFPA->y->coeff[i][j] = (toFPA->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
    295 
    296             // take the next derivative wrt y, catch output (is NULL on last pass)
    297             tmp = psPolynomial2D_dY(NULL, xPy);
    298             psFree (xPy);
    299             xPy = tmp;
    300             tmp = psPolynomial2D_dY(NULL, yPy);
    301             psFree (yPy);
    302             yPy = tmp;
    303         }
    304         // take the next derivative wrt x, catch output (is NULL on last pass)
    305         tmp = psPolynomial2D_dX(NULL, xPx);
    306         psFree (xPx);
    307         xPx = tmp;
    308         tmp = psPolynomial2D_dX(NULL, yPx);
    309         psFree (yPx);
    310         yPx = tmp;
    311     }
    312 
    313     // scale from FPA to TPA (microns / pixel)
    314     double pdelt1 = cdelt1*PM_RAD_DEG / plateScale;
    315     double pdelt2 = cdelt2*PM_RAD_DEG / plateScale;
    316     float rX = 1.0;
    317     float rY = 1.0;
    318 
    319     // projection from TPA to SKY
    320     psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, plateScale, plateScale, wcs->toSky->type);
    321 
    322     if (fpa->toSky == NULL) {
    323         fpa->toTPA = psPlaneDistortIdentity (1);
    324         fpa->fromTPA = psPlaneDistortIdentity (1);
    325         fpa->toTPA->x->coeff[1][0][0][0] = pdelt1;
    326         fpa->toTPA->y->coeff[0][1][0][0] = pdelt2;
    327         fpa->fromTPA->x->coeff[1][0][0][0] = 1.0 / pdelt1;
    328         fpa->fromTPA->y->coeff[0][1][0][0] = 1.0 / pdelt2;
    329         fpa->toSky = toSky;
    330     } else {
    331         if (fpa->toTPA == NULL)
    332             psAbort ("wcs", "projection defined, tangent-plane not defined");
    333         if (fpa->fromTPA == NULL)
    334             psAbort ("wcs", "projection defined, tangent-plane not defined");
    335 
    336         // convert from pixels on this chip to pixels on reference chip
    337         // rX has units of refpixels / pixel
    338         rX = pdelt1 / fpa->toTPA->x->coeff[1][0][0][0];
    339         rY = pdelt2 / fpa->toTPA->y->coeff[0][1][0][0];
    340         for (int i = 0; i <= fitOrder; i++) {
    341             for (int j = 0; j <= fitOrder; j++) {
    342                 toFPA->x->coeff[i][j] *= rX;
    343                 toFPA->y->coeff[i][j] *= rY;
    344             }
    345         }
    346 
    347         // adjust for common toSky, toTPA for mosaic:
    348         // find the FPA coordinate of 0,0 for this chip.
    349         psPlane *chip = psPlaneAlloc();
    350         psPlane *fp = psPlaneAlloc();
    351         psPlane *tp = psPlaneAlloc();
    352         psSphere *sky = psSphereAlloc();
    353         chip->x = chip->y = 0;
    354 
    355         psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
    356         psPlaneDistortApply (tp, fpa->toTPA, fp, 0.0, 0.0);
    357         psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
    358         psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
    359         psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0);
    360 
    361         toFPA->x->coeff[0][0] = fp->x;
    362         toFPA->y->coeff[0][0] = fp->y;
    363 
    364         psFree (fp);
    365         psFree (sky);
    366         psFree (chip);
    367         psFree (toSky);
    368     }
    369 
    370     chip->toFPA = toFPA;
    371     // XXX this needs to perform the full (non-linear) inversion
    372     chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA);
    373 
    374     // this can take a very long time...
    375     while (fpa->toSky->R < 0)
    376         fpa->toSky->R += 2.0*M_PI;
    377     while (fpa->toSky->R > 2.0*M_PI)
    378         fpa->toSky->R -= 2.0*M_PI;
    379 
    380     // remove the correction to the common plate scale
    381     // NOTE: this assumes 1) we are reading in headers generated using per-chip astrometry
    382     // and 2) we are going to measure the mosaic distortion in the next step.
    383     // XXX perhaps make this its own function? (I'll need to store rX somewhere).
    384     if (isMosaic) {
    385         chip->toFPA->x->coeff[0][0] /= rX;
    386         chip->toFPA->x->coeff[1][0] /= rX;
    387         chip->toFPA->x->coeff[0][1] /= rX;
    388         chip->toFPA->y->coeff[0][0] /= rY;
    389         chip->toFPA->y->coeff[1][0] /= rY;
    390         chip->toFPA->y->coeff[0][1] /= rY;
    391     }
    392 
    393     psTrace ("psastro", 5, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
    394              chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
    395              chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
    396              chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
    397 
    398     psTrace ("psastro", 5, "frFPA: %f %f  (%f,%f),(%f,%f)\n",
    399              chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
    400              chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
    401              chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
    402 
    403     psFree (wcs->trans);
    404 
    405     return true;
    406 }
    407 
    408 // convert toFPA / toSky components to traditional WCS
    409 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
    410 // this requires toTP to be the identity transformation
    411 bool pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale)
    412 {
    413 
    414     // techinically, we can have a plate scale here (toTPA:dx,dy != 1)
    415     if (!psPlaneDistortIsIdentity (toTPA))
    416         psAbort ("psastro", "invalid TPA transformation");
    417 
    418     // XXX require toFPA->x->nX == toFPA->x->nY
    419     // XXX require toFPA->y->nX == toFPA->y->nY
    420     // XXX require toFPA->x->nX == toFPA->y->nX
    421     // XXX require toFPA->nX == 1,2,3
    422 
    423     switch (toSky->type) {
     236// convert wcs transformations into header WCS keywords (only handles traditional WCS for the moment)
     237// wcs->trans defines the transformation from pixels to degrees.
     238// wcs->cdelt1,2 carries the original pixels scale.
     239// XXX force PC00i00j to be normalized, or use cdelt1,2 to set the scale?
     240// here I've chosen to force the rotation matrix to be normalized
     241bool pmAstromWCStoHeader (psMetadata *header, const pmAstromWCS *wcs)
     242{
     243    char name[16]; // used to store FITS keyword below (always < 8, so 16 should be safe!)
     244
     245    switch (wcs->toSky->type) {
    424246    case PS_PROJ_SIN:
    425247        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---SIN");
    … …  
    439261        break;
    440262    default:
    441         psLogMsg ("psastro", 2, "warning: unknown projection type %d\n", toSky->type);
    442         return false;
    443     }
    444 
    445     # if (0)
    446         // XXX not really right: needs to deal with non-identity to coeffs
    447         // XXX actually, totally wrong.  fix the conversions
    448         // XXX need to handle the plateScale
    449 
    450         /* discussion of the coord transformations:
    451            X,Y: coord on a chip in pixels
    452            L,M: coord on the focal plane (pixels)
    453            P,Q: coord in the tangent plane (microns or mm?)
    454            R,D: coord on the sky
    455          
    456            this function creates WCS terms which convert directly from chip to sky.
    457            this function requires a linear, unrotated toTPA distortion term
    458            toTPA->x,y->coeff[1][0],[0][1] defines the detector scale (microns / pixel)
    459            tpSky->Xs,Ys defines the plate scale (radians / micron)
    460         */
    461 
    462         // solve for CDELT1,2 (degrees / pixel)
    463         cdelt1 = PM_DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
    464     cdelt2 = PM_DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0];
    465 
    466     // L,M = toFPA(X,Y)
    467     // solve for CRPIX1,2 (Xo,Yo) : L,M(Xo,Yo) = 0,0
    468 
    469     // linear solution for Xo,Yo:
    470     xcoeff = toFPA->x->coeff;
    471     ycoeff = toFPA->y->coeff;
    472     R  = (xcoeff[1][0]*ycoeff[0][1] - xcoeff[0][1]*ycoeff[1][0]);
    473     Xo = det*(ycoeff[0][0]*xcoeff[0][1] - xcoeff[0][0]*ycoeff[0][1]);
    474     Yo = det*(xcoeff[0][0]*ycoeff[1][0] - ycoeff[0][0]*xcoeff[1][0]);
    475 
    476     if (toFPA->x->nX > 1) {
    477 
    478         psPolynomial2D *XdX = psPolynomial2D_dX(toFPA->x);
    479         psPolynomial2D *XdY = psPolynomial2D_dY(toFPA->x);
    480 
    481         psPolynomial2D *YdX = psPolynomial2D_dX(toFPA->y);
    482         psPolynomial2D *YdY = psPolynomial2D_dY(toFPA->y);
     263        psLogMsg ("psastro", 2, "warning: unknown projection type %d\n", wcs->toSky->type);
     264        return false;
     265    }
     266
     267    psMetadataAddF64 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", wcs->toSky->R*PM_DEG_RAD);
     268    psMetadataAddF64 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", wcs->toSky->D*PM_DEG_RAD);
     269
     270    psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", wcs->crpix1);
     271    psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", wcs->crpix2);
     272
     273    // XXX make it optional to write out CDi_j terms, or other versions
     274    // solve for CDELT1,2 (degrees / pixel)
     275    double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->y->coeff[1][0]);
     276    double cdelt2 = hypot (wcs->trans->x->coeff[0][1], wcs->trans->y->coeff[0][1]);
     277    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
     278    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
     279
     280    // test the PC00i00j varient:
     281    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
     282    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
     283    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
     284    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
     285
     286    // Elixir-style polynomial terms
     287    // XXX currently, Elixir/DVO cannot accept mixed orders
     288    // XXX need to respect the masks
     289    int fitOrder = wcs->trans->x->nX;
     290    if (fitOrder > 1) {
     291        for (int i = 0; i <= fitOrder; i++) {
     292            for (int j = 0; j <= fitOrder; j++) {
     293                if (i + j < 2)
     294                    continue;
     295                if (i + j > fitOrder)
     296                    continue;
     297                sprintf (name, "PCA1X%1dY%1d", i, j);
     298                psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
     299                sprintf (name, "PCA2X%1dY%1d", i, j);
     300                psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
     301            }
     302        }
     303        psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
     304    }
     305
     306    return (true);
     307}
     308
     309// interpret header WCS (only handles traditional WCS for the moment)
     310// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
     311bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale, bool isMosaic)
     312{
     313    psPlaneTransform *toFPA;
     314
     315    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
     316    if (!wcs) {
     317        return false;
     318    }
     319
     320    /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
     321     * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following
     322     * elements defined:
     323     * type (CTYPE)
     324     * crval1,2 (in RA,DEC degrees)
     325     * crpix1,2
     326     * cdelt1,2 (in degrees / pixel)
     327     * plateScale (radians / physical TPA units)
     328     *
     329     * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
     330     * 1) the output is in pixel (not degrees): divide by cdelt1,2 raised to an appropriate power
     331     * 2) X,Y are applied directly, without an applied Xo,Yo offset
     332     * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
     333     */
     334
     335    // convert wcs->trans to a matrix which yields L,M in pixels
     336    double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     337    double cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
     338    for (int i = 0; i <= wcs->trans->x->nX; i++) {
     339        for (int j = 0; j <= wcs->trans->x->nX; j++) {
     340            wcs->trans->x->coeff[i][j] /= cdelt1;
     341            wcs->trans->y->coeff[i][j] /= cdelt2;
     342        }
     343    }
     344
     345    // validate fit order
     346    int fitOrder = wcs->trans->x->nX;
     347    toFPA = psPlaneTransformAlloc(fitOrder, fitOrder);
     348
     349    /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
     350     * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and
     351     * noting that the constant term in each is the coefficient in the case of L(x,y) and is the
     352     * value of L'(-xo,-yo) in the second case.  in this case, xo,yo = crpix1,2
     353     */
     354
     355    psPolynomial2D *tmp;
     356
     357    psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcs->trans->x);
     358    psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcs->trans->y);
     359
     360    for (int i = 0; i <= fitOrder; i++) {
     361        psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
     362        psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
     363        for (int j = 0; j <= fitOrder; j++) {
     364            toFPA->x->mask[i][j] = wcs->trans->x->mask[i][j];
     365            toFPA->y->mask[i][j] = wcs->trans->y->mask[i][j];
     366            toFPA->x->coeff[i][j] = (toFPA->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
     367            toFPA->y->coeff[i][j] = (toFPA->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1);
     368
     369            // take the next derivative wrt y, catch output (is NULL on last pass)
     370            tmp = psPolynomial2D_dY(NULL, xPy);
     371            psFree (xPy);
     372            xPy = tmp;
     373            tmp = psPolynomial2D_dY(NULL, yPy);
     374            psFree (yPy);
     375            yPy = tmp;
     376        }
     377        // take the next derivative wrt x, catch output (is NULL on last pass)
     378        tmp = psPolynomial2D_dX(NULL, xPx);
     379        psFree (xPx);
     380        xPx = tmp;
     381        tmp = psPolynomial2D_dX(NULL, yPx);
     382        psFree (yPx);
     383        yPx = tmp;
     384    }
     385
     386    // scale from FPA to TPA (microns / pixel)
     387    double pdelt1 = cdelt1*PM_RAD_DEG / plateScale;
     388    double pdelt2 = cdelt2*PM_RAD_DEG / plateScale;
     389    float rX = 1.0;
     390    float rY = 1.0;
     391
     392    // projection from TPA to SKY
     393    psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, plateScale, plateScale, wcs->toSky->type);
     394
     395    if (fpa->toSky == NULL) {
     396        fpa->toTPA = psPlaneDistortIdentity (1);
     397        fpa->fromTPA = psPlaneDistortIdentity (1);
     398        fpa->toTPA->x->coeff[1][0][0][0] = pdelt1;
     399        fpa->toTPA->y->coeff[0][1][0][0] = pdelt2;
     400        fpa->fromTPA->x->coeff[1][0][0][0] = 1.0 / pdelt1;
     401        fpa->fromTPA->y->coeff[0][1][0][0] = 1.0 / pdelt2;
     402        fpa->toSky = toSky;
     403    } else {
     404        if (fpa->toTPA == NULL)
     405            psAbort ("wcs", "projection defined, tangent-plane not defined");
     406        if (fpa->fromTPA == NULL)
     407            psAbort ("wcs", "projection defined, tangent-plane not defined");
     408
     409        // convert from pixels on this chip to pixels on reference chip
     410        // rX has units of refpixels / pixel
     411        rX = pdelt1 / fpa->toTPA->x->coeff[1][0][0][0];
     412        rY = pdelt2 / fpa->toTPA->y->coeff[0][1][0][0];
     413        for (int i = 0; i <= fitOrder; i++) {
     414            for (int j = 0; j <= fitOrder; j++) {
     415                toFPA->x->coeff[i][j] *= rX;
     416                toFPA->y->coeff[i][j] *= rY;
     417            }
     418        }
     419
     420        // adjust for common toSky, toTPA for mosaic:
     421        // find the FPA coordinate of 0,0 for this chip.
     422        psPlane *chip = psPlaneAlloc();
     423        psPlane *fp = psPlaneAlloc();
     424        psPlane *tp = psPlaneAlloc();
     425        psSphere *sky = psSphereAlloc();
     426        chip->x = chip->y = 0;
     427
     428        psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
     429        psPlaneDistortApply (tp, fpa->toTPA, fp, 0.0, 0.0);
     430        psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
     431        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     432        psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0);
     433
     434        toFPA->x->coeff[0][0] = fp->x;
     435        toFPA->y->coeff[0][0] = fp->y;
     436
     437        psFree (fp);
     438        psFree (sky);
     439        psFree (chip);
     440        psFree (toSky);
     441    }
     442
     443    chip->toFPA = toFPA;
     444    // XXX this needs to perform the full (non-linear) inversion
     445    // XXX we need to pull the region from the chip metadata
     446    psRegion region = psRegionSet (0, 4000, 0, 4000);
     447    chip->fromFPA = psPlaneTransformInvert(NULL, toFPA, region, 50);
     448
     449    // this can take a very long time...
     450    while (fpa->toSky->R < 0)
     451        fpa->toSky->R += 2.0*M_PI;
     452    while (fpa->toSky->R > 2.0*M_PI)
     453        fpa->toSky->R -= 2.0*M_PI;
     454
     455    // remove the correction to the common plate scale
     456    // NOTE: this assumes 1) we are reading in headers generated using per-chip astrometry
     457    // and 2) we are going to measure the mosaic distortion in the next step.
     458    // XXX perhaps make this its own function? (I'll need to store rX somewhere).
     459    if (isMosaic) {
     460        chip->toFPA->x->coeff[0][0] /= rX;
     461        chip->toFPA->x->coeff[1][0] /= rX;
     462        chip->toFPA->x->coeff[0][1] /= rX;
     463        chip->toFPA->y->coeff[0][0] /= rY;
     464        chip->toFPA->y->coeff[1][0] /= rY;
     465        chip->toFPA->y->coeff[0][1] /= rY;
     466    }
     467
     468    psTrace ("psastro", 5, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
     469             chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
     470             chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
     471             chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
     472
     473    psTrace ("psastro", 5, "frFPA: %f %f  (%f,%f),(%f,%f)\n",
     474             chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
     475             chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
     476             chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
     477
     478    psFree (wcs);
     479
     480    return true;
     481}
     482
     483// convert toFPA / toSky components to pmAstromWCS, then write to the headers
     484// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     485// this requires toTP to be the identity transformation
     486bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip)
     487{
     488    // techinically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
     489    if (!psPlaneDistortIsIdentity (fpa->toTPA))
     490        psAbort ("psastro", "invalid TPA transformation");
     491
     492    // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
     493    // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
     494    // XXX require chip->toFPA->x->nX == chip->toFPA->y->nX
     495    // XXX require chip->toFPA->nX == 1,2,3
     496
     497    int fitOrder = chip->toFPA->x->nX;
     498    pmAstromWCS *wcs = pmAstromWCSAlloc(fitOrder, fitOrder);
     499
     500    // convert projection from TPA to SKY into wcs projection (degrees to radians)
     501    wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, fpa->toSky->type);
     502    wcs->crval1 = fpa->toSky->R*PM_DEG_RAD;
     503    wcs->crval2 = fpa->toSky->D*PM_DEG_RAD;
     504
     505    // crpix1,2 = X,Y(crval1,2)
     506    // start with linear solution for Xo,Yo:
     507    double R  = (chip->toFPA->x->coeff[1][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->x->coeff[1][0]);
     508    double Xo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1])/R;
     509    double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0])/R;
     510
     511    // iterate to actual solution: requires small non-linear terms
     512    if (fitOrder > 1) {
     513        psPolynomial2D *XdX = psPolynomial2D_dX(NULL, chip->toFPA->x);
     514        psPolynomial2D *XdY = psPolynomial2D_dY(NULL, chip->toFPA->x);
     515
     516        psPolynomial2D *YdX = psPolynomial2D_dX(NULL, chip->toFPA->y);
     517        psPolynomial2D *YdY = psPolynomial2D_dY(NULL, chip->toFPA->y);
    483518
    484519        psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32);
    … …  
    494529            Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
    495530
    496             Beta->data.F32[0] = psPolynomial2DEval (toFPA->x, Xo, Yo);
    497             Beta->data.F32[1] = psPolynomial2DEval (toFPA->y, Xo, Yo);
     531            Beta->data.F32[0] = psPolynomial2DEval (chip->toFPA->x, Xo, Yo);
     532            Beta->data.F32[1] = psPolynomial2DEval (chip->toFPA->y, Xo, Yo);
    498533
    499534            psMatrixGJSolveF32 (Alpha, Beta);
    … …  
    502537            Yo += Beta->data.F32[1];
    503538        }
    504     }
    505 
    506     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", Xo);
    507     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", Yo);
    508 
    509     psPolynomial2D *xWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->x->nX, toFPA->x->nY);
    510     psPolynomial2D *yWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->y->nX, toFPA->y->nY);
    511 
    512     psPolynomial2D *xPx = psPolynomial2DCopy (toFPA->x);
    513     psPolynomial2D *yPx = psPolynomial2DCopy (toFPA->y);
    514 
    515     // skip the zero order terms
    516     // XXX double check that these relationships are correct
    517     for (int i = 0; i < toFPA->x->nX; i++) {
    518         psPolynomial2D *xPy = psPolynomial2DCopy (xPx);
    519         psPolynomial2D *yPy = psPolynomial2DCopy (yPx);
    520         for (int j = 0; j < toFPA->x->nY; j++) {
    521             xWCS->coords[i][j] = psPolynomial2DEval (xPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j);
    522             yWCS->coords[i][j] = psPolynomial2DEval (yPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j);
    523             psPolynomial2D_dY(xPy, xPy);
    524             psPolynomial2D_dY(yPy, yPy);
    525         }
    526         psPolynomial2D_dX(xPx, xPx);
    527         psPolynomial2D_dX(yPx, yPx);
    528     }
    529 
    530     while (coords[0].crval1 < 0)
    531         coords[0].crval1 += 360.0;
    532     while (coords[0].crval1 > 360.0)
    533         coords[0].crval1 -= 360.0;
    534 
    535     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PM_DEG_RAD);
    536     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PM_DEG_RAD);
    537     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", Xo);
    538     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", Yo);
    539     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", cdelt1);
    540     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", cdelt2);
    541 
    542     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", xWCS->coeff[1][0]);
    543     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", xWCS->coeff[0][1]);
    544     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", yWCS->coeff[1][0]);
    545     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", yWCS->coeff[0][1]);
    546 
    547     // XXX respect the masks
    548     for (int i = 0; i < xWCS->nX; i++) {
    549         for (int j = 0; j < xWCS->nX; j++) {
    550             if (i + j < 2)
    551                 continue;
    552             sprintf (name, "PCA1dX%1dY%1d", i, j);
    553             psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", xWCS->coeff[i][j]);
    554         }
    555     }
    556     for (int i = 0; i < yWCS->nX; i++) {
    557         for (int j = 0; j < yWCS->nX; j++) {
    558             if (i + j < 2)
    559                 continue;
    560             sprintf (name, "PCA2dX%1dY%1d", i, j);
    561             psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", yWCS->coeff[i][j]);
    562         }
    563     }
    564 
    565     # else
    566 
    567         psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PM_DEG_RAD);
    568     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PM_DEG_RAD);
    569     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
    570     psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
    571     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PM_DEG_RAD*plateScale);
    572     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PM_DEG_RAD*plateScale);
    573 
    574     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
    575     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
    576     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
    577     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
    578 
    579     # endif
    580 
    581     // alternative representations use
    582     // CD1_1 = PC001001*CDELT1, etc
    583     // make these representations optional
    584 
     539        psFree (Alpha);
     540        psFree (Beta);
     541        psFree (XdX);
     542        psFree (XdY);
     543        psFree (YdX);
     544        psFree (YdY);
     545    }
     546    wcs->crpix1 = Xo;
     547    wcs->crpix2 = Yo;
     548
     549    // convert the chip->toFPA polynomials (with 0,0 ref) into wcs polynomials, with Xo,Yo ref
     550    // chip->toFPA(x,y) = wcs->trans(x-xo,y-yo) -- see comment in pmAstromReadWCS
     551    psPolynomial2D *tmp;
     552    psPolynomial2D *xPx = psPolynomial2DCopy (NULL, chip->toFPA->x);
     553    psPolynomial2D *yPx = psPolynomial2DCopy (NULL, chip->toFPA->y);
     554
     555    for (int i = 0; i <= fitOrder; i++) {
     556        psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
     557        psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
     558        for (int j = 0; j <= fitOrder; j++) {
     559            wcs->trans->x->mask[i][j] = chip->toFPA->x->mask[i][j];
     560            wcs->trans->y->mask[i][j] = chip->toFPA->y->mask[i][j];
     561            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);
     562            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);
     563
     564            // take the next derivative wrt y, catch output (is NULL on last pass)
     565            tmp = psPolynomial2D_dY(NULL, xPy);
     566            psFree (xPy);
     567            xPy = tmp;
     568            tmp = psPolynomial2D_dY(NULL, yPy);
     569            psFree (yPy);
     570            yPy = tmp;
     571        }
     572        // take the next derivative wrt x, catch output (is NULL on last pass)
     573        tmp = psPolynomial2D_dX(NULL, xPx);
     574        psFree (xPx);
     575        xPx = tmp;
     576        tmp = psPolynomial2D_dX(NULL, yPx);
     577        psFree (yPx);
     578        yPx = tmp;
     579    }
     580
     581    // cdelt1,2 convert from pixels->degrees
     582    double cdelt1 = fpa->toTPA->x->coeff[1][0][0][0]*fpa->toSky->Xs*PM_DEG_RAD;
     583    double cdelt2 = fpa->toTPA->y->coeff[0][1][0][0]*fpa->toSky->Ys*PM_DEG_RAD;
     584
     585    // convert wcs->trans to a matrix which yields L,M in pixels
     586    for (int i = 0; i <= wcs->trans->x->nX; i++) {
     587        for (int j = 0; j <= wcs->trans->x->nX; j++) {
     588            wcs->trans->x->coeff[i][j] *= cdelt1;
     589            wcs->trans->y->coeff[i][j] *= cdelt2;
     590        }
     591    }
     592
     593    pmAstromWCStoHeader (header, wcs);
     594
     595    psFree (wcs);
    585596    return true;
    586597}
    … …  
    841852# endif
    842853
     854/* discussion of the coord transformations:
     855   X,Y: coord on a chip in pixels
     856   L,M: coord on the focal plane (pixels)
     857   P,Q: coord in the tangent plane (microns or mm?)
     858   R,D: coord on the sky
     859 
     860   this function creates WCS terms which convert directly from chip to sky.
     861   this function requires a linear, unrotated toTPA distortion term
     862   toTPA->x,y->coeff[1][0],[0][1] defines the detector scale (microns / pixel)
     863   tpSky->Xs,Ys defines the plate scale (radians / micron)
     864*/
     865

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

    r10612 r10712  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-12-10 18:30:07 $
     9*  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-12-14 06:26:15 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    3232bool pmAstromWCStoSky (psSphere *sky, pmAstromWCS *wcs, psPlane *chip);
    3333bool pmAstromWCStoChip (psPlane *chip, pmAstromWCS *wcs, psSphere *sky);
    34 pmAstromWCS *pmAstromWCSfromHeader (psMetadata *header);
     34pmAstromWCS *pmAstromWCSfromHeader (const psMetadata *header);
     35bool pmAstromWCStoHeader (psMetadata *header, const pmAstromWCS *wcs);
    3536
    36 bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic);
    37 bool pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale);
     37bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale, bool isMosaic);
     38bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip);
    3839bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale);
    3940psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale);
Note: See TracChangeset for help on using the changeset viewer.