Index: trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10643) +++ trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10712) @@ -7,6 +7,6 @@ * @author EAM, IfA * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-12 08:00:44 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-14 06:26:15 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -107,5 +107,5 @@ // interpret header WCS keywords (only handles traditional WCS for the moment) -pmAstromWCS *pmAstromWCSfromHeader (psMetadata *header) +pmAstromWCS *pmAstromWCSfromHeader (const psMetadata *header) { psProjectionType type; @@ -234,192 +234,14 @@ } -// 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, - * wcs->trans, 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) - * plateScale (radians / physical TPA units) - * - * 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 - * 3) there is an allowed Lo,Mo term ([0][0] coefficients) - */ - - // 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 - * 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 *tmp; - - psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcs->trans->x); - psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcs->trans->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] = 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*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 (wcs->toSky->R, wcs->toSky->D, plateScale, plateScale, wcs->toSky->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); - 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); - - 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]); - - psFree (wcs->trans); - - return true; -} - -// convert toFPA / toSky components to traditional WCS -// plateScale is nominal physical scale on tangent plane (microns / arcsecond) -// this requires toTP to be the identity transformation -bool pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale) -{ - - // techinically, we can have a plate scale here (toTPA:dx,dy != 1) - if (!psPlaneDistortIsIdentity (toTPA)) - psAbort ("psastro", "invalid TPA transformation"); - - // XXX require toFPA->x->nX == toFPA->x->nY - // XXX require toFPA->y->nX == toFPA->y->nY - // XXX require toFPA->x->nX == toFPA->y->nX - // XXX require toFPA->nX == 1,2,3 - - switch (toSky->type) { +// convert wcs transformations into header WCS keywords (only handles traditional WCS for the moment) +// wcs->trans defines the transformation from pixels to degrees. +// wcs->cdelt1,2 carries the original pixels scale. +// XXX force PC00i00j to be normalized, or use cdelt1,2 to set the scale? +// here I've chosen to force the rotation matrix to be normalized +bool pmAstromWCStoHeader (psMetadata *header, const pmAstromWCS *wcs) +{ + char name[16]; // used to store FITS keyword below (always < 8, so 16 should be safe!) + + switch (wcs->toSky->type) { case PS_PROJ_SIN: psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---SIN"); @@ -439,46 +261,259 @@ break; default: - psLogMsg ("psastro", 2, "warning: unknown projection type %d\n", toSky->type); - return false; - } - - # if (0) - // XXX not really right: needs to deal with non-identity to coeffs - // XXX actually, totally wrong. fix the conversions - // XXX need to handle the plateScale - - /* 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) - */ - - // solve for CDELT1,2 (degrees / pixel) - 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) - // solve for CRPIX1,2 (Xo,Yo) : L,M(Xo,Yo) = 0,0 - - // linear solution for Xo,Yo: - xcoeff = toFPA->x->coeff; - ycoeff = toFPA->y->coeff; - R = (xcoeff[1][0]*ycoeff[0][1] - xcoeff[0][1]*ycoeff[1][0]); - Xo = det*(ycoeff[0][0]*xcoeff[0][1] - xcoeff[0][0]*ycoeff[0][1]); - Yo = det*(xcoeff[0][0]*ycoeff[1][0] - ycoeff[0][0]*xcoeff[1][0]); - - if (toFPA->x->nX > 1) { - - psPolynomial2D *XdX = psPolynomial2D_dX(toFPA->x); - psPolynomial2D *XdY = psPolynomial2D_dY(toFPA->x); - - psPolynomial2D *YdX = psPolynomial2D_dX(toFPA->y); - psPolynomial2D *YdY = psPolynomial2D_dY(toFPA->y); + psLogMsg ("psastro", 2, "warning: unknown projection type %d\n", wcs->toSky->type); + return false; + } + + psMetadataAddF64 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", wcs->toSky->R*PM_DEG_RAD); + psMetadataAddF64 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", wcs->toSky->D*PM_DEG_RAD); + + psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", wcs->crpix1); + psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", wcs->crpix2); + + // XXX make it optional to write out CDi_j terms, or other versions + // solve for CDELT1,2 (degrees / pixel) + double cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->y->coeff[1][0]); + double cdelt2 = hypot (wcs->trans->x->coeff[0][1], wcs->trans->y->coeff[0][1]); + psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1); + psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2); + + // test the PC00i00j varient: + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1 + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2 + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1 + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2 + + // Elixir-style polynomial terms + // XXX currently, Elixir/DVO cannot accept mixed orders + // XXX need to respect the masks + int fitOrder = wcs->trans->x->nX; + if (fitOrder > 1) { + for (int i = 0; i <= fitOrder; i++) { + for (int j = 0; j <= fitOrder; j++) { + if (i + j < 2) + continue; + if (i + j > fitOrder) + continue; + sprintf (name, "PCA1X%1dY%1d", i, j); + psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j)); + sprintf (name, "PCA2X%1dY%1d", i, j); + psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(wcs->cdelt1, i) / pow(wcs->cdelt2, j)); + } + } + psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder); + } + + return (true); +} + +// 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, const 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, + * wcs->trans, 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) + * plateScale (radians / physical TPA units) + * + * 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 + * 3) there is an allowed Lo,Mo term ([0][0] coefficients) + */ + + // 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 + * 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 *tmp; + + psPolynomial2D *xPx = psPolynomial2DCopy (NULL, wcs->trans->x); + psPolynomial2D *yPx = psPolynomial2DCopy (NULL, wcs->trans->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] = 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*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 (wcs->toSky->R, wcs->toSky->D, plateScale, plateScale, wcs->toSky->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); + 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); + + 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 + // XXX we need to pull the region from the chip metadata + psRegion region = psRegionSet (0, 4000, 0, 4000); + chip->fromFPA = psPlaneTransformInvert(NULL, toFPA, region, 50); + + // 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]); + + psFree (wcs); + + return true; +} + +// convert toFPA / toSky components to pmAstromWCS, then write to the headers +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +// this requires toTP to be the identity transformation +bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip) +{ + // techinically, we can have a plate scale here (fpa->toTPA:dx,dy != 1) + if (!psPlaneDistortIsIdentity (fpa->toTPA)) + psAbort ("psastro", "invalid TPA transformation"); + + // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY + // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY + // XXX require chip->toFPA->x->nX == chip->toFPA->y->nX + // XXX require chip->toFPA->nX == 1,2,3 + + int fitOrder = chip->toFPA->x->nX; + pmAstromWCS *wcs = pmAstromWCSAlloc(fitOrder, fitOrder); + + // convert projection from TPA to SKY into wcs projection (degrees to radians) + wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, fpa->toSky->type); + wcs->crval1 = fpa->toSky->R*PM_DEG_RAD; + wcs->crval2 = fpa->toSky->D*PM_DEG_RAD; + + // crpix1,2 = X,Y(crval1,2) + // start with linear solution for Xo,Yo: + double R = (chip->toFPA->x->coeff[1][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->x->coeff[1][0]); + double Xo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1])/R; + double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0])/R; + + // iterate to actual solution: requires small non-linear terms + if (fitOrder > 1) { + psPolynomial2D *XdX = psPolynomial2D_dX(NULL, chip->toFPA->x); + psPolynomial2D *XdY = psPolynomial2D_dY(NULL, chip->toFPA->x); + + psPolynomial2D *YdX = psPolynomial2D_dX(NULL, chip->toFPA->y); + psPolynomial2D *YdY = psPolynomial2D_dY(NULL, chip->toFPA->y); psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32); @@ -494,6 +529,6 @@ Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo); - Beta->data.F32[0] = psPolynomial2DEval (toFPA->x, Xo, Yo); - Beta->data.F32[1] = psPolynomial2DEval (toFPA->y, Xo, Yo); + Beta->data.F32[0] = psPolynomial2DEval (chip->toFPA->x, Xo, Yo); + Beta->data.F32[1] = psPolynomial2DEval (chip->toFPA->y, Xo, Yo); psMatrixGJSolveF32 (Alpha, Beta); @@ -502,85 +537,61 @@ Yo += Beta->data.F32[1]; } - } - - psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", Xo); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", Yo); - - psPolynomial2D *xWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->x->nX, toFPA->x->nY); - psPolynomial2D *yWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->y->nX, toFPA->y->nY); - - psPolynomial2D *xPx = psPolynomial2DCopy (toFPA->x); - psPolynomial2D *yPx = psPolynomial2DCopy (toFPA->y); - - // skip the zero order terms - // XXX double check that these relationships are correct - for (int i = 0; i < toFPA->x->nX; i++) { - psPolynomial2D *xPy = psPolynomial2DCopy (xPx); - psPolynomial2D *yPy = psPolynomial2DCopy (yPx); - for (int j = 0; j < toFPA->x->nY; j++) { - xWCS->coords[i][j] = psPolynomial2DEval (xPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j); - yWCS->coords[i][j] = psPolynomial2DEval (yPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j); - psPolynomial2D_dY(xPy, xPy); - psPolynomial2D_dY(yPy, yPy); - } - psPolynomial2D_dX(xPx, xPx); - psPolynomial2D_dX(yPx, yPx); - } - - while (coords[0].crval1 < 0) - coords[0].crval1 += 360.0; - while (coords[0].crval1 > 360.0) - coords[0].crval1 -= 360.0; - - 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); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2); - - psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", xWCS->coeff[1][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", xWCS->coeff[0][1]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", yWCS->coeff[1][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", yWCS->coeff[0][1]); - - // XXX respect the masks - for (int i = 0; i < xWCS->nX; i++) { - for (int j = 0; j < xWCS->nX; j++) { - if (i + j < 2) - continue; - sprintf (name, "PCA1dX%1dY%1d", i, j); - psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", xWCS->coeff[i][j]); - } - } - for (int i = 0; i < yWCS->nX; i++) { - for (int j = 0; j < yWCS->nX; j++) { - if (i + j < 2) - continue; - sprintf (name, "PCA2dX%1dY%1d", i, j); - psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", yWCS->coeff[i][j]); - } - } - - # else - - 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*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); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale); - - # endif - - // alternative representations use - // CD1_1 = PC001001*CDELT1, etc - // make these representations optional - + psFree (Alpha); + psFree (Beta); + psFree (XdX); + psFree (XdY); + psFree (YdX); + psFree (YdY); + } + wcs->crpix1 = Xo; + wcs->crpix2 = Yo; + + // convert the chip->toFPA polynomials (with 0,0 ref) into wcs polynomials, with Xo,Yo ref + // chip->toFPA(x,y) = wcs->trans(x-xo,y-yo) -- see comment in pmAstromReadWCS + psPolynomial2D *tmp; + psPolynomial2D *xPx = psPolynomial2DCopy (NULL, chip->toFPA->x); + psPolynomial2D *yPx = psPolynomial2DCopy (NULL, chip->toFPA->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++) { + wcs->trans->x->mask[i][j] = chip->toFPA->x->mask[i][j]; + wcs->trans->y->mask[i][j] = chip->toFPA->y->mask[i][j]; + wcs->trans->x->coeff[i][j] = (wcs->trans->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, wcs->crpix1, wcs->crpix2) / tgamma(i+1) / tgamma(j+1); + wcs->trans->y->coeff[i][j] = (wcs->trans->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; + } + + // cdelt1,2 convert from pixels->degrees + double cdelt1 = fpa->toTPA->x->coeff[1][0][0][0]*fpa->toSky->Xs*PM_DEG_RAD; + double cdelt2 = fpa->toTPA->y->coeff[0][1][0][0]*fpa->toSky->Ys*PM_DEG_RAD; + + // convert wcs->trans to a matrix which yields L,M in pixels + 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; + } + } + + pmAstromWCStoHeader (header, wcs); + + psFree (wcs); return true; } @@ -841,2 +852,14 @@ # endif +/* 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) +*/ +