Index: /trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10611) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10612) @@ -7,6 +7,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-10 04:14:45 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-10 18:30:07 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -21,12 +21,92 @@ #include "pmAstrometryWCS.h" -// interpret header WCS (only handles traditional WCS for the moment) -// plateScale is nominal physical scale on tangent plane (radians / TPA physical units) -bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) +// the following functions support coordinate transformations direcly related to the FITS WCS +// keywords. The FITS WCS allows for only a single level of transformation, thus it is not +// appropriate for mosaic astrometry consisting of telescope distortion plus chip terms. +// Below, we support the Elixir convention of using two connected FITS headers to define two +// levels of coordinate transformation. In the pmFPA structure, the projection, distortion, +// and FPA-to-Chip transformations are carried independently. NOTE: The FITS WCS keywords do +// not represent a simply polynomial. Instead, they have no constant term, and the coordinates +// are corrected to a reference pixel before the polynomial transformation is applied. + +static void pmAstromWCSFree (pmAstromWCS *wcs) +{ + + if (!wcs) + return; + psFree (wcs->trans); +} + +pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder) +{ + + pmAstromWCS *wcs = (pmAstromWCS *) psAlloc(sizeof(pmAstromWCS)); + psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree); + + wcs->trans = psPlaneTransformAlloc (nXorder, nYorder); + + memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE); + memset (wcs->ctype2, 0, PM_ASTROM_WCS_TYPE_SIZE); + return wcs; +} + +bool pmAstromWCStoSky (psSphere *sky, pmAstromWCS *wcs, psPlane *chip) +{ + + if (chip == NULL) + return false; + if (sky == NULL) + return false; + if (wcs == NULL) + return false; + + psPlane *Chip = psPlaneAlloc(); + psPlane *FP = psPlaneAlloc(); + + Chip->x = chip->x - wcs->crpix1; + Chip->y = chip->y - wcs->crpix2; + + psPlaneTransformApply (FP, wcs->trans, Chip); + psDeproject (sky, FP, wcs->toSky); // find the RA,DEC coord of the focal-plane coordinate + + psFree (Chip); + psFree (FP); + return true; +} + +bool pmAstromWCStoChip (psPlane *chip, pmAstromWCS *wcs, psSphere *sky) +{ + + if (chip == NULL) + return false; + if (sky == NULL) + return false; + if (wcs == NULL) + return false; + + psError(PS_ERR_UNKNOWN, true, "not yet implemented: needs to invert the transformation"); + return false; + + psPlane *Chip = psPlaneAlloc(); + psPlane *FP = psPlaneAlloc(); + + psProject (FP, sky, wcs->toSky); // find the RA,DEC coord of the focal-plane coordinate + + // I need the inverse of wcs->transform at this point + psPlaneTransformApply (Chip, wcs->trans, FP); + + chip->x = Chip->x + wcs->crpix1; + chip->y = Chip->y + wcs->crpix2; + + psFree (Chip); + psFree (FP); + return true; +} + +// interpret header WCS keywords (only handles traditional WCS for the moment) +pmAstromWCS *pmAstromWCSfromHeader (psMetadata *header) { psProjectionType type; - psPlaneTransform *toFPA; bool status, pcKeys, cdKeys, isPoly; - float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2; char name[16]; // used to store FITS keyword below (always < 8, so 16 should be safe!) @@ -35,5 +115,5 @@ if (!status) { psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n"); - return false; + return NULL; } @@ -51,5 +131,5 @@ if (type == PS_PROJ_NTYPE) { psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype); - return false; + return NULL; } @@ -67,33 +147,35 @@ psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms"); // XXX we could default here to RA, DEC, ROTANGLE - return false; + return NULL; } if (isPoly && !pcKeys) { psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms"); - return false; + return NULL; if (fitOrder == 0) fitOrder = 1; if ((fitOrder > 3) || (fitOrder < 1)) { psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder); - return false; + return NULL; } } else { fitOrder = 1; } + + pmAstromWCS *wcs = pmAstromWCSAlloc (fitOrder, fitOrder); // construct a transformation from X,Y in pixels to L,M in pixels // NOTE that the WCS keywords convert X,Y to degrees first (using cdelt1,2) // and then define a transformation from degrees to degrees - psPlaneTransform *wcsTrans = psPlaneTransformAlloc (fitOrder, fitOrder); - - crval1 = psMetadataLookupF32 (&status, header, "CRVAL1"); - crval2 = psMetadataLookupF32 (&status, header, "CRVAL2"); - crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1"); - crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2"); + + wcs->crval1 = psMetadataLookupF32 (&status, header, "CRVAL1"); + wcs->crval2 = psMetadataLookupF32 (&status, header, "CRVAL2"); + wcs->crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1"); + wcs->crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2"); + wcs->toSky = psProjectionAlloc (wcs->crval1*PM_RAD_DEG, wcs->crval2*PM_RAD_DEG, PM_RAD_DEG, PM_RAD_DEG, type); // test the CDELTi varient if (pcKeys) { - cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1"); - cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2"); + wcs->cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1"); + wcs->cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2"); // test the CROTAi varient: @@ -101,21 +183,22 @@ double rotate = psMetadataLookupF32 (&status, header, "CROTA2"); if (status) { - wcsTrans->x->coeff[1][0] = +cdelt1 * cos(rotate*PS_RAD_DEG); // == PC1_1 - wcsTrans->x->coeff[0][1] = -cdelt2 * sin(rotate*PS_RAD_DEG); // == PC1_2 - wcsTrans->y->coeff[1][0] = +cdelt1 * sin(rotate*PS_RAD_DEG); // == PC2_1 - wcsTrans->y->coeff[1][0] = +cdelt2 * cos(rotate*PS_RAD_DEG); // == PC2_2 - goto got_matrix; + wcs->trans->x->coeff[1][0] = +wcs->cdelt1 * cos(rotate*PM_RAD_DEG); // == PC1_1 + wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2 + wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1 + wcs->trans->y->coeff[1][0] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2 + return wcs; } // test the PC00i00j varient: - wcsTrans->x->coeff[1][0] = cdelt1 * psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1 - wcsTrans->x->coeff[0][1] = cdelt2 * psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2 - wcsTrans->y->coeff[1][0] = cdelt1 * psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1 - wcsTrans->y->coeff[0][1] = cdelt2 * psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2 + wcs->trans->x->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1 + wcs->trans->x->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2 + wcs->trans->y->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1 + wcs->trans->y->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2 if (isPoly) { // Elixir-style polynomial terms - for (int i = 0; i < wcsTrans->x->nX; i++) { - for (int j = 0; j < wcsTrans->x->nX; j++) { + // XXX currently, Elixir/DVO cannot accept mixed orders + for (int i = 0; i < fitOrder; i++) { + for (int j = 0; j < fitOrder; j++) { if (i + j < 2) continue; @@ -123,37 +206,40 @@ continue; sprintf (name, "PCA1dX%1dY%1d", i, j); - wcsTrans->x->coeff[i][j] = pow(cdelt1, i) * pow(cdelt2, j) * psMetadataLookupF32 (&status, header, name); + wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name); + sprintf (name, "PCA2dX%1dY%1d", i, j); + wcs->trans->y->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name); } } - // XXX currently, Elixir/DVO cannot accept mixed orders - for (int i = 0; i < wcsTrans->y->nX; i++) { - for (int j = 0; j < wcsTrans->y->nX; j++) { - if (i + j < 2) - continue; - if (i + j > fitOrder) - continue; - sprintf (name, "PCA2dX%1dY%1d", i, j); - wcsTrans->y->coeff[i][j] = pow(cdelt1, i) * pow(cdelt2, j) * psMetadataLookupF32 (&status, header, name); - } - } - } - goto got_matrix; + } + return wcs; } // test the CDi_j varient if (cdKeys) { - wcsTrans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD1_1"); // == PC1_1 - wcsTrans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD1_2"); // == PC1_2 - wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1 - wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2 - goto got_matrix; + wcs->cdelt1 = 1.0; + wcs->cdelt2 = 1.0; + wcs->trans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD1_1"); // == PC1_1 + wcs->trans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD1_2"); // == PC1_2 + wcs->trans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1 + wcs->trans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2 + return wcs; } psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); - return false; - -got_matrix: + return NULL; +} + +// interpret header WCS (only handles traditional WCS for the moment) +// plateScale is nominal physical scale on tangent plane (radians / TPA physical units) +bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) +{ + psPlaneTransform *toFPA; + + pmAstromWCS *wcs = pmAstromWCSfromHeader (header); + if (!wcs) { + return false; + } /* at this point, we have extracted from the header the WCS terms in the form of a polynomial, - * wcsTrans, which will convert X,Y in pixels to L,M in degrees. we also have the following + * wcs->trans, which will convert X,Y in pixels to L,M in degrees. we also have the following * elements defined: * type (CTYPE) @@ -163,5 +249,5 @@ * plateScale (radians / physical TPA units) * - * now we convert wcsTrans to toFPA, which is different from wcsTrans in 3 important ways: + * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways: * 1) the output is in pixel (not degrees): divide by cdelt1,2 raised to an appropriate power * 2) X,Y are applied directly, without an applied Xo,Yo offset @@ -169,17 +255,28 @@ */ - /*** XXXX need to extend these formulae to higher-order terms ***/ - - // XXX free an existing toFPA? + // convert wcs->trans to a matrix which yields L,M in pixels + double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]); + double cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]); + for (int i = 0; i <= wcs->trans->x->nX; i++) { + for (int j = 0; j <= wcs->trans->x->nX; j++) { + wcs->trans->x->coeff[i][j] /= cdelt1; + wcs->trans->y->coeff[i][j] /= cdelt2; + } + } + + // validate fit order + int fitOrder = wcs->trans->x->nX; toFPA = psPlaneTransformAlloc(fitOrder, fitOrder); /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j - * we can transform L(x,y) into L'(x+xo,y+yo) by taking the derivatives of both sides and + * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and * noting that the constant term in each is the coefficient in the case of L(x,y) and is the - * value of L'(xo,yo) in the second case. in this case, xo,yo = crpix1,2 + * value of L'(-xo,-yo) in the second case. in this case, xo,yo = crpix1,2 */ - psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcsTrans->x); - psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcsTrans->y); + psPolynomial2D *tmp; + + psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcs->trans->x); + psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcs->trans->y); for (int i = 0; i <= fitOrder; i++) { @@ -187,44 +284,34 @@ psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx); for (int j = 0; j <= fitOrder; j++) { - toFPA->x->mask[i][j] = wcsTrans->x->mask[i][j]; - toFPA->y->mask[i][j] = wcsTrans->y->mask[i][j]; - toFPA->x->coeff[i][j] = (toFPA->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, crpix1, crpix2) / tgamma(i+1) / tgamma(j+1) / cdelt1; - toFPA->y->coeff[i][j] = (toFPA->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, crpix1, crpix2) / tgamma(i+1) / tgamma(j+1) / cdelt2; - - // take the next derivative wrt y - psPolynomial2D_dY(xPy, xPy); - psPolynomial2D_dY(yPy, yPy); - } - psFree (xPy); - psFree (yPy); - // take the next derivative wrt x - psPolynomial2D_dX(xPx, xPx); - psPolynomial2D_dX(yPx, yPx); - } - psFree (xPx); - psFree (yPx); - - // save until we verify the transformation - # if (0) - // basic transformation from chip to FPA (FPA in pixels) - toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2); - toFPA->x->coeff[1][0] = pc1_1; - toFPA->x->coeff[0][1] = pc1_2; - toFPA->x->mask[1][1] = 1; - - toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2); - toFPA->y->coeff[1][0] = pc2_1; - toFPA->y->coeff[0][1] = pc2_2; - toFPA->y->mask[1][1] = 1; - # endif + toFPA->x->mask[i][j] = wcs->trans->x->mask[i][j]; + toFPA->y->mask[i][j] = wcs->trans->y->mask[i][j]; + toFPA->x->coeff[i][j] = (toFPA->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1); + toFPA->y->coeff[i][j] = (toFPA->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, -wcs->crpix1, -wcs->crpix2) / tgamma(i+1) / tgamma(j+1); + + // take the next derivative wrt y, catch output (is NULL on last pass) + tmp = psPolynomial2D_dY(NULL, xPy); + psFree (xPy); + xPy = tmp; + tmp = psPolynomial2D_dY(NULL, yPy); + psFree (yPy); + yPy = tmp; + } + // take the next derivative wrt x, catch output (is NULL on last pass) + tmp = psPolynomial2D_dX(NULL, xPx); + psFree (xPx); + xPx = tmp; + tmp = psPolynomial2D_dX(NULL, yPx); + psFree (yPx); + yPx = tmp; + } // scale from FPA to TPA (microns / pixel) - double pdelt1 = cdelt1*PS_RAD_DEG / plateScale; - double pdelt2 = cdelt2*PS_RAD_DEG / plateScale; + double pdelt1 = cdelt1*PM_RAD_DEG / plateScale; + double pdelt2 = cdelt2*PM_RAD_DEG / plateScale; float rX = 1.0; float rY = 1.0; // projection from TPA to SKY - psProjection *toSky = psProjectionAlloc (crval1*PS_RAD_DEG, crval2*PS_RAD_DEG, plateScale, plateScale, type); + psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, plateScale, plateScale, wcs->toSky->type); if (fpa->toSky == NULL) { @@ -263,6 +350,6 @@ psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate psPlaneDistortApply (tp, fpa->toTPA, fp, 0.0, 0.0); - p_psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate - p_psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection + psDeproject (sky, tp, toSky); // find the RA,DEC coord of the focal-plane coordinate + psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0); @@ -309,9 +396,5 @@ chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]); - psLogMsg ("psastro", 3, "field center: %f,%f, plate scale: %f,%f (arcsec/pixel)\n", - PS_DEG_RAD*fpa->toSky->R, PS_DEG_RAD*fpa->toSky->D, - 3600*PS_DEG_RAD*fpa->toSky->Xs, 3600*PS_DEG_RAD*fpa->toSky->Ys); - - psFree (wcsTrans); + psFree (wcs->trans); return true; @@ -373,6 +456,6 @@ // solve for CDELT1,2 (degrees / pixel) - cdelt1 = PS_DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0]; - cdelt2 = PS_DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0]; + cdelt1 = PM_DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0]; + cdelt2 = PM_DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0]; // L,M = toFPA(X,Y) @@ -445,6 +528,6 @@ coords[0].crval1 -= 360.0; - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*PS_DEG_RAD); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*PS_DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*PM_DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*PM_DEG_RAD); psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", Xo); psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", Yo); @@ -477,10 +560,10 @@ # else - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*PS_DEG_RAD); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*PS_DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*PM_DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*PM_DEG_RAD); psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", toFPA->x->coeff[0][0]); psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", toFPA->y->coeff[0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*PM_DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*PM_DEG_RAD*plateScale); psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale); @@ -566,10 +649,10 @@ // XXX need to handle the plateScale?? - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*PS_DEG_RAD); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*PS_DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*PM_DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*PM_DEG_RAD); psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", 0.0); psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", 0.0); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*PS_DEG_RAD*plateScale); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*PS_DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*PM_DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*PM_DEG_RAD*plateScale); if (toTP->x->nX != toTP->x->nY) @@ -738,2 +821,17 @@ *****/ + +// save until we verify the transformation +# if (0) + // basic transformation from chip to FPA (FPA in pixels) + toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2); +toFPA->x->coeff[1][0] = pc1_1; +toFPA->x->coeff[0][1] = pc1_2; +toFPA->x->mask[1][1] = 1; + +toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2); +toFPA->y->coeff[1][0] = pc2_1; +toFPA->y->coeff[0][1] = pc2_2; +toFPA->y->mask[1][1] = 1; +# endif + Index: /trunk/psModules/src/astrom/pmAstrometryWCS.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10611) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10612) @@ -7,6 +7,6 @@ * @author EAM, IfA * -* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-10 02:06:47 $ +* @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-10 18:30:07 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -15,4 +15,22 @@ #ifndef PM_ASTROMETRY_WCS_H #define PM_ASTROMETRY_WCS_H + +#define PM_ASTROM_WCS_TYPE_SIZE 80 +typedef struct +{ + char ctype1[PM_ASTROM_WCS_TYPE_SIZE]; + char ctype2[PM_ASTROM_WCS_TYPE_SIZE]; + double crval1, crval2; + double crpix1, crpix2; + double cdelt1, cdelt2; + psProjection *toSky; + psPlaneTransform *trans; +} +pmAstromWCS; + +pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder); +bool pmAstromWCStoSky (psSphere *sky, pmAstromWCS *wcs, psPlane *chip); +bool pmAstromWCStoChip (psPlane *chip, pmAstromWCS *wcs, psSphere *sky); +pmAstromWCS *pmAstromWCSfromHeader (psMetadata *header); bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic); @@ -25,6 +43,11 @@ bool psPlaneDistortIsIdentity (psPlaneDistort *distort); -# define PS_DEG_RAD 57.295779513082322 -# define PS_RAD_DEG 0.017453292519943 +# define PM_DEG_RAD 57.295779513082322 +# define PM_RAD_DEG 0.017453292519943 #endif // PM_ASTROMETRY_WCS_H + +/* + * the wcs->trans component defines a polynomial which converts (x-crpix1),(y-crpix2) to + * L,M in degrees + */ Index: /trunk/psModules/test/astrom/tap_pmAstrometryWCS.c =================================================================== --- /trunk/psModules/test/astrom/tap_pmAstrometryWCS.c (revision 10611) +++ /trunk/psModules/test/astrom/tap_pmAstrometryWCS.c (revision 10612) @@ -83,4 +83,84 @@ } + { + diag("test pmAstromReadWCS"); + psMemId id = psMemGetId(); + + psPlaneTransform *wcsTrans = psPlaneTransformAlloc (1, 1); + wcsTrans->x->coeff[0][0] = 0.0; + wcsTrans->y->coeff[0][0] = 0.0; + wcsTrans->x->coeff[1][0] = 1.0; + wcsTrans->x->coeff[0][1] = 0.0; + wcsTrans->y->coeff[1][0] = 0.0; + wcsTrans->y->coeff[0][1] = 1.0; + + + + // construct a header with a simple set of WCS values + // convert to pmFPA components and check + + psMetadata *header = psMetadataAlloc(); + + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--TAN"); + + // center coords (R,D) + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", 0.0); + + // center coords (X,Y) + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", 10.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", 10.0); + + // degrees per pixel + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", 1.0/3600.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", 1.0/3600.0); + + // rotation matrix + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", 1.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", 1.0); + + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (NULL, "test"); + + // toFPA carries pixel scale (pixels per micron) + // toTPA carries plate scale (microns per arcsecond) + bool status = pmAstromReadWCS (fpa, chip, header, 25.0*PS_RAD_DEG/3600.0, false); + ok (status, "converted WCS keywords to FPA astrometry"); + skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + ok(fpa->toSky->type == PS_PROJ_TAN, "correct projection (TAN)"); + + // make these tests double + ok_float(fpa->toSky->R*PS_DEG_RAD, 0.0, "projection center RA %f", fpa->toSky->R*PS_DEG_RAD); + ok_float(fpa->toSky->R*PS_DEG_RAD, 0.0, "projection center DEC %f", fpa->toSky->R*PS_DEG_RAD); + + ok_float(fpa->toSky->Xs, 25.0*PS_RAD_DEG/3600.0, "projection X scale %f", fpa->toSky->Xs); + ok_float(fpa->toSky->Ys, 25.0*PS_RAD_DEG/3600.0, "projection X scale %f", fpa->toSky->Ys); + + ok_float(fpa->toTPA->x->coeff[1][0][0][0], 0.04, "TP scale (mm per pixel): %f", fpa->toTPA->x->coeff[1][0][0][0]); + ok_float(fpa->toTPA->y->coeff[0][1][0][0], 0.04, "TP scale (mm per pixel): %f", fpa->toTPA->x->coeff[1][0][0][0]); + + ok_float(chip->toFPA->x->coeff[0][0], -10.0, "ref pixel X: %f", chip->toFPA->x->coeff[0][0]); + ok_float(chip->toFPA->y->coeff[0][0], -10.0, "ref pixel Y: %f", chip->toFPA->y->coeff[0][0]); + + ok_float(chip->toFPA->x->coeff[1][0], 1.0, "CD1_1: %f", chip->toFPA->x->coeff[1][0]); + ok_float(chip->toFPA->x->coeff[0][1], 0.0, "CD1_2: %f", chip->toFPA->x->coeff[0][1]); + ok_float(chip->toFPA->y->coeff[1][0], 0.0, "CD2_1: %f", chip->toFPA->y->coeff[1][0]); + ok_float(chip->toFPA->y->coeff[0][1], 1.0, "CD2_2: %f", chip->toFPA->y->coeff[0][1]); + + // apply both systems to real data: + + + psFree (fpa); + psFree (chip); + psFree (header); + + skip_end(); + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); + } + return exit_status(); }