IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 4, 2006, 7:02:16 AM (20 years ago)
Author:
eugene
Message:

substantial updates

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastroWCS.c

    r7014 r7332  
    33// interpret header WCS (only handles traditional WCS for the moment)
    44// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
    5 bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale) {
    6 
    7     psPlaneTransform *toFPA;
    8     psPlaneDistort *toTPA;
    9     psProjection *toSky;
     5bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) {
     6
    107    psProjectionType type;
    118    bool status;
    … …  
    1310    float pc1_1, pc1_2, pc2_1, pc2_2;
    1411
    15     // *** interpret header data, convert to crval(i), etc
     12    // interpret header data, convert to crval(i), etc
    1613    char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");
    1714    if (!status) {
    … …  
    8885got_matrix:
    8986
    90     // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the
    91     //           toFPA must be modified to match the crval(i), scale(i) and crpix(i)
    92     //           scale = scale(i)/scale(0) (i == chip #)
    93     //           project crval1(0),crval2(0 using projection
    94     toFPA = psPlaneTransformAlloc (1, 1);
     87    /*****
     88
     89    For mosaic astrometry, we need to have a starting set of projection terms in which the
     90    chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
     91    microns).  This set of projections, coupled with an identity toTPA (ie, no distortion) will
     92    result in substantial errors between the observed and predicted star positions on the focal
     93    plane: this is the measurement of the optical distortion in the camera.  At the same time,
     94    we need to carry around the transformations which allow us to make an accurate calculation
     95    of the position of the stars based on the input (per-chip) astrometry.  These
     96    transformations will allow us to match the raw and ref stars robustly.  To convert the
     97    per-chip astrometry (which may have been calculated with a different plate scale for each
     98    chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust
     99    the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
     100    each chip (eg, cdelt1).  Thus, we need to carry around both the
     101
     102    *****/
     103
     104    {
     105        double rX = 1.0;
     106        double rY = 1.0;
     107
     108        // XXX free an existing toFPA?
     109        psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1);
    95110   
    96     double cdelt = hypot (cdelt1, cdelt2) / plateScale;  // degrees / micron (eg, in fact, whatever unit user chooses for focal plane)
    97     cdelt1 /= cdelt;
    98     cdelt2 /= cdelt;
    99 
    100     toFPA->x->coeff[0][0] = -(pc1_1*cdelt1*crpix1 + pc1_2*cdelt2*crpix2);
    101     toFPA->x->coeff[1][0] = +(pc1_1*cdelt1);
    102     toFPA->x->coeff[0][1] = +(pc1_2*cdelt2);
    103     toFPA->x->mask[1][1]  = 1;
    104 
    105     toFPA->y->coeff[0][0] = -(pc2_1*cdelt1*crpix1 + pc2_2*cdelt2*crpix2);
    106     toFPA->y->coeff[1][0] = +(pc2_1*cdelt1);
    107     toFPA->y->coeff[0][1] = +(pc2_2*cdelt2);
    108     toFPA->y->mask[1][1]  = 1;
    109     *toFPAout = toFPA;
    110 
    111     if (*toSkyOut != NULL) {
    112         if (*toTPAout == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
    113 
    114         psSphere sky;
    115         psPlane tpa;
    116 
    117         sky.r = crval1*RAD_DEG;
    118         sky.d = crval2*RAD_DEG;
    119         p_psProject (&tpa, &sky, *toSkyOut);
     111        // basic transformation from chip to FPA
     112        toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
     113        toFPA->x->coeff[1][0] = pc1_1;
     114        toFPA->x->coeff[0][1] = pc1_2;
     115        toFPA->x->mask[1][1]  = 1;
     116
     117        toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2);
     118        toFPA->y->coeff[1][0] = pc2_1;
     119        toFPA->y->coeff[0][1] = pc2_2;
     120        toFPA->y->mask[1][1]  = 1;
     121
     122        // projection from TPA to SKY
     123        psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
     124
     125        if (fpa->toSky == NULL) {
     126            // XXX for now, use the identity for TPA <--> FPA
     127            fpa->toTangentPlane = psPlaneDistortIdentity (1);
     128            fpa->fromTangentPlane = psPlaneDistortIdentity (1);
     129            fpa->toSky = toSky;
     130        } else {
     131            if (fpa->toTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
     132            if (fpa->fromTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
     133
     134            // adjust for common toSky for mosaic:
     135            // ignore the TPA since toTPA is identity
     136            // find the FPA coordinate of 0,0 for this chip.
     137            psPlane *fp = psPlaneAlloc();
     138            psPlane *chip = psPlaneAlloc();
     139            psSphere *sky = psSphereAlloc();
     140            chip->x = chip->y = 0;
    120141       
    121         // XXX for the moment, assume toTPA is the identity transformation
    122         toFPA->x->coeff[0][0] = tpa.x;
    123         toFPA->y->coeff[0][0] = tpa.y;
    124     } else {
    125         // XXX for now, use the identity for TPA <--> FPA
    126         toTPA = psPlaneDistortIdentity ();
    127 
    128         // center of projection is (0,0) coordinate of TPA
    129         toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, RAD_DEG*cdelt, RAD_DEG*cdelt, type);
    130 
    131         *toTPAout = toTPA;
    132         *toSkyOut = toSky;
     142            psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
     143            p_psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate
     144            p_psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     145
     146            toFPA->x->coeff[0][0] = fp->x;
     147            toFPA->y->coeff[0][0] = fp->y;
     148
     149            rX = toSky->Xs/fpa->toSky->Xs;
     150            rY = toSky->Ys/fpa->toSky->Ys;
     151
     152            // correct to common plate scale
     153            toFPA->x->coeff[1][0] *= rX;
     154            toFPA->x->coeff[0][1] *= rX;
     155            toFPA->y->coeff[1][0] *= rY;
     156            toFPA->y->coeff[0][1] *= rY;
     157
     158            psFree (fp);
     159            psFree (sky);
     160            psFree (chip);
     161            psFree (toSky);
     162        }
     163
     164        chip->toFPA = toFPA;
     165        chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA);
     166
     167        // remove the correction to the common plate scale
     168        if (isMosaic) {
     169            chip->toFPA->x->coeff[1][0] /= rX;
     170            chip->toFPA->x->coeff[0][1] /= rX;
     171            chip->toFPA->y->coeff[1][0] /= rY;
     172            chip->toFPA->y->coeff[0][1] /= rY;
     173        }
     174
     175        fprintf (stderr, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
     176                 chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
     177                 chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
     178                 chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
     179
     180        fprintf (stderr, "frFPA: %f %f  (%f,%f),(%f,%f)\n",
     181                 chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
     182                 chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
     183                 chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
     184
     185        fprintf (stderr, "field: %f,%f  %f,%f\n",
     186                 DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D,
     187                 3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys);
     188
    133189    }
    134190    return true;
    135191}
    136 
    137192
    138193// convert toFPA / toSky components to traditional WCS
    … …  
    169224    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
    170225    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
    171     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
    172     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
    173 
    174     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
    175     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
    176     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
    177     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
     226    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
     227    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
     228
     229    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
     230    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
     231    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
     232    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
    178233
    179234    // alternative representations use
    … …  
    184239}
    185240
    186 psPlaneDistort *psPlaneDistortIdentity () {
     241// convert toFPA / toSky components to traditional WCS
     242// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     243bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale) {
     244
     245    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---WRP");
     246    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--WRP");
     247
     248    // XXX not really right: needs to deal with non-identity toTP coeffs & plateScale
     249    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",   PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
     250    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",   PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
     251    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", 0.0);
     252    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", 0.0);
     253    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", 1.0);
     254    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", 1.0);
     255
     256    if (toFPA->x->nX != toFPA->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)");
     257    if (toFPA->x->nX != toFPA->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)");
     258    if (toFPA->x->nX != toFPA->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)");
     259
     260    switch (toFPA->x->nX) {
     261      case 3:
     262        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toFPA->x->coeff[3][0]);
     263        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toFPA->x->coeff[2][1]);
     264        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toFPA->x->coeff[1][2]);
     265        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toFPA->x->coeff[0][3]);
     266
     267        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toFPA->y->coeff[3][0]);
     268        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toFPA->y->coeff[2][1]);
     269        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toFPA->y->coeff[1][2]);
     270        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toFPA->y->coeff[0][3]);
     271
     272      case 2:
     273        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toFPA->x->coeff[2][0]);
     274        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toFPA->x->coeff[1][1]);
     275        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toFPA->x->coeff[0][2]);
     276
     277        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toFPA->y->coeff[2][0]);
     278        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toFPA->y->coeff[1][1]);
     279        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toFPA->y->coeff[0][2]);
     280
     281      case 1:
     282        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
     283        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
     284        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
     285        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
     286        break;
     287
     288      case 0:
     289        psAbort ("psastro", "invalid tangent plane order");
     290    }
     291    return true;
     292}
     293
     294// convert toFPA / toSky components to traditional WCS
     295// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     296psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale) {
     297
     298    psMetadata *header = psMetadataAlloc ();
     299
     300    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---DIS");
     301    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--DIS");
     302
     303    // XXX need to handle the plateScale??
     304    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",   PS_META_REPLACE, "", toSky->R*DEG_RAD);
     305    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",   PS_META_REPLACE, "", toSky->D*DEG_RAD);
     306    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", 0.0);
     307    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", 0.0);
     308    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
     309    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
     310
     311    if (toTP->x->nX != toTP->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)");
     312    if (toTP->x->nX != toTP->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)");
     313    if (toTP->x->nX != toTP->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)");
     314
     315    switch (toTP->x->nX) {
     316      case 3:
     317        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toTP->x->coeff[3][0][0][0]);
     318        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toTP->x->coeff[2][1][0][0]);
     319        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toTP->x->coeff[1][2][0][0]);
     320        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toTP->x->coeff[0][3][0][0]);
     321
     322        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toTP->y->coeff[3][0][0][0]);
     323        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toTP->y->coeff[2][1][0][0]);
     324        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toTP->y->coeff[1][2][0][0]);
     325        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toTP->y->coeff[0][3][0][0]);
     326
     327      case 2:
     328        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toTP->x->coeff[2][0][0][0]);
     329        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toTP->x->coeff[1][1][0][0]);
     330        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toTP->x->coeff[0][2][0][0]);
     331
     332        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toTP->y->coeff[2][0][0][0]);
     333        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toTP->y->coeff[1][1][0][0]);
     334        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toTP->y->coeff[0][2][0][0]);
     335
     336      case 1:
     337        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toTP->x->coeff[1][0][0][0]);
     338        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toTP->x->coeff[0][1][0][0]);
     339        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toTP->y->coeff[1][0][0][0]);
     340        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toTP->y->coeff[0][1][0][0]);
     341        break;
     342
     343      case 0:
     344        psAbort ("psastro", "invalid tangent plane order");
     345    }
     346    return header;
     347}
     348
     349// XXX for the moment, we mimic the limited Elixir mosastro orders
     350psPlaneDistort *psPlaneDistortIdentity (int order) {
    187351
    188352    psPlaneDistort *distort;
    189353
    190     distort = psPlaneDistortAlloc (1, 1, 0, 0);
    191     distort->x->coeff[0][0][0][0] = 0;
     354    if (order < 1) psAbort ("psastro", "invalid order");
     355    if (order > 3) psAbort ("psastro", "invalid order");
     356   
     357    // all coeffs and masks initially set to 0
     358    distort = psPlaneDistortAlloc (order, order, 0, 0);
     359
     360    for (int i = 0; i <= order; i++) {
     361        for (int j = 0; j <= order; j++) {
     362            if (i + j > order) {
     363                distort->x->mask [1][1][0][0] = 1;
     364                distort->y->mask [1][1][0][0] = 1;
     365            }
     366        }
     367    }
    192368    distort->x->coeff[1][0][0][0] = 1;
    193     distort->x->coeff[0][1][0][0] = 0;
    194     distort->x->mask [1][1][0][0] = 1;
    195 
    196     distort->y->coeff[0][0][0][0] = 0;
    197     distort->y->coeff[1][0][0][0] = 0;
    198369    distort->y->coeff[0][1][0][0] = 1;
    199     distort->y->mask [1][1][0][0] = 1;
    200370
    201371    return distort;
Note: See TracChangeset for help on using the changeset viewer.