Index: /trunk/ppstamp/src/ppstampCleanup.c =================================================================== --- /trunk/ppstamp/src/ppstampCleanup.c (revision 15317) +++ /trunk/ppstamp/src/ppstampCleanup.c (revision 15318) @@ -20,5 +20,5 @@ // psMemCheckLeaks (0, NULL, stderr, false); - // fprintf (stderr, "Found %d leaks in %s\n", psMemCheckLeaks (0, NULL, NULL, false), "ppstamp"); + fprintf (stderr, "Found %d leaks in %s\n", psMemCheckLeaks (0, NULL, NULL, false), "ppstamp"); return; Index: /trunk/ppstamp/src/ppstampMakeStamp.c =================================================================== --- /trunk/ppstamp/src/ppstampMakeStamp.c (revision 15317) +++ /trunk/ppstamp/src/ppstampMakeStamp.c (revision 15318) @@ -5,11 +5,12 @@ #include "ppstamp.h" #include "pmFPAAstrometry.h" - -#define USE_MOSAIC +#include "pmAstrometryUtils.h" + +#define WCS_NONLIN_TOL 0.001 // Non-linear tolerance for header WCS typedef enum { - PPSTAMP_OFF_CHIP, - PPSTAMP_PARTIALLY_ON_CHIP, - PPSTAMP_ON_CHIP, + PPSTAMP_OFF, + PPSTAMP_PARTIALLY_ON, + PPSTAMP_ON, PPSTAMP_ERROR } ppstampOverlap; @@ -17,37 +18,43 @@ static bool regionContainsPoint(psRegion *region, psPlane *point); static bool regionContainsRegion(psRegion *outer, psRegion *inner); -static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip); - -static pmFPAfile *buildMosaic(pmConfig *config, pmFPAfile *input, pmFPAview *view, pmFPAview *mosaicView) -{ - bool status; - - pmFPAfile *mosaic = psMetadataLookupPtr(&status, config->files, "PPSTAMP.CHIP"); - if (!status) { - psErrorStackPrint(stderr, "can't find mosaic i/o file\n"); - exit(EXIT_FAILURE); - } - - pmChip *mChip = pmFPAviewThisChip(mosaicView, mosaic->fpa); - pmChip *inChip = pmFPAviewThisChip(view, input->fpa); - if (!mChip->hdu && !mChip->parent->hdu) { - const char *name = psMetadataLookupStr(&status, input->fpa->concepts, "FPA.NAME"); // Name of FPA - pmFPAAddSourceFromView(mosaic->fpa, name, mosaicView, mosaic->format); - } - - psMaskType blankMask = pmConfigMask("BLANK", config); - - status = pmChipMosaic(mChip, inChip, true, blankMask); - if (status) { - return mosaic; - } else { - return NULL; - } +static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell, + psRegion *roi); + + +static ppstampOverlap findCell(pmFPAview *view, ppstampOptions *options, pmFPAfile *input, + pmAstromObj *center, pmCell **ppCell) +{ + pmCell *cell; + + view->cell = -1; + while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) { + psRegion *cellExtent = pmCellExtent(cell); + if (regionContainsPoint(cellExtent, center->chip)) { + if (regionContainsRegion(cellExtent, &options->roi)) { + *ppCell = cell; + psFree(cellExtent); + return PPSTAMP_ON; + } else { + fprintf(stderr, "ROI must be confined to single cell. Partial overlap cell %d\n", + view->cell); + psFree(cellExtent); + return PPSTAMP_PARTIALLY_ON; + } + } + psFree(cellExtent); + } + + // This shouldn't ever happen since ROI is on the chip, it must at least partially overlap + // one of the cells + psError(PS_ERR_PROGRAMMING, false, "findCell couldn't find a cell containing the center\n"); + + return PPSTAMP_OFF; } static bool makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input, - pmChip *chip, pmFPAview *view) + pmChip *inChip, pmCell *inCell, pmFPAview *view) { int status = false; + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTAMP.OUTPUT"); if (!output) { @@ -66,23 +73,14 @@ outview = NULL; + // psMetadataPrint(stderr, inChip->concepts, 0); + // psMetadataPrint(stderr, inCell->concepts, 0); // these default to zero would that be ok? psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE, "Binning in x", 1); psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, "Binning in y", 1); - // we extract our postage stamp from a mosaic so that - // pmChipMosaic can handle all of the complexities of parity, binning, etc. - pmFPAview *mosaicView = pmFPAviewAlloc(0); - mosaicView->chip = 0; - pmFPAfile *mosaic = buildMosaic(config, input, view, mosaicView); - if (mosaic == NULL) { - psFree(mosaicView); - return false; - } - - mosaicView->cell = 0; pmReadout *readout; - while ((readout = pmFPAviewNextReadout (mosaicView, mosaic->fpa, 1)) != NULL) { + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { if (!readout->data_exists) { - // HMM, this probably means something is wrong + // XXX if this happens something is wrong continue; } @@ -94,5 +92,5 @@ } - // define a subset of the mosaic's readout + // XXX Do we have to do somethign to deal with negative parity psImage *subsetImage = psImageSubset(readout->image, options->roi); if (subsetImage) { @@ -101,6 +99,6 @@ psFree(subsetImage); if (outReadout->image) { - // TODO: put this into a function in psImage - // (What I really needed was a function to "create a copy of this portion of an image) + // XXX put this into a function in psImage + // (What I really needed was a function to "create a copy of this portion of an image") P_PSIMAGE_SET_COL0(outReadout->image, 0); P_PSIMAGE_SET_ROW0(outReadout->image, 0); @@ -121,13 +119,83 @@ } - psFree(mosaicView); - if (status) { - status = copyMetadata(output, input, chip); + status = copyMetadata(output, input, inChip, inCell, &options->roi); } return status; } -static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip) + +static bool copyWCS(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) + +#ifdef BRUTE_FORCE + outChip->toFPA = inChip->toFPA; + outChip->fromFPA = inChip->fromFPA; + inChip->toFPA = 0; + inChip->fromFPA = 0; + + // This gets us within one pixel of the correct answer + psMetadataItemSupplement(outHDU->header, hdu->header, "CD1_1"); + psMetadataItemSupplement(outHDU->header, hdu->header, "CD1_2"); + psMetadataItemSupplement(outHDU->header, hdu->header, "CD2_1"); + psMetadataItemSupplement(outHDU->header, hdu->header, "CD2_2"); + psMetadataItemSupplement(outHDU->header, hdu->header, "CRVAL1"); + psMetadataItemSupplement(outHDU->header, hdu->header, "CRVAL2"); + + double crpix1 = psMetadataLookupF64(NULL, hdu->header, "CRPIX1"); + double crpix2 = psMetadataLookupF64(NULL, hdu->header, "CRPIX2"); + + crpix1 -= roi->x0; + crpix2 -= roi->y0; + psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", crpix1); + psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", crpix2); + + +#else + + bool combineTransform = false; + if (combineTransform) { + // translation from output chip coordinates to input chip coordinates + psPlaneTransform *trans = psPlaneTransformAlloc(1, 1); + + trans->x->coeff[0][0] = roi->x0; + trans->x->coeff[0][1] = 0; + trans->x->coeff[1][0] = 1.0; + trans->x->coeff[1][1] = 0; + + trans->y->coeff[0][0] = roi->y0; + trans->y->coeff[0][1] = 1.0; + trans->y->coeff[1][0] = 0; + trans->y->coeff[1][1] = 0; + + // combine the translation with the input's transform + // Note: the region and magic number are not actually used by psPlaneTransformCombine + outChip->toFPA = psPlaneTransformCombine(NULL, trans, inChip->toFPA, *roi, 50); + + // for completeness we set fromFPA, but since we don't use it ... + outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50); + + // psFree(region); + psFree(trans); + } else { + outChip->toFPA = psPlaneTransformSetCenter(NULL, inChip->toFPA, roi->x0, roi->y0); + outChip->fromFPA = psPlaneTransformInvert(NULL, outChip->toFPA, *roi, 50); + } + + + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_NONLIN_TOL)) { + psError(PS_ERR_UNKNOWN, false, "Failed to write WCS\n"); + return false; + } +#endif + + return true; +} + + +static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell, psRegion *roi) { pmChip *outChip; @@ -141,5 +209,5 @@ outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts); - pmHDU *inHDU = pmHDUGetHighest(input->fpa, inChip, NULL); + pmHDU *inHDU = pmHDUGetHighest(input->fpa, inChip, NULL); pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL); @@ -149,29 +217,42 @@ outHDU->header = psMetadataAlloc(); } - // TODO set up the astrometry related information + + // 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 references + input->fpa->toTPA = 0; + input->fpa->fromTPA = 0; + input->fpa->toSky = 0; + + if (!copyWCS(outHDU, output->fpa, outChip, inChip, roi)) { + return false; + } + return true; } -ppstampOverlap findROI(ppstampOptions *options, pmFPAfile *input, pmChip *chip, - pmAstromObj *center, psRegion *roi) +// findROI: find whether the region of interest is completely contained in a cell on a given chip. +// If so return the cell +ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view, pmFPAfile *input, pmChip *chip, + pmAstromObj *center, psRegion *roi, pmCell **ppCell) { psRegion *chipExtent = pmChipPixels(chip); bool onChip = false; - ppstampOverlap returnval = PPSTAMP_OFF_CHIP; + ppstampOverlap returnval = PPSTAMP_OFF; psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); - if (options->celestialCenter || options->celestialRange) { - // set up the astrometry - pmHDU *hdu = pmHDUFromChip(chip); - PS_ASSERT_PTR_NON_NULL(hdu, 1) - PS_ASSERT_PTR_NON_NULL(hdu->header, 1) - - if (!pmAstromReadWCS(input->fpa, chip, hdu->header, 1.0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS\n"); - psFree(chipExtent); - psFree(center); - return PPSTAMP_ERROR; - } + // set up the astrometry + pmHDU *hdu = pmHDUFromChip(chip); + PS_ASSERT_PTR_NON_NULL(hdu, 1) + PS_ASSERT_PTR_NON_NULL(hdu->header, 1) + + if (!pmAstromReadWCS(input->fpa, chip, hdu->header, 1.0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS\n"); + psFree(chipExtent); + psFree(center); + return PPSTAMP_ERROR; } @@ -191,5 +272,6 @@ if (regionContainsPoint(chipExtent, center->chip)) { - psLogMsg("ppstampMakeStamp", 3, "Found center on chip: %s\n", chipName); + psLogMsg("ppstampMakeStamp", 3, "Found center (%f %f) on chip: %s\n", + center->chip->x, center->chip->y, chipName); onChip = true; } @@ -227,7 +309,11 @@ if (regionContainsRegion(chipExtent, roi)) { - returnval = PPSTAMP_ON_CHIP; + returnval = findCell(view, options, input, center, ppCell); } else { - returnval = PPSTAMP_PARTIALLY_ON_CHIP; + fprintf(stderr, "Partial Overlap chip %s\n", chipName); + fprintf(stderr, "ChipExtent: %s\n", psRegionToString(*chipExtent)); + fprintf(stderr, "ROI: %s\n", psRegionToString(*roi)); + + returnval = PPSTAMP_PARTIALLY_ON; } } @@ -255,5 +341,4 @@ // files associated with the science image if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - // TODO log error psError(PS_ERR_UNKNOWN, false, "Failed to load input."); psFree(center); @@ -272,5 +357,4 @@ if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - // TODO log error psError(PS_ERR_UNKNOWN, false, "failed to load chip"); status = false; @@ -278,13 +362,14 @@ } - ppstampOverlap overlap = findROI(options, input, chip, center, &options->roi); - - if (overlap == PPSTAMP_ON_CHIP) { - - returnval = makeStamp(config, options, input, chip, view); + pmCell *cell; + ppstampOverlap overlap = findROI(options, view, input, chip, center, &options->roi, &cell); + + if (overlap == PPSTAMP_ON) { + + returnval = makeStamp(config, options, input, chip, cell, view); allDone = true; - } else if (overlap == PPSTAMP_PARTIALLY_ON_CHIP) { - psError(PS_ERR_UNKNOWN, false, "Region of interest must be confined to a single chip\n"); + } else if (overlap == PPSTAMP_PARTIALLY_ON) { + psError(PS_ERR_UNKNOWN, false, "Region of interest must be confined to a single cell\n"); // XXX: perhaps print a helpful message about what coordinates would be valid returnval = false; @@ -301,4 +386,16 @@ view->cell = -1; break; + } else { + // Remove the transformations from the fpa so that they will be + // recalculated from scratch for the next chip. + // If I don't do this the calculation of the astrometry parameters for the output file + // don't come out correctly. See the code on the false side of the test + // if (fpa->toSky == NULL) in the function pmAstromWCStoFPA + psFree(input->fpa->fromTPA); + psFree(input->fpa->toTPA); + psFree(input->fpa->toSky); + input->fpa->fromTPA = NULL; + input->fpa->toTPA = NULL; + input->fpa->toSky = NULL; } }