IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 9, 2006, 4:06:47 PM (20 years ago)
Author:
magnier
Message:

moving WCS to modules, adding high-order terms

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r10599 r10603  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-12-09 21:06:44 $
     9*  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-12-10 02:06:47 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    2222
    2323// interpret header WCS (only handles traditional WCS for the moment)
    24 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     24// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
    2525bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic)
    2626{
    2727
    28     psProjectionType type;
     28    # if (0)
     29        psProjectionType type;
    2930    bool status, pcKeys, cdKeys;
    3031    float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
    … …  
    5556
    5657    // what type of WCS keywords are available?
     58    // XXX add check for CROTA2
     59    int fitOrder = psMetadataLookupS32 (&isPoly, header, "NPLYTERM");
    5760    psMetadataLookupF32 (&pcKeys, header, "PC001001");
    5861    psMetadataLookupF32 (&cdKeys, header, "CD1_1");
    5962
    6063    if (cdKeys && pcKeys) {
    61         psLogMsg ("psastro", 2, "warning: both CDi_j and PC00i00j defined in headers, using CDi_j terms\n");
     64        // XXX make this an option
     65        psLogMsg ("psastro", 2, "warning: both CDi_j and PC00i00j defined in headers, using PC00i00j terms\n");
    6266    }
    6367    if (!cdKeys && !pcKeys) {
    … …  
    6670        return false;
    6771    }
     72    if (isPoly && !pcKeys) {
     73        psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms");
     74        return false;
     75        if (fitOrder = 0)
     76            fitOrder = 1;
     77        if ((fitOrder > 3) || (fitOrder < 1)) {
     78            psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder);
     79            return false;
     80        }
     81    } else {
     82        fitOrder = 1;
     83    }
     84
     85    // construct a transformation from X,Y in pixels to L,M in pixels
     86    // NOTE that the WCS keywords convert X,Y to degrees first (using cdelt1,2)
     87    // and then define a transformation from degrees to degrees
     88    psPlaneTransform *wcsTrans = psPlaneTransformAlloc (fitOrder, fitOrder);
    6889
    6990    crval1 = psMetadataLookupF32 (&status, header, "CRVAL1");
    … …  
    7899
    79100        // test the CROTAi varient:
     101        // XXX double check lambda..
    80102        double rotate = psMetadataLookupF32 (&status, header, "CROTA2");
    81103        if (status) {
    82104            double Lambda = cdelt2 / cdelt1;
    83             pc1_1 =  cos(rotate*RAD_DEG);
    84             pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
    85             pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
    86             pc2_2 =  cos(rotate*RAD_DEG);
     105            wcsTrans->x->coeff[1][0] = +cos(rotate*PS_RAD_DEG); // == PC1_1
     106            wcsTrans->x->coeff[0][1] = -sin(rotate*PS_RAD_DEG) * Lambda; // == PC1_2
     107            wcsTrans->y->coeff[1][0] = +sin(rotate*PS_RAD_DEG) / Lambda; // == PC2_1
     108            wcsTrans->y->coeff[1][0] = +cos(rotate*PS_RAD_DEG); // == PC2_2
    87109            goto got_matrix;
    88110        }
    89111
    90112        // test the PC00i00j varient:
    91         pc1_1 = psMetadataLookupF32 (&status, header, "PC001001");
    92         if (status) {
    93             pc1_2 = psMetadataLookupF32 (&status, header, "PC001002");
    94             pc2_1 = psMetadataLookupF32 (&status, header, "PC002001");
    95             pc2_2 = psMetadataLookupF32 (&status, header, "PC002002");
    96 
    97             // XXX EAM : add Elixir polynomial terms here eventually
     113        wcsTrans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1
     114        wcsTrans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2
     115        wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1
     116        wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2
     117
     118        if (isPoly) {
     119            // Elixir-style polynomial terms
     120            for (int i = 0; i < wcsTrans->x->nX; i++) {
     121                for (int j = 0; j < wcsTrans->x->nX; j++) {
     122                    if (i + j < 2)
     123                        continue;
     124                    if (i + j > fitOrder)
     125                        continue;
     126                    sprintf (name, "PCA1dX%1dY%1d", i, j);
     127                    wcsTrans->x->coeff[i][j] = psMetadataLookupF32 (&status, header, name);
     128                }
     129            }
     130            // XXX currently, Elixir/DVO cannot accept mixed orders
     131            for (int i = 0; i < wcsTrans->y->nX; i++) {
     132                for (int j = 0; j < wcsTrans->y->nX; j++) {
     133                    if (i + j < 2)
     134                        continue;
     135                    if (i + j > fitOrder)
     136                        continue;
     137                    sprintf (name, "PCA2dX%1dY%1d", i, j);
     138                    wcsTrans->y->coeff[i][j] = psMetadataLookupF32 (&status, header, name);
     139                }
     140            }
    98141            goto got_matrix;
    99142        }
    … …  
    104147    // test the CDi_j varient
    105148    if (cdKeys) {
    106         pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1");
    107         pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2");
    108         pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1");
    109         pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2");
    110 
    111         // renormalize to cdelt1, cdelt2, etc
    112         double scale = hypot (pc1_1, pc1_2);
     149        wcsTrans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD1_1"); // == PC1_1
     150        wcsTrans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD1_2"); // == PC1_2
     151        wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1
     152        wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2
     153
     154        // normalize rotation matrix, generate cdelt1, cdelt2
     155        double scale = hypot (wcsTrans->x->coeff[1][0], wcsTrans->y->coeff[0][1]);
    113156        cdelt1 = cdelt2 = scale;
    114         pc1_1 /= scale;
    115         pc1_2 /= scale;
    116         pc2_1 /= scale;
    117         pc2_2 /= scale;
     157        wcsTrans->x->coeff[1][0] /= scale;
     158        wcsTrans->x->coeff[0][1] /= scale;
     159        wcsTrans->y->coeff[1][0] /= scale;
     160        wcsTrans->y->coeff[0][1] /= scale;
    118161        goto got_matrix;
    119162    }
    … …  
    140183    *****/
    141184
     185    /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
     186       wcsTrans, which will con
     187       along with the addition
     188       PC00i00j representation:
     189     * type (CTYPE)
     190     * crval1,2 (in RA,DEC degrees)
     191     * crpix1,2
     192     * cdelt1,2 (in degrees / pixel)
     193     * pci,j (normalized)
     194     *
     195     * we also have plateScale (radians / physical TPA units)
     196     * now we convert to pmFPA terms
     197     */
     198
     199    /*** XXXX need to extend these formulae to higher-order terms ***/
     200
    142201    {
    143         double rX = 1.0;
    144         double rY = 1.0;
    145 
    146202        // XXX free an existing toFPA?
    147203        psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1);
    148204
    149         // basic transformation from chip to FPA
     205        //
     206
     207        // basic transformation from chip to FPA (FPA in pixels)
    150208        toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
    151209        toFPA->x->coeff[1][0] = pc1_1;
    … …  
    158216        toFPA->y->mask[1][1]  = 1;
    159217
     218        // scale from FPA to TPA (microns / pixel)
     219        double pdelt1 = cdelt1*PS_RAD_DEG / plateScale;
     220        double pdelt2 = cdelt2*PS_RAD_DEG / plateScale;
     221        float rX = 1.0;
     222        float rY = 1.0;
     223
    160224        // projection from TPA to SKY
    161         psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
     225        psProjection *toSky = psProjectionAlloc (crval1*PS_RAD_DEG, crval2*PS_RAD_DEG, plateScale, plateScale, type);
    162226
    163227        if (fpa->toSky == NULL)
    164228        {
    165             // XXX for now, use the identity for TPA <--> FPA
    166             fpa->toTangentPlane = psPlaneDistortIdentity (1);
    167             fpa->fromTangentPlane = psPlaneDistortIdentity (1);
     229            fpa->toTPA = psPlaneDistortIdentity (1);
     230            fpa->fromTPA = psPlaneDistortIdentity (1);
     231            fpa->toTPA->x->coeff[1][0][0][0] = pdelt1;
     232            fpa->toTPA->y->coeff[0][1][0][0] = pdelt2;
     233            fpa->fromTPA->x->coeff[1][0][0][0] = 1.0 / pdelt1;
     234            fpa->fromTPA->y->coeff[0][1][0][0] = 1.0 / pdelt2;
    168235            fpa->toSky = toSky;
    169236        } else
    170237        {
    171             if (fpa->toTangentPlane == NULL)
     238            if (fpa->toTPA == NULL)
    172239                psAbort ("wcs", "projection defined, tangent-plane not defined");
    173             if (fpa->fromTangentPlane == NULL)
     240            if (fpa->fromTPA == NULL)
    174241                psAbort ("wcs", "projection defined, tangent-plane not defined");
    175242
    176             // adjust for common toSky for mosaic:
    177             // ignore the TPA since toTPA is identity
     243            // convert from pixels on this chip to pixels on reference chip
     244            // rX has units of refpixels / pixel
     245            rX = pdelt1 / fpa->toTPA->x->coeff[1][0][0][0];
     246            rY = pdelt2 / fpa->toTPA->y->coeff[0][1][0][0];
     247
     248            // correct to common plate scale (output is in refpixel units)
     249            toFPA->x->coeff[0][0] *= rX;
     250            toFPA->x->coeff[1][0] *= rX;
     251            toFPA->x->coeff[0][1] *= rX;
     252            toFPA->y->coeff[0][0] *= rY;
     253            toFPA->y->coeff[1][0] *= rY;
     254            toFPA->y->coeff[0][1] *= rY;
     255
     256            // adjust for common toSky, toTPA for mosaic:
    178257            // find the FPA coordinate of 0,0 for this chip.
     258            psPlane *chip = psPlaneAlloc();
    179259            psPlane *fp = psPlaneAlloc();
    180             psPlane *chip = psPlaneAlloc();
     260            psPlane *tp = psPlaneAlloc();
    181261            psSphere *sky = psSphereAlloc();
    182262            chip->x = chip->y = 0;
    183263
    184264            psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
    185             p_psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate
    186             p_psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     265            psPlaneDistortApply (tp, fpa->toTPA, fp, 0.0, 0.0);
     266            p_psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate
     267            p_psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     268            psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0);
    187269
    188270            toFPA->x->coeff[0][0] = fp->x;
    189271            toFPA->y->coeff[0][0] = fp->y;
    190 
    191             rX = toSky->Xs/fpa->toSky->Xs;
    192             rY = toSky->Ys/fpa->toSky->Ys;
    193 
    194             // correct to common plate scale
    195             toFPA->x->coeff[1][0] *= rX;
    196             toFPA->x->coeff[0][1] *= rX;
    197             toFPA->y->coeff[1][0] *= rY;
    198             toFPA->y->coeff[0][1] *= rY;
    199272
    200273            psFree (fp);
    … …  
    207280        chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA);
    208281
    209         while (fpa->toSky->R <        0)
     282        // this can take a very long time...
     283        while (fpa->toSky->R < 0)
    210284            fpa->toSky->R += 2.0*M_PI;
    211285        while (fpa->toSky->R > 2.0*M_PI)
    … …  
    213287
    214288        // remove the correction to the common plate scale
     289        // NOTE: this assumes 1) we are reading in headers generated using per-chip astrometry
     290        // and 2) we are going to measure the mosaic distortion in the next step.
     291        // XXX perhaps make this its own function? (I'll need to store rX somewhere).
    215292        if (isMosaic)
    216293        {
     294            chip->toFPA->x->coeff[0][0] /= rX;
    217295            chip->toFPA->x->coeff[1][0] /= rX;
    218296            chip->toFPA->x->coeff[0][1] /= rX;
     297            chip->toFPA->y->coeff[0][0] /= rY;
    219298            chip->toFPA->y->coeff[1][0] /= rY;
    220299            chip->toFPA->y->coeff[0][1] /= rY;
    … …  
    232311
    233312        psLogMsg ("psastro", 3, "field center: %f,%f, plate scale: %f,%f (arcsec/pixel)\n",
    234                   DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D,
    235                   3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys);
    236     }
     313                  PS_DEG_RAD*fpa->toSky->R, PS_DEG_RAD*fpa->toSky->D,
     314                  3600*PS_DEG_RAD*fpa->toSky->Xs, 3600*PS_DEG_RAD*fpa->toSky->Ys);
     315    }
     316    # endif
    237317    return true;
    238318}
    … …  
    293373
    294374        // solve for CDELT1,2 (degrees / pixel)
    295         cdelt1 = DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
    296     cdelt2 = DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0];
     375        cdelt1 = PS_DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
     376    cdelt2 = PS_DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0];
    297377
    298378    // L,M = toFPA(X,Y)
    … …  
    365445        coords[0].crval1 -= 360.0;
    366446
    367     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*DEG_RAD);
    368     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*DEG_RAD);
     447    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
     448    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
    369449    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", Xo);
    370450    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", Yo);
    … …  
    397477    # else
    398478
    399         psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*DEG_RAD);
    400     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*DEG_RAD);
     479        psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
     480    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
    401481    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
    402482    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
    403     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
    404     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
     483    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale);
     484    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale);
    405485
    406486    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
    … …  
    486566
    487567    // XXX need to handle the plateScale??
    488     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*DEG_RAD);
    489     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*DEG_RAD);
     568    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",  PS_META_REPLACE, "", toSky->R*PS_DEG_RAD);
     569    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",  PS_META_REPLACE, "", toSky->D*PS_DEG_RAD);
    490570    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",  PS_META_REPLACE, "", 0.0);
    491571    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",  PS_META_REPLACE, "", 0.0);
    492     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
    493     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
     572    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",  PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale);
     573    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",  PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale);
    494574
    495575    if (toTP->x->nX != toTP->x->nY)
Note: See TracChangeset for help on using the changeset viewer.