Index: trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 27460) +++ trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 27461) @@ -865,8 +865,9 @@ } +#define noCOMPARE_TRANS #ifdef COMPARE_TRANS // compare the computed coordintes from this transform with the original psPlane *new = psPlaneAlloc(); - printf(" i chip_x fpa_x fpa_x_fit dx chip_y fpa_y fpa_y_fit dy dx > 0.5 || dy > 0.5\n"); + printf(" i chip_x tpa_x tpa_x_fit dx chip_y tpa_y tpa_y_fit dy dx > 0.5 || dy > 0.5\n"); for (int i=0; iy - d->y; bool bigerr = (fabs(xerr) > .5) || (fabs(yerr) > .5); - printf("%4d %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %s\n" + printf("%4d %9.2f %9.2f %9.2f %9.4f %9.2f %9.2f %9.2f %9.4f %s\n" , i, s->x, new->x, d->x, xerr, s->y, new->y, d->y, yerr, bigerr ? "BIGERR" : ""); } @@ -891,56 +892,175 @@ } -bool pmAstromLinearizeTransforms(pmFPA *fpa, pmChip *chip) -{ - psRegion *chipBounds = pmChipPixels(chip); - -#ifdef TESTING_CAMRUN_14728 - chipBounds->y0 = 1874.; - chipBounds->x1 = 2387.; -#endif - - // First combine the chip to FPA and FPA to TPA into a single transformation - psPlaneTransform *chipToTPA = psPlaneTransformCombine(NULL, chip->toFPA, fpa->toTPA, *chipBounds, 50); +bool pmAstromLinearizeTransforms(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *outputBounds, double offset_x, double offset_y) +{ + PS_ASSERT_PTR_NON_NULL(inFPA, NULL); + PS_ASSERT_PTR_NON_NULL(inChip, NULL); + + if (outFPA == NULL) { + outFPA = inFPA; + } + if (outChip == NULL) { + outChip = inChip; + } + if (outputBounds == NULL) { + outputBounds = pmChipPixels(outChip); + } + + // First combine the "chip to FPA" and "FPA to TPA" into a single transformation + psPlaneTransform *chipToTPA = psPlaneTransformCombine(NULL, inChip->toFPA, inFPA->toTPA, *outputBounds, 50); if (!chipToTPA) { - psFree(chipBounds); psError(PS_ERR_UNKNOWN, false, "failed to create chipToTPA"); return false; } - // Next do a linear fit - psPlaneTransform *newToFPA = linearFitToTransform(chipToTPA, chipBounds); + // Next do the linear fit within the output boundary pixels + psPlaneTransform *chipToFPA = linearFitToTransform(chipToTPA, outputBounds); psFree(chipToTPA); - if (!newToFPA) { - psFree(chipBounds); + if (!chipToFPA) { psError(PS_ERR_UNKNOWN, false, "linear fit of chip to TPA transform failed"); return false; } - psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *chipBounds, 50); - psFree(chipBounds); + // if requested, change the center + psPlaneTransform *outToFPA; + if (offset_x != 0. && offset_y != 0.) { + outToFPA = psPlaneTransformSetCenter(NULL, chipToFPA, offset_x, offset_y); + psFree(chipToFPA); + } else { + outToFPA = chipToFPA; + } + + psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 50); + if (!outFromFPA) { + psFree(outToFPA); + psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed"); + return false; + } + + // Success. Now set the fpa's toTPA and fromTPA to identity and replace the chip's transforms. + + psFree(outFPA->toTPA); + outFPA->toTPA = psPlaneTransformIdentity(1); + + psFree(outFPA->fromTPA); + outFPA->fromTPA = psPlaneTransformIdentity(1); + + psFree(outChip->toFPA); + outChip->toFPA = outToFPA; + + psFree(outChip->fromFPA); + outChip->fromFPA = outFromFPA; + + // Finally, change the type for the projection. + outFPA->toSky->type = PS_PROJ_TAN; + + return true; +} + +bool pmAstromLinearizeToSky(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *bounds) +{ + PS_ASSERT_PTR_NON_NULL(inFPA, NULL); + PS_ASSERT_PTR_NON_NULL(inChip, NULL); + PS_ASSERT_PTR_NON_NULL(outFPA, NULL); + PS_ASSERT_PTR_NON_NULL(outChip, NULL); + PS_ASSERT_PTR_NON_NULL(bounds, NULL); + + // outFPA projection must be defined as the goal + + // the output transformations are: + // chip -> FPA : standard linear trans with needed rotation, etc + // FPA -> TPA : identidy + + int nSamples = 10; // 10 samples in each dimension + + double deltaX = (bounds->x1 - bounds->x0); + double deltaY = (bounds->y1 - bounds->y0); + + psArray *src = psArrayAllocEmpty(nSamples * nSamples); + psArray *dst = psArrayAllocEmpty(nSamples * nSamples); + + psPlane srcFP, srcTP; + + for (int j = 0; j < nSamples; j++) { + double y = bounds->y0 + (j * deltaY / nSamples); + for (int i = 0; i < nSamples; i++) { + + psSphere srcSky; + psPlane *srcChip = psPlaneAlloc(); + psPlane *dstTP = psPlaneAlloc(); + + srcChip->x = bounds->x0 + (i * deltaX / nSamples); + srcChip->y = y; + + psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip); + psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP); + psDeproject (&srcSky, &srcTP, inFPA->toSky); + + // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD); + + psProject (dstTP, &srcSky, outFPA->toSky); + + srcChip->x -= bounds->x0; + srcChip->y -= bounds->y0; + psArrayAdd (src, 100, srcChip); + psArrayAdd (dst, 100, dstTP); + + psFree(srcChip); // drop our refs to s and d + psFree(dstTP); + } + } + + psPlaneTransform *newToFPA = psPlaneTransformAlloc(1, 1); + newToFPA->x->coeffMask[1][1] = 1; + newToFPA->y->coeffMask[1][1] = 1; + + if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) { + psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed"); + psFree(src); + psFree(dst); + return NULL; + } + +# if (0) + for (int i = 0; i < src->n; i++) { + + psSphere srcSky, dstSky; + psPlane *srcChip = src->data[i]; + psPlane *dstTP = dst->data[i]; + + psPlaneTransformApply (&srcFP, newToFPA, srcChip); + psDeproject (&srcSky, &srcFP, outFPA->toSky); + psDeproject (&dstSky, dstTP, outFPA->toSky); + + double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0; + double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0; + fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD); + + } +# endif + + psFree(src); + psFree(dst); + + // this is a linear transformation + psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *bounds, 1); if (!newFromFPA) { psFree(newToFPA); - psError(PS_ERR_UNKNOWN, false, "inversion of fit of chip to TPA transform failed"); + psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed"); return false; } // Success. Now set the fpa's toTPA and fromTPA to identity and replace the chip's transforms. - - psPlaneTransform *newToTPA = psPlaneTransformIdentity(1); - psFree(fpa->toTPA); - fpa->toTPA = newToTPA; - - psPlaneTransform *newFromTPA = psPlaneTransformIdentity(1); - psFree(fpa->fromTPA); - fpa->fromTPA = newFromTPA; - - psFree(chip->toFPA); - chip->toFPA = newToFPA; - - psFree(chip->fromFPA); - chip->fromFPA = newFromFPA; - - // Finally change the type for the projection. - fpa->toSky->type = PS_PROJ_TAN; + psFree(outChip->toFPA); + outChip->toFPA = newToFPA; + + psFree(outChip->fromFPA); + outChip->fromFPA = newFromFPA; + + psFree(outFPA->toTPA); + outFPA->toTPA = psPlaneTransformIdentity(1); + + psFree(outFPA->fromTPA); + outFPA->fromTPA = psPlaneTransformIdentity(1); return true; Index: trunk/psModules/src/astrom/pmAstrometryWCS.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 27460) +++ trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 27461) @@ -63,5 +63,6 @@ bool pmAstromWriteBilevelMosaic (psMetadata *header, const pmFPA *fpa, double tol); -bool pmAstromLinearizeTransforms(pmFPA *fpa, pmChip *); +bool pmAstromLinearizeTransforms(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *outputBounds, double offset_x, double offset_y); +bool pmAstromLinearizeToSky(pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *bounds); // move to pslib Index: trunk/pstamp/src/ppstampMakeStamp.c =================================================================== --- trunk/pstamp/src/ppstampMakeStamp.c (revision 27460) +++ trunk/pstamp/src/ppstampMakeStamp.c (revision 27461) @@ -17,26 +17,58 @@ static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip); +static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip); static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance); // convert the input chip's transforms to the output -static bool convertWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi) -{ - pmHDU *hdu = pmHDUFromChip(inChip); - PS_ASSERT_PTR_NON_NULL(hdu, 1) - PS_ASSERT_PTR_NON_NULL(hdu->header, 1) - - // Change the reference pixel to account for the change in origin between the stamp and - // the original image - // XXX: shouldn't we be using the center of the stamp instead of the corner? - outChip->toFPA = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0); - - outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50); +static bool convertWCS(pmHDU *outHDU, pmFPA *inFPA, pmChip *inChip, pmFPA *outFPA, pmChip *outChip, psRegion *roi) +{ + PS_ASSERT_PTR_NON_NULL(outHDU, 1); + PS_ASSERT_PTR_NON_NULL(outHDU->header, 1); + + // Center of stamp in chip coordinates + double center_x = (roi->x0 + roi->x1) / 2; + double center_y = (roi->y0 + roi->y1) / 2; + pmAstromObj *check = pmAstromObjAlloc(); + check->chip->x = center_x; + check->chip->y = center_y; + chipToSky(check, inFPA, inChip); + double r_before = check->sky->r; + double d_before = check->sky->d; + psLogMsg("ppstampMakeStamp", 2, "Before fit Center Pixel RA: %f DEC: %f (degrees)\n", + r_before * PS_DEG_RAD, d_before * PS_DEG_RAD); if (outFPA->toSky->type != PS_PROJ_TAN) { - if (!pmAstromLinearizeTransforms(outFPA, outChip)) { + // we need to make astrometry terms for the output which are centered on the output chip center, + // but we keep the original plate scale + psFree(outFPA->toSky); + outFPA->toSky = psProjectionAlloc (r_before, d_before, inFPA->toSky->Xs, inFPA->toSky->Ys, PS_PROJ_TAN); + + if (!pmAstromLinearizeToSky(inFPA, inChip, outFPA, outChip, roi)) { + psFree(check); psError(PS_ERR_UNKNOWN, false, "Failed to linearize astrometry\n"); return false; } - } + check->chip->x = (roi->x1 - roi->x0) / 2; + check->chip->y = (roi->y1 - roi->y0) / 2; + chipToSky(check, outFPA, outChip); + double r_after = check->sky->r; + double d_after = check->sky->d; + + psLogMsg("ppstampMakeStamp", 2, "After fit: Center Pixel RA: %f DEC: %f (degrees)\n", + r_after * PS_DEG_RAD, d_after * PS_DEG_RAD); + psLogMsg("ppstampMakeStamp", 2, "Error in fit to astrometry %.2f %.2f (arcseconds)\n", + (r_after - r_before) * PS_DEG_RAD * 3600, (d_after - d_before) * PS_DEG_RAD * 3600); + + // XXX: should we fail if the fit is bad ?? + + } else { + outChip->toFPA = psPlaneTransformSetCenter(NULL, inChip->toFPA, (int) roi->x0, (int) roi->y0); + outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50); + if (!outChip->fromFPA) { + psError(PS_ERR_UNKNOWN, false, "inversion of toFPA transform failed"); + return false; + } + } + psFree(check); if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) { @@ -92,15 +124,12 @@ // If input had WCS convert it for stamp if (input->fpa->toSky) { - // steal the input fpa's transforms - output->fpa->toTPA = input->fpa->toTPA; - output->fpa->fromTPA = input->fpa->fromTPA; - output->fpa->toSky = input->fpa->toSky; - - // drop the references - input->fpa->toTPA = 0; - input->fpa->fromTPA = 0; - input->fpa->toSky = 0; - - if (!convertWCS(outHDU, output->fpa, outChip, inChip, &options->roi)) { + // copy the input fpa's transforms + // XXX: if we change to making multiple stamps per process we'll need to copy + // the transformations + output->fpa->toTPA = psMemIncrRefCounter(input->fpa->toTPA); + output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA); + output->fpa->toSky = psMemIncrRefCounter(input->fpa->toSky); + + if (!convertWCS(outHDU, input->fpa, inChip, output->fpa, outChip, &options->roi)) { return false; }