IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 5, 2007, 10:20:06 AM (19 years ago)
Author:
magnier
Message:

various fixes to get FPA/TPA scales correct and consistent

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r10873 r10921  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-01-01 21:05:46 $
     9 *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-01-05 20:20:06 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    236236    wcs->crpix2 = psMetadataLookupF64 (&status, header, "CRPIX2");
    237237    wcs->toSky = psProjectionAlloc (wcs->crval1*PM_RAD_DEG, wcs->crval2*PM_RAD_DEG, PM_RAD_DEG, PM_RAD_DEG, type);
     238    // XXX I think this is wrong for linear proj
    238239
    239240    // test the CDELTi varient
    … …  
    253254        }
    254255
    255         // test the PC00i00j varient:
     256        // FITS WCS PCi,j has units of unity
     257        // wcs->trans has units of degrees/pixel
    256258        wcs->trans->x->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF64 (&status, header, "PC001001"); // == PC1_1
    257259        wcs->trans->x->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF64 (&status, header, "PC001002"); // == PC1_2
    … …  
    283285    // test the CDi_j varient
    284286    if (cdKeys) {
    285         wcs->cdelt1 = 1.0;
    286         wcs->cdelt2 = 1.0;
    287287        wcs->trans->x->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD1_1"); // == PC1_1
    288288        wcs->trans->x->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD1_2"); // == PC1_2
    289289        wcs->trans->y->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD2_1"); // == PC2_1
    290290        wcs->trans->y->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD2_2"); // == PC2_2
     291        wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     292        wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
    291293        return wcs;
    292294    }
    … …  
    321323
    322324    // XXX make it optional to write out CDi_j terms, or other versions
    323     // solve for CDELT1,2 (degrees / pixel)
    324     double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->y->coeff[1][0]);
    325     double cdelt2 = hypot (wcs->trans->x->coeff[0][1], wcs->trans->y->coeff[0][1]);
     325    // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity
     326    double cdelt1 = wcs->cdelt1;
     327    double cdelt2 = wcs->cdelt2;
    326328    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
    327329    psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
    … …  
    336338    // XXX currently, Elixir/DVO cannot accept mixed orders
    337339    // XXX need to respect the masks
     340    // XXX is wcs->cdelt1,2 always consistent?
    338341    int fitOrder = wcs->trans->x->nX;
    339342    if (fitOrder > 1) {
    … …  
    345348                    continue;
    346349                sprintf (name, "PCA1X%1dY%1d", i, j);
    347                 psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
     350                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
    348351                sprintf (name, "PCA2X%1dY%1d", i, j);
    349                 psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j));
     352                psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
    350353            }
    351354        }
    … …  
    357360
    358361// interpret header WCS (only handles traditional WCS for the moment)
    359 // pixelScale is the pixel size in microns per pixel
     362// pixelScale is the pixel size in microns/pixel
    360363bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double pixelScale)
    361364{
    362365    psPlaneTransform *toFPA;
    363366
    364     /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
    365      * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following
    366      * elements defined:
    367      * type (projection type)
    368      * crval1,2 (in RA,DEC degrees)
    369      * crpix1,2
    370      * cdelt1,2 (in degrees / pixel)
    371      * pixelScale (microns / pixel)
    372      *
    373      * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
    374      * 1) the output is in microns (not degrees): divide by cdelt1,2
    375      * 2) X,Y are applied directly, without an applied Xo,Yo offset
    376      * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
    377      */
    378 
    379     // create transformation with 0,0 reference pixel
     367    // create transformation with 0,0 reference pixel and units of degrees/pixel
    380368    toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
    381369
    382     // modify scale of toFPA to yield L,M in microns
    383     double cdelt1 = hypot (toFPA->x->coeff[1][0], toFPA->x->coeff[0][1]);
    384     double cdelt2 = hypot (toFPA->y->coeff[1][0], toFPA->y->coeff[0][1]);
     370    // modify scale of toFPA to have units of microns/pixel
     371    // cdelt1,2 has units of degree/pixel
    385372    for (int i = 0; i <= toFPA->x->nX; i++) {
    386373        for (int j = 0; j <= toFPA->x->nX; j++) {
    387             toFPA->x->coeff[i][j] *= pixelScale/cdelt1;
    388             toFPA->y->coeff[i][j] *= pixelScale/cdelt2;
    389         }
    390     }
    391 
    392     // scale from TPA (microns) to degrees
    393     double pdelt1 = cdelt1 / pixelScale;
    394     double pdelt2 = cdelt2 / pixelScale;
     374            toFPA->x->coeff[i][j] *= pixelScale/wcs->cdelt1;
     375            toFPA->y->coeff[i][j] *= pixelScale/wcs->cdelt2;
     376        }
     377    }
     378
     379    // pdelt1,2 has units of degree/micron
     380    double pdelt1 = wcs->cdelt1 / pixelScale;
     381    double pdelt2 = wcs->cdelt2 / pixelScale;
    395382
    396383    // projection from TPA (linear microns) to SKY (radians)
    … …  
    485472     */
    486473
    487     // create transformation with 0,0 reference pixel
     474    // create transformation with 0,0 reference pixel and units of microns/pixel
    488475    chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
    489476
    … …  
    502489bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs)
    503490{
    504     // projection from TPA to SKY
    505     fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->toSky->Xs, wcs->toSky->Ys, wcs->toSky->type);
     491    // projection from TPA (microns) to SKY (radians)
     492    // cdelt1,2 has units of degrees/micron
     493    fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->cdelt1*PM_RAD_DEG, wcs->cdelt2*PM_RAD_DEG, wcs->toSky->type);
    506494
    507495    // create transformation with 0,0 reference pixel
    508496    fpa->toTPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
     497
     498    // convert fpa->toTPA to units of unity (microns/micron)
     499    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
     500        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
     501            fpa->toTPA->x->coeff[i][j] /= wcs->cdelt1;
     502            fpa->toTPA->y->coeff[i][j] /= wcs->cdelt2;
     503        }
     504    }
    509505
    510506    // XXX this needs to perform the full (non-linear) inversion
    … …  
    525521
    526522    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
     523    // XXX not really: toTPA needs to have unity scale for distortion fitting function
    527524    if (!psPlaneTransformIsDiagonal (fpa->toTPA))
    528525        psAbort ("psastro", "invalid TPA transformation");
    … …  
    535532    pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
    536533
    537     // convert projection from TPA to SKY into wcs projection (degrees to radians)
     534    // convert projection from FPA to SKY into wcs projection (degrees to radians)
    538535    wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, fpa->toSky->type);
    539536    wcs->crval1 = fpa->toSky->R*PM_DEG_RAD;
    … …  
    543540    psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
    544541
    545     // create wcs transform from toFPA, result converts pixels to microns
     542    // create wcs transform from toFPA, resulting transformation has units of microns/pixel
    546543    // adjust wcs transform to use center as reference coordinate
    547544    psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
    … …  
    552549    psFree (center);
    553550
    554     // pdelt1,2 converts from microns->degrees
     551    // pdelt1,2 has units of degrees/micron
    555552    double pdelt1 = fpa->toSky->Xs * PM_DEG_RAD;
    556553    double pdelt2 = fpa->toSky->Ys * PM_DEG_RAD;
    557554
    558     // cdelt1,2 converts from pixels->degrees
    559     wcs->cdelt1 = pdelt1 * hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]);
    560     wcs->cdelt2 = pdelt2 * hypot (chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
    561 
    562     // convert wcs->trans to a matrix which yields L,M in degrees
     555    // convert wcs->trans to a matrix with units of degrees/pixel
    563556    for (int i = 0; i <= wcs->trans->x->nX; i++) {
    564         for (int j = 0; j <= wcs->trans->x->nX; j++) {
     557        for (int j = 0; j <= wcs->trans->x->nY; j++) {
    565558            wcs->trans->x->coeff[i][j] *= pdelt1;
    566559            wcs->trans->y->coeff[i][j] *= pdelt2;
    567560        }
    568561    }
     562
     563    // cdelt1,2 has units of degrees/pixel
     564    wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     565    wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
     566
    569567    return wcs;
    570568}
    … …  
    588586
    589587    // Chip to FPA transformation is a Cartesian 'projection'
     588    // reference pixel for FPA is 0.0, 0.0
    590589    wcs->toSky = psProjectionAlloc (0.0, 0.0, 1.0, 1.0, PS_PROJ_WRP);
    591 
    592     // reference pixel for FPA is 0.0, 0.0
    593590    wcs->crval1 = 0.0;
    594591    wcs->crval2 = 0.0;
    … …  
    598595
    599596    // adjust wcs transform to use center as reference coordinate
     597    // resulting transformation has units of microns/pixel
    600598    psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
    601599
    … …  
    605603    psFree (center);
    606604
    607     // output coordinates are in FPA pixels
    608     wcs->cdelt1 = 1.0;
    609     wcs->cdelt2 = 1.0;
     605    // output coordinates are in microns : CDELT1,2 has units of microns/pixel
     606    wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     607    wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
    610608
    611609    return wcs;
    … …  
    633631
    634632    // adjust wcs transform to use center as reference coordinate
     633    // resulting transformation has units of unity (microns/micron)
    635634    psPlaneTransformSetCenter (wcs->trans, fpa->toTPA, center->x, center->y);
    636635
    … …  
    640639    psFree (center);
    641640
    642     // output coordinates are in FPA pixels
    643     wcs->cdelt1 = 1.0;
    644     wcs->cdelt2 = 1.0;
     641    // pdelt1,2 has units of degrees/micron
     642    double pdelt1 = fpa->toSky->Xs * PM_DEG_RAD;
     643    double pdelt2 = fpa->toSky->Ys * PM_DEG_RAD;
     644
     645    // convert wcs->trans to units of degree/micron
     646    for (int i = 0; i <= wcs->trans->x->nX; i++) {
     647        for (int j = 0; j <= wcs->trans->x->nY; j++) {
     648            wcs->trans->x->coeff[i][j] *= pdelt1;
     649            wcs->trans->y->coeff[i][j] *= pdelt2;
     650        }
     651    }
     652
     653    // cdelt1,2 has units of degrees/micron
     654    wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     655    wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
    645656
    646657    return wcs;
    … …  
    699710*/
    700711
     712/* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
     713 * wcs->trans, which will convert X,Y in pixels to L,M in degrees.  we also have the following
     714 * elements defined:
     715 * type (projection type)
     716 * crval1,2 (in RA,DEC degrees)
     717 * crpix1,2
     718 * cdelt1,2 (in degrees / pixel)
     719 * pixelScale (microns / pixel)
     720 *
     721 * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways:
     722 * 1) the output is in microns (not degrees): divide by cdelt1,2
     723 * 2) X,Y are applied directly, without an applied Xo,Yo offset
     724 * 3) there is an allowed Lo,Mo term ([0][0] coefficients)
     725 */
     726
Note: See TracChangeset for help on using the changeset viewer.