Index: /trunk/psModules/src/astrom/Makefile.am =================================================================== --- /trunk/psModules/src/astrom/Makefile.am (revision 10824) +++ /trunk/psModules/src/astrom/Makefile.am (revision 10825) @@ -6,4 +6,5 @@ pmAstrometryObjects.c \ pmAstrometryDistortion.c \ + pmAstrometryUtils.c \ pmAstrometryWCS.c @@ -11,4 +12,5 @@ pmAstrometryObjects.h \ pmAstrometryDistortion.h \ + pmAstrometryUtils.h \ pmAstrometryWCS.h Index: /trunk/psModules/src/astrom/pmAstrometryDistortion.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryDistortion.c (revision 10824) +++ /trunk/psModules/src/astrom/pmAstrometryDistortion.c (revision 10825) @@ -7,6 +7,6 @@ * @author EAM, IfA * -* @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-10 02:06:47 $ +* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-22 21:23:06 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -187,5 +187,5 @@ for (int i = 0; i < fpa->toTPA->x->nX; i++) { for (int j = 0; j < fpa->toTPA->x->nY; j++) { - if (fpa->toTPA->x->mask[i][j][0][0]) { + if (fpa->toTPA->x->mask[i][j]) { localX->mask[i-1][j] = 1; localY->mask[i][j-1] = 1; @@ -203,19 +203,19 @@ // update fpa->toTP distortion terms - fpa->toTPA->x->coeff[0][0][0][0] = 0; + fpa->toTPA->x->coeff[0][0] = 0; for (int i = 1; i < fpa->toTPA->x->nX; i++) { - if (!fpa->toTPA->x->mask[i][0][0][0]) { - fpa->toTPA->x->coeff[i][0][0][0] = localX->coeff[i-1][0] / i; + if (!fpa->toTPA->x->mask[i][0]) { + fpa->toTPA->x->coeff[i][0] = localX->coeff[i-1][0] / i; } } for (int j = 1; j < fpa->toTPA->x->nY; j++) { - if (!fpa->toTPA->x->mask[0][j][0][0]) { - fpa->toTPA->x->coeff[0][j][0][0] = localY->coeff[0][j-1] / j; + if (!fpa->toTPA->x->mask[0][j]) { + fpa->toTPA->x->coeff[0][j] = localY->coeff[0][j-1] / j; } } for (int i = 1; i < fpa->toTPA->x->nX; i++) { for (int j = 1; j < fpa->toTPA->x->nY; j++) { - if (!fpa->toTPA->x->mask[i][j][0][0]) { - fpa->toTPA->x->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j); + if (!fpa->toTPA->x->mask[i][j]) { + fpa->toTPA->x->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j); } } @@ -232,5 +232,5 @@ for (int i = 0; i < fpa->toTPA->y->nX; i++) { for (int j = 0; j < fpa->toTPA->y->nY; j++) { - if (fpa->toTPA->y->mask[i][j][0][0]) { + if (fpa->toTPA->y->mask[i][j]) { localX->mask[i-1][j] = 1; localY->mask[i][j-1] = 1; @@ -259,19 +259,19 @@ // update fpa->toTP distortion terms - fpa->toTPA->y->coeff[0][0][0][0] = 0; + fpa->toTPA->y->coeff[0][0] = 0; for (int i = 1; i < fpa->toTPA->y->nX; i++) { - if (!fpa->toTPA->y->mask[i][0][0][0]) { - fpa->toTPA->y->coeff[i][0][0][0] = localX->coeff[i-1][0] / i; + if (!fpa->toTPA->y->mask[i][0]) { + fpa->toTPA->y->coeff[i][0] = localX->coeff[i-1][0] / i; } } for (int j = 1; j < fpa->toTPA->y->nY; j++) { - if (!fpa->toTPA->y->mask[0][j][0][0]) { - fpa->toTPA->y->coeff[0][j][0][0] = localY->coeff[0][j-1] / j; + if (!fpa->toTPA->y->mask[0][j]) { + fpa->toTPA->y->coeff[0][j] = localY->coeff[0][j-1] / j; } } for (int i = 1; i < fpa->toTPA->y->nX; i++) { for (int j = 1; j < fpa->toTPA->y->nY; j++) { - if (!fpa->toTPA->y->mask[i][j][0][0]) { - fpa->toTPA->y->coeff[i][j][0][0] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j); + if (!fpa->toTPA->y->mask[i][j]) { + fpa->toTPA->y->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j); } } Index: /trunk/psModules/src/astrom/pmAstrometryUtils.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 10825) +++ /trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 10825) @@ -0,0 +1,339 @@ +/** @file pmAstrometryUtils.c + * + * @brief utility functions for transform and distort functions + * + * @ingroup Astrometry + * + * @author EAM, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-22 21:23:06 $ + * + * Copyright 2006 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include "pmFPA.h" +#include "pmAstrometryUtils.h" + +// given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y +// for which L,M = 0,0. tol is the allowed error on x,y. +psPlane *psPlaneTransformGetCenter (psPlaneTransform *trans, float tol) +{ + + // crpix1,2 = X,Y(crval1,2) + // start with linear solution for Xo,Yo: + double R = (trans->x->coeff[1][0]*trans->y->coeff[0][1] - trans->x->coeff[0][1]*trans->y->coeff[1][0]); + double Xo = (trans->y->coeff[0][0]*trans->x->coeff[0][1] - trans->x->coeff[0][0]*trans->y->coeff[0][1])/R; + double Yo = (trans->x->coeff[0][0]*trans->y->coeff[1][0] - trans->y->coeff[0][0]*trans->x->coeff[1][0])/R; + + // iterate to actual solution: requires small non-linear terms + if (trans->x->nX > 1) { + psPolynomial2D *XdX = psPolynomial2D_dX(NULL, trans->x); + psPolynomial2D *XdY = psPolynomial2D_dY(NULL, trans->x); + + psPolynomial2D *YdX = psPolynomial2D_dX(NULL, trans->y); + psPolynomial2D *YdY = psPolynomial2D_dY(NULL, trans->y); + + psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32); + psVector *Beta = psVectorAlloc (2, PS_DATA_F32); + + /* this loop uses the Newton-Raphson method to solve for Xo,Yo + * it needs the high order terms to be small + * Xo,Yo are in pixels; + */ + float dPos = tol + 1; + for (int i = 0; (dPos > tol) && (i < 10); i++) { + // NOTE: order for Alpha is: [y][x] + Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo); + Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo); + Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo); + Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo); + + Beta->data.F32[0] = psPolynomial2DEval (trans->x, Xo, Yo); + Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo); + + psMatrixGJSolveF32 (Alpha, Beta); + + Xo -= Beta->data.F32[0]; + Yo -= Beta->data.F32[1]; + dPos = hypot(Xo,Yo); + } + psFree (Alpha); + psFree (Beta); + psFree (XdX); + psFree (XdY); + psFree (YdX); + psFree (YdY); + } + psPlane *center = psPlaneAlloc (); + center->x = Xo; + center->y = Yo; + + return center; +} + +// convert a transformation L(x,y) to L'(x-xo,y-yo) +psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, float Xo, float Yo) +{ + + // validate fit order + + if (output == NULL) { + output = psPlaneTransformAlloc(input->x->nX, input->x->nY); + } + + /* 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. + */ + + psPolynomial2D *tmp; + + psPolynomial2D *xPx = psPolynomial2DCopy (NULL, input->x); + psPolynomial2D *yPx = psPolynomial2DCopy (NULL, input->y); + + for (int i = 0; i <= input->x->nX; i++) { + psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx); + psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx); + for (int j = 0; j <= input->x->nY; j++) { + output->x->mask[i][j] = input->x->mask[i][j]; + output->y->mask[i][j] = input->y->mask[i][j]; + output->x->coeff[i][j] = (output->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1); + output->y->coeff[i][j] = (output->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, Xo, Yo) / 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; + } + return output; +} + +// construct a psPlaneTransform which is the identify transformation +psPlaneTransform *psPlaneTransformIdentity (int order) +{ + + psPlaneTransform *transform; + + if (order < 1) + psAbort ("psastro", "invalid order"); + if (order > 3) + psAbort ("psastro", "invalid order"); + + // all coeffs and masks initially set to 0 + transform = psPlaneTransformAlloc (order, order); + + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + transform->x->mask [i][j] = 1; + transform->y->mask [i][j] = 1; + } + } + } + transform->x->coeff[1][0] = 1; + transform->y->coeff[0][1] = 1; + transform->x->mask[1][0] = 0; + transform->y->mask[0][1] = 0; + + return transform; +} + +// check that the given psPlaneTransform is the identity * (Xs,Ys) +bool psPlaneTransformIsDiagonal (psPlaneTransform *transform) +{ + + int order; + bool status; + + // we currently only support up to 3rd order polynomials + if (transform->x->nX < 1) + return false; + if (transform->x->nY < 1) + return false; + if (transform->y->nX < 1) + return false; + if (transform->y->nY < 1) + return false; + + if (transform->x->nX != transform->x->nY) + return false; + if (transform->y->nX != transform->y->nY) + return false; + + // these are not actually valid tests + if (transform->x->nX > 3) + return false; + if (transform->x->nY > 3) + return false; + if (transform->y->nX > 3) + return false; + if (transform->y->nY > 3) + return false; + + status = true; + order = transform->x->nX; + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + // high-order cross terms must be masked (eg, x^3 y^2) + status &= transform->x->mask[i][j]; + } else { + status &= !transform->x->mask[i][j]; + if ((i == 1) && (i + j == 1)) { + // linear, diagonal terms must be non-zero + status &= (fabs(transform->x->coeff[i][j]) > FLT_EPSILON); + } else { + // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) + status &= (fabs(transform->x->coeff[i][j]) < FLT_EPSILON); + } + } + } + } + + order = transform->y->nX; + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + // high-order cross terms must be masked (eg, x^3 y^2) + status &= transform->y->mask[i][j]; + } else { + status &= !transform->y->mask[i][j]; + if ((j == 1) && (i + j == 1)) { + // linear, diagonal terms must be 1.0 + status &= (fabs(transform->y->coeff[i][j]) > FLT_EPSILON); + } else { + // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) + status &= (fabs(transform->y->coeff[i][j]) < FLT_EPSILON); + } + } + } + } + return status; +} + +// construct a psPlaneDistort which is the identify transformation +psPlaneDistort *psPlaneDistortIdentity (int order) +{ + + psPlaneDistort *distort; + + if (order < 1) + psAbort ("psastro", "invalid order"); + if (order > 3) + psAbort ("psastro", "invalid order"); + + // all coeffs and masks initially set to 0 + distort = psPlaneDistortAlloc (order, order, 0, 0); + + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + distort->x->mask [i][j][0][0] = 1; + distort->y->mask [i][j][0][0] = 1; + } + } + } + distort->x->coeff[1][0][0][0] = 1; + distort->y->coeff[0][1][0][0] = 1; + + return distort; +} + +// check that the given psPlaneDistort is the identity * (Xs,Ys) +bool psPlaneDistortIsDiagonal (psPlaneDistort *distort) +{ + + int order; + bool status; + + // we currently only support up to 3rd order polynomials + if (distort->x->nX < 1) + return false; + if (distort->x->nY < 1) + return false; + if (distort->y->nX < 1) + return false; + if (distort->y->nY < 1) + return false; + + if (distort->x->nX > 3) + return false; + if (distort->x->nY > 3) + return false; + if (distort->y->nX > 3) + return false; + if (distort->y->nY > 3) + return false; + + if (distort->x->nZ > 0) + return false; + if (distort->x->nT > 0) + return false; + if (distort->y->nZ > 0) + return false; + if (distort->y->nT > 0) + return false; + + if (distort->x->nX != distort->x->nY) + return false; + if (distort->y->nX != distort->y->nY) + return false; + + status = true; + order = distort->x->nX; + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + // high-order cross terms must be masked (eg, x^3 y^2) + status &= distort->x->mask[i][j][0][0]; + } else { + status &= !distort->x->mask[i][j][0][0]; + if ((i == 1) && (i + j == 1)) { + // linear, diagonal terms must be 1.0 + status &= (fabs(distort->x->coeff[i][j][0][0]) > FLT_EPSILON); + } else { + // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) + status &= (fabs(distort->x->coeff[i][j][0][0]) < FLT_EPSILON); + } + } + } + } + + order = distort->y->nX; + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + // high-order cross terms must be masked (eg, x^3 y^2) + status &= distort->y->mask[i][j][0][0]; + } else { + status &= !distort->y->mask[i][j][0][0]; + if ((j == 1) && (i + j == 1)) { + // linear, diagonal terms must be 1.0 + status &= (fabs(distort->y->coeff[i][j][0][0]) > FLT_EPSILON); + } else { + // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) + status &= (fabs(distort->y->coeff[i][j][0][0]) < FLT_EPSILON); + } + } + } + } + return status; +} Index: /trunk/psModules/src/astrom/pmAstrometryUtils.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryUtils.h (revision 10825) +++ /trunk/psModules/src/astrom/pmAstrometryUtils.h (revision 10825) @@ -0,0 +1,27 @@ +/** @file pmAstrometryUtils.h +* +* @brief utility functions for transform and distort functions +* +* @ingroup Astrometry +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-22 21:23:06 $ +* +* Copyright 2006 Institute for Astronomy, University of Hawaii +*/ + +#ifndef PM_ASTROMETRY_UTILS_H +#define PM_ASTROMETRY_UTILS_H + +psPlane *psPlaneTransformGetCenter (psPlaneTransform *trans, float tol); +psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, float Xo, float Yo); + +psPlaneTransform *psPlaneTransformIdentity (int order); +psPlaneDistort *psPlaneDistortIdentity (int order); + +bool psPlaneTransformIsDiagonal (psPlaneTransform *transform); +bool psPlaneDistortIsDiagonal (psPlaneDistort *distort); + +#endif Index: /trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10824) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10825) @@ -7,6 +7,6 @@ * @author EAM, IfA * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-20 16:37:55 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-22 21:23:06 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -20,4 +20,5 @@ #include "pmFPA.h" #include "pmAstrometryWCS.h" +#include "pmAstrometryUtils.h" // the following functions support coordinate transformations direcly related to the FITS WCS @@ -27,30 +28,96 @@ // 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 +// not represent a simple 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); - psFree (wcs->toSky); -} - -pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder) -{ - - pmAstromWCS *wcs = (pmAstromWCS *) psAlloc(sizeof(pmAstromWCS)); - psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree); - - wcs->trans = psPlaneTransformAlloc (nXorder, nYorder); - wcs->toSky = NULL; - - memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE); - memset (wcs->ctype2, 0, PM_ASTROM_WCS_TYPE_SIZE); - return wcs; -} - +// 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) +{ + pmAstromWCS *wcs = pmAstromWCSfromHeader (header); + if (!wcs) { + return false; + } + + bool status = pmAstromWCStoFPA (fpa, chip, wcs, plateScale); + + psFree (wcs); + return status; +} + +// convert toFPA / toSky components to pmAstromWCS +// tolerance is convergence for inversion of non-linear terms in pixels +bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol) +{ + // 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 + + // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1) + if (!psPlaneTransformIsDiagonal (fpa->toTPA)) + psAbort ("psastro", "invalid TPA transformation"); + + pmAstromWCS *wcs = pmAstromWCSfromFPA(fpa, chip, tol); + pmAstromWCStoHeader (header, wcs); + + psFree (wcs); + return true; +} + +// interpret chip header WCS as bilevel chip components +// plateScale is nominal physical scale on tangent plane (radians / TPA physical units) +bool pmAstromReadBilevelChip (pmChip *chip, const psMetadata *header) +{ + pmAstromWCS *wcs = pmAstromWCSfromHeader (header); + if (!wcs) { + return false; + } + + bool status = pmAstromWCSBileveltoChip (chip, wcs); + + psFree (wcs); + return status; +} + +// convert toFPA / toSky components to traditional WCS +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header) +{ + pmAstromWCS *wcs = pmAstromWCSfromHeader (header); + if (!wcs) { + return false; + } + + bool status = pmAstromWCSBileveltoFPA (fpa, wcs); + + psFree (wcs); + return status; +} + +// convert chip->toFPA components to bilevel WCS +bool pmAstromWriteBilevelChip (psMetadata *header, const pmChip *chip, float tol) +{ + pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (chip, tol); + + pmAstromWCStoHeader (header, wcs); + + psFree (wcs); + return true; +} + + +// convert fpa->toTPA, fpa->toSky components to bilevel WCS +bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, float tol) +{ + pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol); + + pmAstromWCStoHeader (header, wcs); + + psFree (wcs); + return true; +} + +// convert coordinates from chip to sky using a pmAstromWCS structure bool pmAstromWCStoSky (psSphere *sky, pmAstromWCS *wcs, psPlane *chip) { @@ -77,4 +144,5 @@ } +// convert coordinates from sky to chip using a pmAstromWCS structure bool pmAstromWCStoChip (psPlane *chip, pmAstromWCS *wcs, psSphere *sky) { @@ -106,5 +174,5 @@ } -// interpret header WCS keywords (only handles traditional WCS for the moment) +// interpret header WCS keywords (valid for bilevel and traditional WCS) pmAstromWCS *pmAstromWCSfromHeader (const psMetadata *header) { @@ -121,14 +189,7 @@ // determine projection type - // XXX add the Elixir DIS / WRP two-layer projection here - type = PS_PROJ_NTYPE; - if (!strcmp (&ctype[4], "-SIN")) - type = PS_PROJ_SIN; - if (!strcmp (&ctype[4], "-TAN")) - type = PS_PROJ_TAN; - if (!strcmp (&ctype[4], "-AIT")) - type = PS_PROJ_AIT; - if (!strcmp (&ctype[4], "-PAR")) - type = PS_PROJ_PAR; + // XXX there are two indications for higher-order terms: the type (DIS,WRP,PLY,ZPL) and + // the value of NPLYTERM. + type = psProjectTypeFromString (ctype); if (type == PS_PROJ_NTYPE) { psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype); @@ -245,26 +306,13 @@ { 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"); - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--SIN"); - break; - case PS_PROJ_TAN: - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN"); - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--TAN"); - break; - case PS_PROJ_AIT: - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---AIT"); - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--AIT"); - break; - case PS_PROJ_PAR: - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---PAR"); - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--PAR"); - break; - default: - psLogMsg ("psastro", 2, "warning: unknown projection type %d\n", wcs->toSky->type); - return false; - } + char *type; + + type = psProjectTypeToString (wcs->toSky->type, "RA--"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", type); + psFree (type); + + type = psProjectTypeToString (wcs->toSky->type, "DEC-"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", type); + psFree (type); psMetadataAddF64 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", wcs->toSky->R*PM_DEG_RAD); @@ -312,5 +360,5 @@ // interpret header WCS (only handles traditional WCS for the moment) // plateScale is nominal physical scale on tangent plane (radians / TPA physical units) -bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale, bool isMosaic) +bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale) { psPlaneTransform *toFPA; @@ -326,5 +374,5 @@ * * 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 + * 1) the output is in pixels (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) @@ -341,44 +389,6 @@ } - // 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; - } + // create transformation with 0,0 reference pixel + toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2); // scale from FPA to TPA (microns / pixel) @@ -392,12 +402,17 @@ 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->toTPA = psPlaneTransformIdentity (1); + fpa->fromTPA = psPlaneTransformIdentity (1); + fpa->toTPA->x->coeff[1][0] = pdelt1; + fpa->toTPA->y->coeff[0][1] = pdelt2; + fpa->fromTPA->x->coeff[1][0] = 1.0 / pdelt1; + fpa->fromTPA->y->coeff[0][1] = 1.0 / pdelt2; fpa->toSky = toSky; } else { + + // this section allows the loaded chip to be included in an fpa structure in which + // other chips have already been loaded (ie, the fpa->toTPA, fpa->toSky components have + // already been defined). we have to adjust to match the existing transformation. + if (fpa->toTPA == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined"); @@ -407,8 +422,8 @@ // 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++) { + rX = pdelt1 / fpa->toTPA->x->coeff[1][0]; + rY = pdelt2 / fpa->toTPA->y->coeff[0][1]; + for (int i = 0; i <= fpa->toTPA->x->nX; i++) { + for (int j = 0; j <= fpa->toTPA->x->nY; j++) { toFPA->x->coeff[i][j] *= rX; toFPA->y->coeff[i][j] *= rY; @@ -425,8 +440,8 @@ 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); + psPlaneTransformApply (tp, fpa->toTPA, fp); 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); + psPlaneTransformApply (fp, fpa->fromTPA, tp); toFPA->x->coeff[0][0] = fp->x; @@ -454,17 +469,4 @@ 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], @@ -480,11 +482,54 @@ } +// convert a pmAstromWCS structure representing a bilevel chip into corresponding chip elements +bool pmAstromWCSBileveltoChip (pmChip *chip, const pmAstromWCS *wcs) +{ + /* 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) + */ + + // create transformation with 0,0 reference pixel + chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2); + + // XXX this needs to perform the full (non-linear) inversion + // XXX we need to pull the region from the chip metadata + // is chip trimsec defined? do we need to ensure the + // XXX should this function be defined for the CELL, not the CHIP? + // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC"); + psRegion region = psRegionSet (0, 4000, 0, 4000); + chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, region, 50); + + return true; +} + +// convert a pmAstromWCS structure representing a bilevel mosaic into corresponding fpa elements +bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs) +{ + // projection from TPA to SKY + fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->toSky->Xs, wcs->toSky->Ys, wcs->toSky->type); + + // create transformation with 0,0 reference pixel + fpa->toTPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2); + + // XXX this needs to perform the full (non-linear) inversion + // XXX we need to pull the region from the chip metadata + // is chip trimsec defined? do we need to ensure the + // XXX should this function be defined for the CELL, not the CHIP? + // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC"); + psRegion region = psRegionSet (0, 4000, 0, 4000); + fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 50); + + return true; +} + // convert toFPA / toSky components to pmAstromWCS -// tolerance is in pixels +// tolerance is allowed error in center solution in pixels pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, float tol) { // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1) - if (!psPlaneDistortIsDiagonal (fpa->toTPA)) + if (!psPlaneTransformIsDiagonal (fpa->toTPA)) psAbort ("psastro", "invalid TPA transformation"); @@ -494,6 +539,5 @@ // XXX require chip->toFPA->nX == 1,2,3 - int fitOrder = chip->toFPA->x->nX; - pmAstromWCS *wcs = pmAstromWCSAlloc(fitOrder, fitOrder); + pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY); // convert projection from TPA to SKY into wcs projection (degrees to radians) @@ -502,87 +546,18 @@ 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->y->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->y->coeff[1][0]); - double Xo = (chip->toFPA->y->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[0][1])/R; - double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[1][0] - chip->toFPA->y->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); - psVector *Beta = psVectorAlloc (2, PS_DATA_F32); - - // XXX this loop is rather arbitrary in length... - // XXX measure the error and use as criterion - /* this is the Newton-Raphson method to solve for Xo,Yo - it needs the high order terms to be small */ - // Xo,Yo are in pixels; - float dPos = tol + 1; - for (int i = 0; (dPos > tol) && (i < 10); i++) { - // NOTE: order is: [y][x] - Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo); - Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo); - Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo); - Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, 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); - - Xo -= Beta->data.F32[0]; - Yo -= Beta->data.F32[1]; - dPos = hypot(Xo,Yo); - } - 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; - } + // given transformation, solve for coordinates which yields output coordinates of 0,0 + psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol); + + // adjust wcs transform to use center as reference coordinate + psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y); + + // calculated center is crpix1,2 + wcs->crpix1 = center->x; + wcs->crpix2 = center->y; + psFree (center); // 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; + double cdelt1 = fpa->toTPA->x->coeff[1][0]*fpa->toSky->Xs*PM_DEG_RAD; + double cdelt2 = fpa->toTPA->y->coeff[0][1]*fpa->toSky->Ys*PM_DEG_RAD; wcs->cdelt1 = cdelt1; wcs->cdelt2 = cdelt2; @@ -598,22 +573,12 @@ } -// 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) -{ - pmAstromWCS *wcs = pmAstromWCSfromHeader (header); - if (!wcs) { - return false; - } - - bool status = pmAstromWCStoFPA (fpa, chip, wcs, plateScale, isMosaic); - - psFree (wcs); - return status; -} - -// convert toFPA / toSky components to pmAstromWCS, then write to the headers -// tolerance is in pixels -bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol) +/* the bilevel astrometry description consists of a polynomial warping from + chip coordinates to FPA coordinates (coords->ctype = LIN---WRP), followed + by a polynomial representation of the telescope distortion + the projection + (coords->ctype = RA---DIS). +*/ + +// convert the chip-level toFPA to a wcs polynomial transformation +pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, float tol) { // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY @@ -622,263 +587,88 @@ // XXX require chip->toFPA->nX == 1,2,3 - // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1) - if (!psPlaneDistortIsDiagonal (fpa->toTPA)) - psAbort ("psastro", "invalid TPA transformation"); - - pmAstromWCS *wcs = pmAstromWCSfromFPA(fpa, chip, tol); - pmAstromWCStoHeader (header, wcs); - - psFree (wcs); - return true; -} - -// convert toFPA / toSky components to traditional WCS -// plateScale is nominal physical scale on tangent plane (microns / arcsecond) -bool pmAstromWriteBilevelChip (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol) -{ - pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (fpa, chip, tol); - - pmAstromWCStoHeader (header, wcs); - - psFree (wcs); - return true; -} - - -// convert toFPA / toSky components to traditional WCS -// plateScale is nominal physical scale on tangent plane (microns / arcsecond) -bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, float tol) -{ - pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol); - - pmAstromWCStoHeader (header, wcs); - - psFree (wcs); - return true; -} - -pmAstromWCS *pmAstromWCSBilevelChipFromFPA (pmFPA *fpa, float tol) -{ - - // conv - - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---WRP"); - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--WRP"); - - // XXX not really right: needs to deal with non-identity toTP coeffs & plateScale - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toFPA->x->coeff[0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toFPA->y->coeff[0][0]); - 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, "", 1.0); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", 1.0); - - if (toFPA->x->nX != toFPA->x->nY) - psAbort ("psastro", "mis-matched tangent plane orders (1)"); - if (toFPA->x->nX != toFPA->y->nX) - psAbort ("psastro", "mis-matched tangent plane orders (2)"); - if (toFPA->x->nX != toFPA->y->nY) - psAbort ("psastro", "mis-matched tangent plane orders (3)"); - - switch (toFPA->x->nX) { - case 3: - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toFPA->x->coeff[3][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toFPA->x->coeff[2][1]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toFPA->x->coeff[1][2]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toFPA->x->coeff[0][3]); - - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toFPA->y->coeff[3][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toFPA->y->coeff[2][1]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toFPA->y->coeff[1][2]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toFPA->y->coeff[0][3]); - - case 2: - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toFPA->x->coeff[2][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toFPA->x->coeff[1][1]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toFPA->x->coeff[0][2]); - - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toFPA->y->coeff[2][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toFPA->y->coeff[1][1]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toFPA->y->coeff[0][2]); - - case 1: - 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); - break; - - case 0: - psAbort ("psastro", "invalid tangent plane order"); - } - return true; -} - -// convert toFPA / toSky components to traditional WCS -// plateScale is nominal physical scale on tangent plane (microns / arcsecond) -psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale) -{ - - psMetadata *header = psMetadataAlloc (); - - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---DIS"); - psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--DIS"); - - // XXX need to handle the plateScale?? - 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*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) - psAbort ("psastro", "mis-matched tangent plane orders (1)"); - if (toTP->x->nX != toTP->y->nX) - psAbort ("psastro", "mis-matched tangent plane orders (2)"); - if (toTP->x->nX != toTP->y->nY) - psAbort ("psastro", "mis-matched tangent plane orders (3)"); - - switch (toTP->x->nX) { - case 3: - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toTP->x->coeff[3][0][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toTP->x->coeff[2][1][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toTP->x->coeff[1][2][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toTP->x->coeff[0][3][0][0]); - - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toTP->y->coeff[3][0][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toTP->y->coeff[2][1][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toTP->y->coeff[1][2][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toTP->y->coeff[0][3][0][0]); - - case 2: - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toTP->x->coeff[2][0][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toTP->x->coeff[1][1][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toTP->x->coeff[0][2][0][0]); - - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toTP->y->coeff[2][0][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toTP->y->coeff[1][1][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toTP->y->coeff[0][2][0][0]); - - case 1: - psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toTP->x->coeff[1][0][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toTP->x->coeff[0][1][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toTP->y->coeff[1][0][0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toTP->y->coeff[0][1][0][0]); - break; - - case 0: - psAbort ("psastro", "invalid tangent plane order"); - } - return header; -} - -// construct a psPlaneDistort which is the identify transformation -psPlaneDistort *psPlaneDistortIdentity (int order) -{ - - psPlaneDistort *distort; - - if (order < 1) - psAbort ("psastro", "invalid order"); - if (order > 3) - psAbort ("psastro", "invalid order"); - - // all coeffs and masks initially set to 0 - distort = psPlaneDistortAlloc (order, order, 0, 0); - - for (int i = 0; i <= order; i++) { - for (int j = 0; j <= order; j++) { - if (i + j > order) { - distort->x->mask [i][j][0][0] = 1; - distort->y->mask [i][j][0][0] = 1; - } - } - } - distort->x->coeff[1][0][0][0] = 1; - distort->y->coeff[0][1][0][0] = 1; - - return distort; -} - -// check that the given psPlaneDistort is the identity * (Xs,Ys) -bool psPlaneDistortIsDiagonal (psPlaneDistort *distort) -{ - - int order; - bool status; - - // we currently only support up to 3rd order polynomials - if (distort->x->nX < 1) - return false; - if (distort->x->nY < 1) - return false; - if (distort->y->nX < 1) - return false; - if (distort->y->nY < 1) - return false; - - if (distort->x->nX > 3) - return false; - if (distort->x->nY > 3) - return false; - if (distort->y->nX > 3) - return false; - if (distort->y->nY > 3) - return false; - - if (distort->x->nZ > 0) - return false; - if (distort->x->nT > 0) - return false; - if (distort->y->nZ > 0) - return false; - if (distort->y->nT > 0) - return false; - - if (distort->x->nX != distort->x->nY) - return false; - if (distort->y->nX != distort->y->nY) - return false; - - status = true; - order = distort->x->nX; - for (int i = 0; i <= order; i++) { - for (int j = 0; j <= order; j++) { - if (i + j > order) { - // high-order cross terms must be masked (eg, x^3 y^2) - status &= distort->x->mask[i][j][0][0]; - } else { - status &= !distort->x->mask[i][j][0][0]; - if ((i == 1) && (i + j == 1)) { - // linear, diagonal terms must be 1.0 - status &= (fabs(distort->x->coeff[i][j][0][0]) > FLT_EPSILON); - } else { - // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) - status &= (fabs(distort->x->coeff[i][j][0][0]) < FLT_EPSILON); - } - } - } - } - - order = distort->y->nX; - for (int i = 0; i <= order; i++) { - for (int j = 0; j <= order; j++) { - if (i + j > order) { - // high-order cross terms must be masked (eg, x^3 y^2) - status &= distort->y->mask[i][j][0][0]; - } else { - status &= !distort->y->mask[i][j][0][0]; - if ((j == 1) && (i + j == 1)) { - // linear, diagonal terms must be 1.0 - status &= (fabs(distort->y->coeff[i][j][0][0]) > FLT_EPSILON); - } else { - // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) - status &= (fabs(distort->y->coeff[i][j][0][0]) < FLT_EPSILON); - } - } - } - } - return status; + // convert chip->toFPA to wcs format (WRP) + pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY); + + // Chip to FPA transformation is a Cartesian 'projection' + wcs->toSky = psProjectionAlloc (0.0, 0.0, 1.0, 1.0, PS_PROJ_WRP); + + // reference pixel for FPA is 0.0, 0.0 + wcs->crval1 = 0.0; + wcs->crval2 = 0.0; + + // given transformation, solve for coordinates which yields output coordinates of 0,0 + psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol); + + // adjust wcs transform to use center as reference coordinate + psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y); + + // calculated center is crpix1,2 + wcs->crpix1 = center->x; + wcs->crpix2 = center->y; + psFree (center); + + // output coordinates are in FPA pixels + wcs->cdelt1 = 1.0; + wcs->cdelt2 = 1.0; + + return wcs; +} + +// convert the fpa-level toTPA, toSky to a wcs polynomial transformation +pmAstromWCS *pmAstromWCSBilevelMosaicFromFPA (const pmFPA *fpa, float tol) +{ + // XXX require fpa->toTPA->x->nX == fpa->toTPA->x->nY + // XXX require fpa->toTPA->y->nX == fpa->toTPA->y->nY + // XXX require fpa->toTPA->x->nX == fpa->toTPA->y->nX + // XXX require fpa->toTPA->nX == 1,2,3 + // XXX require fpa->toSky->type == PS_PROJ_TAN + + // convert fpa->toTPA + fpa->toSky to wcs format (DIS) + pmAstromWCS *wcs = pmAstromWCSAlloc(fpa->toTPA->x->nX, fpa->toTPA->x->nY); + + // 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, PS_PROJ_DIS); + wcs->crval1 = fpa->toSky->R*PM_DEG_RAD; + wcs->crval2 = fpa->toSky->D*PM_DEG_RAD; + + // given transformation, solve for coordinates which yields output coordinates of 0,0 + psPlane *center = psPlaneTransformGetCenter (fpa->toTPA, tol); + + // adjust wcs transform to use center as reference coordinate + psPlaneTransformSetCenter (wcs->trans, fpa->toTPA, center->x, center->y); + + // calculated center is crpix1,2 + wcs->crpix1 = center->x; + wcs->crpix2 = center->y; + psFree (center); + + // output coordinates are in FPA pixels + wcs->cdelt1 = 1.0; + wcs->cdelt2 = 1.0; + + return wcs; +} + +static void pmAstromWCSFree (pmAstromWCS *wcs) +{ + + if (!wcs) + return; + psFree (wcs->trans); + psFree (wcs->toSky); +} + +pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder) +{ + + pmAstromWCS *wcs = (pmAstromWCS *) psAlloc(sizeof(pmAstromWCS)); + psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree); + + wcs->trans = psPlaneTransformAlloc (nXorder, nYorder); + wcs->toSky = NULL; + + memset (wcs->ctype1, 0, PM_ASTROM_WCS_TYPE_SIZE); + memset (wcs->ctype2, 0, PM_ASTROM_WCS_TYPE_SIZE); + return wcs; } @@ -900,19 +690,4 @@ *****/ - -// 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 - /* discussion of the coord transformations: X,Y: coord on a chip in pixels Index: /trunk/psModules/src/astrom/pmAstrometryWCS.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10824) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10825) @@ -7,6 +7,6 @@ * @author EAM, IfA * -* @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-17 09:46:56 $ +* @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-22 21:23:06 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -29,11 +29,4 @@ pmAstromWCS; -// read wcs terms from the supplied header into the fpa hierarchy components -bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale, bool isMosaic); - -// write the wcs terms from the fpa hierarchy components into the supplied header -// tol is the convergence tolerance for the non-linear solution to the reference pixel -bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol); - // support function for the pmAstromWCS representation pmAstromWCS *pmAstromWCSAlloc (int nXorder, int nYorder); @@ -45,10 +38,28 @@ pmAstromWCS *pmAstromWCSfromHeader (const psMetadata *header); +// convert from wcs terms to chip->toFPA, fpa->toSky,toTPA terms +bool pmAstromWCSBileveltoChip (pmChip *chip, const pmAstromWCS *wcs); +bool pmAstromWCSBileveltoFPA (pmFPA *fpa, const pmAstromWCS *wcs); + +// convert from chip->toFPA, fpa->toSky,toTPA terms to wcs terms +pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, float tol); +pmAstromWCS *pmAstromWCSBilevelMosaicFromFPA (const pmFPA *fpa, float tol); + // convert the pmAstromWCS representation to the FPA representation -bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale, bool isMosaic); +bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale); pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, float tol); -bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale); -psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale); +// read wcs terms from the supplied header into the fpa hierarchy components +bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale); + +// write the wcs terms from the fpa hierarchy components into the supplied header +// tol is the convergence tolerance for the non-linear solution to the reference pixel +bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol); + +bool pmAstromReadBilevelChip (pmChip *chip, const psMetadata *header); +bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header); + +bool pmAstromWriteBilevelChip (psMetadata *header, const pmChip *chip, float tol); +bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, float tol); // move to pslib