Index: trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10603) +++ trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10606) @@ -1,15 +1,15 @@ /** @file pmAstrometryWCS.c -* -* @brief functions to convert FITS WCS keywords to / from pmFPA structures -* -* @ingroup Astrometry -* -* @author EAM, IfA -* -* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-10 02:06:47 $ -* -* Copyright 2006 Institute for Astronomy, University of Hawaii -*/ + * + * @brief functions to convert FITS WCS keywords to / from pmFPA structures + * + * @ingroup Astrometry + * + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-10 04:14:45 $ + * + * Copyright 2006 Institute for Astronomy, University of Hawaii + */ #ifdef HAVE_CONFIG_H @@ -25,10 +25,9 @@ bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) { - - # if (0) - psProjectionType type; - bool status, pcKeys, cdKeys; + psProjectionType type; + psPlaneTransform *toFPA; + bool status, pcKeys, cdKeys, isPoly; float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2; - float pc1_1, pc1_2, pc2_1, pc2_2; + char name[16]; // used to store FITS keyword below (always < 8, so 16 should be safe!) // interpret header data, convert to crval(i), etc @@ -73,5 +72,5 @@ psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms"); return false; - if (fitOrder = 0) + if (fitOrder == 0) fitOrder = 1; if ((fitOrder > 3) || (fitOrder < 1)) { @@ -102,17 +101,16 @@ double rotate = psMetadataLookupF32 (&status, header, "CROTA2"); if (status) { - double Lambda = cdelt2 / cdelt1; - wcsTrans->x->coeff[1][0] = +cos(rotate*PS_RAD_DEG); // == PC1_1 - wcsTrans->x->coeff[0][1] = -sin(rotate*PS_RAD_DEG) * Lambda; // == PC1_2 - wcsTrans->y->coeff[1][0] = +sin(rotate*PS_RAD_DEG) / Lambda; // == PC2_1 - wcsTrans->y->coeff[1][0] = +cos(rotate*PS_RAD_DEG); // == PC2_2 + 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; } // test the PC00i00j varient: - wcsTrans->x->coeff[1][0] = psMetadataLookupF32 (&status, header, "PC001001"); // == PC1_1 - wcsTrans->x->coeff[0][1] = psMetadataLookupF32 (&status, header, "PC001002"); // == PC1_2 - wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "PC002001"); // == PC2_1 - wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "PC002002"); // == PC2_2 + 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 if (isPoly) { @@ -125,5 +123,5 @@ continue; sprintf (name, "PCA1dX%1dY%1d", i, j); - wcsTrans->x->coeff[i][j] = psMetadataLookupF32 (&status, header, name); + wcsTrans->x->coeff[i][j] = pow(cdelt1, i) * pow(cdelt2, j) * psMetadataLookupF32 (&status, header, name); } } @@ -136,11 +134,9 @@ continue; sprintf (name, "PCA2dX%1dY%1d", i, j); - wcsTrans->y->coeff[i][j] = psMetadataLookupF32 (&status, header, name); + wcsTrans->y->coeff[i][j] = pow(cdelt1, i) * pow(cdelt2, j) * psMetadataLookupF32 (&status, header, name); } } - goto got_matrix; - } - psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); - return false; + } + goto got_matrix; } @@ -151,12 +147,4 @@ wcsTrans->y->coeff[1][0] = psMetadataLookupF32 (&status, header, "CD2_1"); // == PC2_1 wcsTrans->y->coeff[0][1] = psMetadataLookupF32 (&status, header, "CD2_2"); // == PC2_2 - - // normalize rotation matrix, generate cdelt1, cdelt2 - double scale = hypot (wcsTrans->x->coeff[1][0], wcsTrans->y->coeff[0][1]); - cdelt1 = cdelt2 = scale; - wcsTrans->x->coeff[1][0] /= scale; - wcsTrans->x->coeff[0][1] /= scale; - wcsTrans->y->coeff[1][0] /= scale; - wcsTrans->y->coeff[0][1] /= scale; goto got_matrix; } @@ -166,153 +154,165 @@ got_matrix: - /***** - - For mosaic astrometry, we need to have a starting set of projection terms in which the - chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or - microns). This set of projections, coupled with an identity toTPA (ie, no distortion) will - result in substantial errors between the observed and predicted star positions on the focal - plane: this is the measurement of the optical distortion in the camera. At the same time, - we need to carry around the transformations which allow us to make an accurate calculation - of the position of the stars based on the input (per-chip) astrometry. These - transformations will allow us to match the raw and ref stars robustly. To convert the - per-chip astrometry (which may have been calculated with a different plate scale for each - chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust - the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for - each chip (eg, cdelt1). Thus, we need to carry around both the - - *****/ - /* at this point, we have extracted from the header the WCS terms in the form of a polynomial, - wcsTrans, which will con - along with the addition - PC00i00j representation: + * wcsTrans, which will convert X,Y in pixels to L,M in degrees. we also have the following + * elements defined: * type (CTYPE) * crval1,2 (in RA,DEC degrees) * crpix1,2 * cdelt1,2 (in degrees / pixel) - * pci,j (normalized) + * plateScale (radians / physical TPA units) * - * we also have plateScale (radians / physical TPA units) - * now we convert to pmFPA terms + * now we convert wcsTrans to toFPA, which is different from wcsTrans 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 + * 3) there is an allowed Lo,Mo term ([0][0] coefficients) */ /*** XXXX need to extend these formulae to higher-order terms ***/ - { - // XXX free an existing toFPA? - psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1); - - // - + // XXX free an existing toFPA? + 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 + * 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 + */ + + psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcsTrans->x); + psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcsTrans->y); + + for (int i = 0; i <= fitOrder; i++) { + psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx); + 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; - - // scale from FPA to TPA (microns / pixel) - double pdelt1 = cdelt1*PS_RAD_DEG / plateScale; - double pdelt2 = cdelt2*PS_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); - - if (fpa->toSky == NULL) - { - fpa->toTPA = psPlaneDistortIdentity (1); - fpa->fromTPA = psPlaneDistortIdentity (1); - fpa->toTPA->x->coeff[1][0][0][0] = pdelt1; - fpa->toTPA->y->coeff[0][1][0][0] = pdelt2; - fpa->fromTPA->x->coeff[1][0][0][0] = 1.0 / pdelt1; - fpa->fromTPA->y->coeff[0][1][0][0] = 1.0 / pdelt2; - fpa->toSky = toSky; - } else - { - if (fpa->toTPA == NULL) - psAbort ("wcs", "projection defined, tangent-plane not defined"); - if (fpa->fromTPA == NULL) - psAbort ("wcs", "projection defined, tangent-plane not defined"); - - // convert from pixels on this chip to pixels on reference chip - // rX has units of refpixels / pixel - rX = pdelt1 / fpa->toTPA->x->coeff[1][0][0][0]; - rY = pdelt2 / fpa->toTPA->y->coeff[0][1][0][0]; - - // correct to common plate scale (output is in refpixel units) - toFPA->x->coeff[0][0] *= rX; - toFPA->x->coeff[1][0] *= rX; - toFPA->x->coeff[0][1] *= rX; - toFPA->y->coeff[0][0] *= rY; - toFPA->y->coeff[1][0] *= rY; - toFPA->y->coeff[0][1] *= rY; - - // adjust for common toSky, toTPA for mosaic: - // find the FPA coordinate of 0,0 for this chip. - psPlane *chip = psPlaneAlloc(); - psPlane *fp = psPlaneAlloc(); - psPlane *tp = psPlaneAlloc(); - psSphere *sky = psSphereAlloc(); - chip->x = chip->y = 0; - - 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 - psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0); - - toFPA->x->coeff[0][0] = fp->x; - toFPA->y->coeff[0][0] = fp->y; - - psFree (fp); - psFree (sky); - psFree (chip); - psFree (toSky); - } - - chip->toFPA = toFPA; - chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA); - - // this can take a very long time... - while (fpa->toSky->R < 0) - fpa->toSky->R += 2.0*M_PI; - while (fpa->toSky->R > 2.0*M_PI) - fpa->toSky->R -= 2.0*M_PI; - - // remove the correction to the common plate scale - // NOTE: this assumes 1) we are reading in headers generated using per-chip astrometry - // and 2) we are going to measure the mosaic distortion in the next step. - // XXX perhaps make this its own function? (I'll need to store rX somewhere). - if (isMosaic) - { - chip->toFPA->x->coeff[0][0] /= rX; - chip->toFPA->x->coeff[1][0] /= rX; - chip->toFPA->x->coeff[0][1] /= rX; - chip->toFPA->y->coeff[0][0] /= rY; - chip->toFPA->y->coeff[1][0] /= rY; - chip->toFPA->y->coeff[0][1] /= rY; - } - - psTrace ("psastro", 5, "toFPA: %f %f (%f,%f),(%f,%f)\n", - chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0], - chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1], - chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]); - - psTrace ("psastro", 5, "frFPA: %f %f (%f,%f),(%f,%f)\n", - chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0], - chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1], - 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); - } + 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 + + // scale from FPA to TPA (microns / pixel) + double pdelt1 = cdelt1*PS_RAD_DEG / plateScale; + double pdelt2 = cdelt2*PS_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); + + if (fpa->toSky == NULL) { + fpa->toTPA = psPlaneDistortIdentity (1); + fpa->fromTPA = psPlaneDistortIdentity (1); + fpa->toTPA->x->coeff[1][0][0][0] = pdelt1; + fpa->toTPA->y->coeff[0][1][0][0] = pdelt2; + fpa->fromTPA->x->coeff[1][0][0][0] = 1.0 / pdelt1; + fpa->fromTPA->y->coeff[0][1][0][0] = 1.0 / pdelt2; + fpa->toSky = toSky; + } else { + if (fpa->toTPA == NULL) + psAbort ("wcs", "projection defined, tangent-plane not defined"); + if (fpa->fromTPA == NULL) + psAbort ("wcs", "projection defined, tangent-plane not defined"); + + // convert from pixels on this chip to pixels on reference chip + // rX has units of refpixels / pixel + rX = pdelt1 / fpa->toTPA->x->coeff[1][0][0][0]; + rY = pdelt2 / fpa->toTPA->y->coeff[0][1][0][0]; + for (int i = 0; i <= fitOrder; i++) { + for (int j = 0; j <= fitOrder; j++) { + toFPA->x->coeff[i][j] *= rX; + toFPA->y->coeff[i][j] *= rY; + } + } + + // adjust for common toSky, toTPA for mosaic: + // find the FPA coordinate of 0,0 for this chip. + psPlane *chip = psPlaneAlloc(); + psPlane *fp = psPlaneAlloc(); + psPlane *tp = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + chip->x = chip->y = 0; + + 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 + psPlaneDistortApply (fp, fpa->fromTPA, tp, 0.0, 0.0); + + toFPA->x->coeff[0][0] = fp->x; + toFPA->y->coeff[0][0] = fp->y; + + psFree (fp); + psFree (sky); + psFree (chip); + psFree (toSky); + } + + chip->toFPA = toFPA; + // XXX this needs to perform the full (non-linear) inversion + chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA); + + // this can take a very long time... + while (fpa->toSky->R < 0) + fpa->toSky->R += 2.0*M_PI; + while (fpa->toSky->R > 2.0*M_PI) + fpa->toSky->R -= 2.0*M_PI; + + // remove the correction to the common plate scale + // NOTE: this assumes 1) we are reading in headers generated using per-chip astrometry + // and 2) we are going to measure the mosaic distortion in the next step. + // XXX perhaps make this its own function? (I'll need to store rX somewhere). + if (isMosaic) { + chip->toFPA->x->coeff[0][0] /= rX; + chip->toFPA->x->coeff[1][0] /= rX; + chip->toFPA->x->coeff[0][1] /= rX; + chip->toFPA->y->coeff[0][0] /= rY; + chip->toFPA->y->coeff[1][0] /= rY; + chip->toFPA->y->coeff[0][1] /= rY; + } + + psTrace ("psastro", 5, "toFPA: %f %f (%f,%f),(%f,%f)\n", + chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0], + chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1], + chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]); + + psTrace ("psastro", 5, "frFPA: %f %f (%f,%f),(%f,%f)\n", + chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0], + chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1], + 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); + return true; } @@ -361,13 +361,13 @@ /* discussion of the coord transformations: - X,Y: coord on a chip in pixels - L,M: coord on the focal plane (pixels) - P,Q: coord in the tangent plane (microns or mm?) - R,D: coord on the sky - - this function creates WCS terms which convert directly from chip to sky. - this function requires a linear, unrotated toTPA distortion term - toTPA->x,y->coeff[1][0],[0][1] defines the detector scale (microns / pixel) - tpSky->Xs,Ys defines the plate scale (radians / micron) + X,Y: coord on a chip in pixels + L,M: coord on the focal plane (pixels) + P,Q: coord in the tangent plane (microns or mm?) + R,D: coord on the sky + + this function creates WCS terms which convert directly from chip to sky. + this function requires a linear, unrotated toTPA distortion term + toTPA->x,y->coeff[1][0],[0][1] defines the detector scale (microns / pixel) + tpSky->Xs,Ys defines the plate scale (radians / micron) */ @@ -720,2 +720,20 @@ return status; } + +/***** + +For mosaic astrometry, we need to have a starting set of projection terms in which the +chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or +microns). This set of projections, coupled with an identity toTPA (ie, no distortion) will +result in substantial errors between the observed and predicted star positions on the focal +plane: this is the measurement of the optical distortion in the camera. At the same time, +we need to carry around the transformations which allow us to make an accurate calculation +of the position of the stars based on the input (per-chip) astrometry. These +transformations will allow us to match the raw and ref stars robustly. To convert the +per-chip astrometry (which may have been calculated with a different plate scale for each +chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust +the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for +each chip (eg, cdelt1). Thus, we need to carry around both the + +*****/ +