Index: /trunk/psastro/src/psastroAstromGuess.c =================================================================== --- /trunk/psastro/src/psastroAstromGuess.c (revision 12535) +++ /trunk/psastro/src/psastroAstromGuess.c (revision 12536) @@ -47,4 +47,17 @@ pmFPA *fpa = input->fpa; + // load mosaic-level astrometry? + bool bilevelAstrometry = false; + pmHDU *phu = pmFPAviewThisPHU (view, fpa); + if (phu) { + char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } + } + if (bilevelAstrometry) { + pmAstromReadBilevelMosaic (fpa, phu->header); + } + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); @@ -53,6 +66,10 @@ // read WCS data from the corresponding header pmHDU *hdu = pmFPAviewThisHDU (view, fpa); + if (bilevelAstrometry) { + pmAstromReadBilevelChip (chip, hdu->header); + } else { + pmAstromReadWCS (fpa, chip, hdu->header, pixelScale); + } - pmAstromReadWCS (fpa, chip, hdu->header, pixelScale); if (newFPA) { newFPA = false; Index: /trunk/psastro/src/psastroLoadRefstars.c =================================================================== --- /trunk/psastro/src/psastroLoadRefstars.c (revision 12535) +++ /trunk/psastro/src/psastroLoadRefstars.c (revision 12536) @@ -83,4 +83,10 @@ unlink (tempFile); + + if (table == NULL) { + psError(PSASTRO_ERR_REFSTARS, true, "failure to load astrometric reference\n"); + return NULL; + } + psLogMsg ("psastro", 3, "read getstar output table : %f sec\n", psTimerMark ("psastro")); Index: /trunk/psastro/src/psastroMosaicAstrom.c =================================================================== --- /trunk/psastro/src/psastroMosaicAstrom.c (revision 12535) +++ /trunk/psastro/src/psastroMosaicAstrom.c (revision 12536) @@ -4,6 +4,4 @@ // XXX require this fpa to have multiple chip extensions and a PHU? bool psastroMosaicAstrom (pmConfig *config, psArray *refs) { - - bool status; // select the current recipe @@ -13,11 +11,4 @@ return false; } - - // physical pixel scale in microns per pixel - double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); - if (!status) { - psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); - return false; - } // select the input data sources @@ -101,22 +92,6 @@ if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.6.dat"); } - // save WCS and analysis metadata in update header XXX this is still wrong: the pmFPAExtent - // function returns the FP dimensions in pixel units relative to the 0,0 corner of chip - // 0,0. I need to have a function which loops over all cells, determines the pixel - // dimensions for each cell, converts them to chip pixels, then uses the fpa astrometry - // structures to get the total dimensions of the fpa in fpa units. equivalent to - // pmFPAExtent - psRegion *region = pmFPAExtent (fpa); - region->x0 *= pixelScale; - region->x1 *= pixelScale; - region->y0 *= pixelScale; - region->y1 *= pixelScale; - - // XXX also need to add DATE/TIME info and NAXIS1, NAXIS2 + // save WCS and analysis metadata in update header. psMetadata *updates = psMetadataAlloc(); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "fpa dimensions (mm)", region->x1 - region->x0); - psMetadataAddS32 (updates, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "fpa dimensions (mm)", region->y1 - region->y0); - psFree (region); - if (!pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL)) { psAbort ("failed to save header terms"); Index: unk/psastro/src/psastroMosaicRescaleChips.c =================================================================== --- /trunk/psastro/src/psastroMosaicRescaleChips.c (revision 12535) +++ (revision ) @@ -1,28 +1,0 @@ -# include "psastro.h" - -// XXX this file is no longer used - -// XXX what is this doing? -// shouldn't it be doing toFPA -> fromFPA - -bool psastroMosaicRescaleChips (pmFPA *fpa) { - - pmChip *chip = NULL; - pmFPAview *view = pmFPAviewAlloc (0); - - // this loop selects the matched stars for all chips - - // psRegion region = psMetadataLookupXXX (chip->concepts, "CHIP.TRIMSEC"); - // psRegion region = psRegionSet (0, 4000, 0, 4000); - // chip->fromFPA = psPlaneTransformInvert(NULL, toFPA, region, 50); - - while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { - psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); - if (!chip->process || !chip->file_exists) { continue; } - - // XXX make '50' be a fraction of the chip size? - chip->fromFPA = psPlaneTransformInvert(chip->fromFPA, chip->toFPA, region, 50); - } - psFree (view); - return true; -} Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 12535) +++ /trunk/psastro/src/psastroUtils.c (revision 12536) @@ -79,11 +79,8 @@ bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) { - psRegion *region = pmChipExtent (chip); - region->x1 -= region->x0; - region->y1 -= region->y0; - region->x0 = 0; - region->y0 = 0; + psRegion *region = pmChipPixels (chip); + psFree (chip->fromFPA); - chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, *region, 20); + chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, *region, 50); psFree (region); @@ -134,148 +131,4 @@ } return true; -} - -bool psastroNormFPA (pmFPA *fpa, psMetadata *config) { - - // save the raw astrometry for later reference - pmFPA *raw = pmFPACopyAstrom (fpa); - - // first pass: measure the per-chip solutions, modify the chip.toFPA terms - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - psastroChipAstrom (chip, config); - } - - // second stage: re-normalize the chip terms, passing the - // average rotation and offset values to the fpa.toSky - if (RENORM) { - - // this code is needed for the mosastro stage, with multiple chip solutions - - double dX, dY, dT, dN; - dX = dY = dT = dN = 0; - - psPlane origin, P1, P2; - origin.x = 0; - origin.y = 0; - - // calculate the average rotation and boresite offset relative to raw - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *iChip = raw->chips->data[i]; - pmChip *oChip = fpa->chips->data[i]; - - // offset of chip - psCoordChipToFPA (&P1, &origin, iChip); - psCoordChipToFPA (&P2, &origin, oChip); - dX += (P2.x - P1.x); - dY += (P2.y - P1.y); - - // get parity-independent rotations for old and new solutions - double T1 = psPlaneTransformGetRotation (iChip->toFPA); - double T2 = psPlaneTransformGetRotation (oChip->toFPA); - dT += T2 - T1; - dN ++; - } - - dT /= dN; - dX /= dN; - dY /= dN; - - // R(T) - double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0]; - double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0]; - double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0]; - double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0]; - - // R(dT) - double dPC1_1 = +cos (dT); - double dPC1_2 = +sin (dT); - double dPC2_1 = -sin (dT); - double dPC2_2 = +cos (dT); - - // R'(T) = R(T) * R(dT) - double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1; - double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2; - double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1; - double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2; - - double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1); - - // R'(-T) (matrix inverse, not just rotation inverse -- keeps parity) - double pi1_1 = +pc2_2 * det; - double pi1_2 = -pc1_2 * det; - double pi2_1 = -pc2_1 * det; - double pi2_2 = +pc1_1 * det; - - // apply the new modifcations in rotation and boresite - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *oChip = fpa->chips->data[i]; - - // r(T) - double pr1_1 = oChip->toFPA->x->coeff[1][0]; - double pr1_2 = oChip->toFPA->x->coeff[0][1]; - double pr2_1 = oChip->toFPA->y->coeff[1][0]; - double pr2_2 = oChip->toFPA->y->coeff[0][1]; - - // ri'(T) = R(T) r(t) - double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1; - double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2; - double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1; - double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2; - - // r'(T) = R'(-T) ri'(T) - oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1; - oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2; - oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1; - oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2; - - double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX; - double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY; - - oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy; - oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy; - } - - fpa->toTPA->x->coeff[0][0][0][0] -= dX; - fpa->toTPA->y->coeff[0][0][0][0] -= dY; - - fpa->toTPA->x->coeff[1][0][0][0] = pc1_1; - fpa->toTPA->x->coeff[0][1][0][0] = pc1_2; - fpa->toTPA->y->coeff[1][0][0][0] = pc2_1; - fpa->toTPA->y->coeff[0][1][0][0] = pc2_2; - } - return true; -} - -psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input) { - - psPolynomial2D *output = psPolynomial2DAlloc (input->nX, input->nY, input->type); - - for (int i = 0; i < input->nX; i++) { - for (int j = 0; j < input->nY; j++) { - output->mask[i][j] = input->mask[i][j]; - output->coeff[i][j] = input->coeff[i][j]; - output->coeffErr[i][j] = input->coeffErr[i][j]; - } - } - return (output); -} - -psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input) { - - psPolynomial4D *output = psPolynomial4DAlloc (input->nX, input->nY, input->nZ, input->nT, input->type); - - for (int i = 0; i < input->nX; i++) { - for (int j = 0; j < input->nY; j++) { - for (int k = 0; k < input->nZ; k++) { - for (int m = 0; m < input->nT; m++) { - output->mask[i][j][k][m] = input->mask[i][j][k][m]; - output->coeff[i][j][k][m] = input->coeff[i][j][k][m]; - output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m]; - } - } - } - } - return (output); } @@ -328,34 +181,4 @@ } -// very crude distortion inversion: assumes 0 order in z and t, linear in x and y: -psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort) { - PS_ASSERT_PTR_NON_NULL(distort, 0); - PS_ASSERT_PTR_NON_NULL(distort->x, 0); - PS_ASSERT_PTR_NON_NULL(distort->y, 0); - - psPlaneDistort *out = psPlaneDistortAlloc(1, 1, 0, 0); - - /* simple matrix inversion code */ - - psF64 r11 = distort->x->coeff[1][0][0][0]; - psF64 r12 = distort->x->coeff[0][1][0][0]; - psF64 r21 = distort->y->coeff[1][0][0][0]; - psF64 r22 = distort->y->coeff[0][1][0][0]; - psF64 xo = distort->x->coeff[0][0][0][0]; - psF64 yo = distort->y->coeff[0][0][0][0]; - - psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); // Inverse of the determinant - - out->x->coeff[1][0][0][0] = +invDet * r22; - out->x->coeff[0][1][0][0] = -invDet * r12; - out->y->coeff[1][0][0][0] = -invDet * r21; - out->y->coeff[0][1][0][0] = +invDet * r11; - - out->x->coeff[0][0][0][0] = - invDet * (r22 * xo - r12 * yo); - out->y->coeff[0][0][0][0] = - invDet * (r11 * yo - r21 * xo); - - return(out); -} - // returns the rotation term, forcing positive parity double psPlaneTransformGetRotation (psPlaneTransform *map) { Index: unk/psastro/src/psastroWCS.c =================================================================== --- /trunk/psastro/src/psastroWCS.c (revision 12535) +++ (revision ) @@ -1,577 +1,0 @@ -# include "psastro.h" - -// interpret header WCS (only handles traditional WCS for the moment) -// plateScale is nominal physical scale on tangent plane (microns / arcsecond) -bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) { - - psProjectionType type; - bool status, pcKeys, cdKeys; - float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2; - float pc1_1, pc1_2, pc2_1, pc2_2; - - // interpret header data, convert to crval(i), etc - char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2"); - if (!status) { - psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n"); - return false; - } - - // 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; - if (type == PS_PROJ_NTYPE) { - psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype); - return false; - } - - // what type of WCS keywords are available? - psMetadataLookupF32 (&pcKeys, header, "PC001001"); - psMetadataLookupF32 (&cdKeys, header, "CD1_1"); - - if (cdKeys && pcKeys) { - psLogMsg ("psastro", 2, "warning: both CDi_j and PC00i00j defined in headers, using CDi_j terms\n"); - } - if (!cdKeys && !pcKeys) { - psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms"); - // XXX we could default here to RA, DEC, ROTANGLE - return false; - } - - crval1 = psMetadataLookupF32 (&status, header, "CRVAL1"); - crval2 = psMetadataLookupF32 (&status, header, "CRVAL2"); - crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1"); - crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2"); - - // test the CDELTi varient - if (pcKeys) { - cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1"); - cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2"); - - // test the CROTAi varient: - double rotate = psMetadataLookupF32 (&status, header, "CROTA2"); - if (status) { - double Lambda = cdelt2 / cdelt1; - pc1_1 = cos(rotate*RAD_DEG); - pc1_2 = -sin(rotate*RAD_DEG) * Lambda; - pc2_1 = sin(rotate*RAD_DEG) / Lambda; - pc2_2 = cos(rotate*RAD_DEG); - goto got_matrix; - } - - // test the PC00i00j varient: - pc1_1 = psMetadataLookupF32 (&status, header, "PC001001"); - if (status) { - pc1_2 = psMetadataLookupF32 (&status, header, "PC001002"); - pc2_1 = psMetadataLookupF32 (&status, header, "PC002001"); - pc2_2 = psMetadataLookupF32 (&status, header, "PC002002"); - - // XXX EAM : add Elixir polynomial terms here eventually - goto got_matrix; - } - psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); - return false; - } - - // test the CDi_j varient - if (cdKeys) { - pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1"); - pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2"); - pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1"); - pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2"); - - // renormalize to cdelt1, cdelt2, etc - double scale = hypot (pc1_1, pc1_2); - cdelt1 = cdelt2 = scale; - pc1_1 /= scale; - pc1_2 /= scale; - pc2_1 /= scale; - pc2_2 /= scale; - goto got_matrix; - } - psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); - return false; - -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 - - *****/ - - { - double rX = 1.0; - double rY = 1.0; - - // XXX free an existing toFPA? - psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1); - - // basic transformation from chip to FPA - 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; - - // projection from TPA to SKY - psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type); - - if (fpa->toSky == NULL) { - // XXX for now, use the identity for TPA <--> FPA - fpa->toTPA = psPlaneDistortIdentity (1); - fpa->fromTPA = psPlaneDistortIdentity (1); - fpa->toSky = toSky; - } else { - if (fpa->toTPA == NULL) psAbort("projection defined, tangent-plane not defined"); - if (fpa->fromTPA == NULL) psAbort("projection defined, tangent-plane not defined"); - - // adjust for common toSky for mosaic: - // ignore the TPA since toTPA is identity - // find the FPA coordinate of 0,0 for this chip. - psPlane *fp = psPlaneAlloc(); - psPlane *chip = 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 - psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate - psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection - - toFPA->x->coeff[0][0] = fp->x; - toFPA->y->coeff[0][0] = fp->y; - - rX = toSky->Xs/fpa->toSky->Xs; - rY = toSky->Ys/fpa->toSky->Ys; - - // correct to common plate scale - toFPA->x->coeff[1][0] *= rX; - toFPA->x->coeff[0][1] *= rX; - toFPA->y->coeff[1][0] *= rY; - toFPA->y->coeff[0][1] *= rY; - - psFree (fp); - psFree (sky); - psFree (chip); - psFree (toSky); - } - - chip->toFPA = toFPA; - chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA); - - 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 - if (isMosaic) { - chip->toFPA->x->coeff[1][0] /= rX; - chip->toFPA->x->coeff[0][1] /= rX; - 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", - DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D, - 3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys); - } - 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("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) { - 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", 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 = DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0]; - cdelt2 = 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); - - 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 - for (int i = 0; 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 (toFPA->x, Xo, Yo); - Beta->data.F32[1] = psPolynomial2DEval (toFPA->y, Xo, Yo); - - psMatrixGJSolveF32 (Alpha, Beta); - - Xo += Beta->data.F32[0]; - 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*DEG_RAD); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*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*DEG_RAD); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*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*DEG_RAD*plateScale); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*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 - - return true; -} - -// convert toFPA / toSky components to traditional WCS -// plateScale is nominal physical scale on tangent plane (microns / arcsecond) -bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale) { - - 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("mis-matched tangent plane orders (1)"); - if (toFPA->x->nX != toFPA->y->nX) psAbort("mis-matched tangent plane orders (2)"); - if (toFPA->x->nX != toFPA->y->nY) psAbort("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("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*DEG_RAD); - psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*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*DEG_RAD*plateScale); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale); - - if (toTP->x->nX != toTP->x->nY) psAbort("mis-matched tangent plane orders (1)"); - if (toTP->x->nX != toTP->y->nX) psAbort("mis-matched tangent plane orders (2)"); - if (toTP->x->nX != toTP->y->nY) psAbort("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("invalid tangent plane order"); - } - return header; -} - -// construct a psPlaneDistort which is the identify transformation -psPlaneDistort *psPlaneDistortIdentity (int order) { - - psPlaneDistort *distort; - - if (order < 1) psAbort("invalid order"); - if (order > 3) psAbort("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 -bool psPlaneDistortIsIdentity (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 { - if (i + j != 1) { - // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) - status &= (distort->x->coeff[i][j][0][0] == 0.0); - } else { - // linear, diagonal terms must be 1.0 - status &= (distort->x->coeff[i][j][0][0] == 1.0); - } - } - } - } - - 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 { - if (i + j != 1) { - // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) - status &= (distort->y->coeff[i][j][0][0] == 0.0); - } else { - // linear, diagonal terms must be 1.0 - status &= (distort->y->coeff[i][j][0][0] == 1.0); - } - } - } - } - return status; -}