Index: /trunk/ppstamp/src/Makefile.am =================================================================== --- /trunk/ppstamp/src/Makefile.am (revision 15322) +++ /trunk/ppstamp/src/Makefile.am (revision 15323) @@ -17,7 +17,8 @@ ppstampCleanup.c \ ppstampErrorCodes.c \ + ppstampMakeStamp.c \ ppstampOptions.c \ ppstampParseCamera.c \ - ppstampMakeStamp.c \ + ppstampRegion.c \ ppstampVersion.c Index: /trunk/ppstamp/src/ppstamp.h =================================================================== --- /trunk/ppstamp/src/ppstamp.h (revision 15322) +++ /trunk/ppstamp/src/ppstamp.h (revision 15323) @@ -16,5 +16,4 @@ #include "ppStats.h" -#define RECIPE_NAME "PPSTAMP" // Name of the recipe to use #define TIMERNAME "ppstamp" // Name of timer @@ -29,23 +28,20 @@ typedef struct { // input arguments - double centerX; - double centerY; - double dX; - double dY; + int centerX; + int centerY; + int dX; + int dY; double centerRA; double centerDEC; double dRA; double dDEC; - bool celestialCenter; // true if center is in RA/dec - bool celestialRange; // true if range is in RA/dec + bool celestialCenter; // true if center is in RA/dec + bool celestialRange; // true if range is in RA/dec psString chipName; // - // Calculated Parameters + // Calculated Values // - int chip; - int cell; psRegion roi; // roi in chip coordinates - // psPlane fpaCenter; // center of ROI in FPA coordinates - // ppstampROI fpaROI; // region of interest in FPA coordinates + } ppstampOptions; @@ -64,19 +60,6 @@ bool ppstampMakeStamp(pmConfig *config, ppstampOptions *); -// Loop over the input -bool ppstampLoop(pmConfig *config, ppstampOptions *options); - -bool ppstampReadout(pmConfig *config, pmFPAview *view); - // free memory, check for leaks void ppstampCleanup (pmConfig *config, ppstampOptions *options); - -#ifdef notyet -// calculate stats, including MD5 -bool ppstampStats (pmConfig *config, pmChip *chip, const pmFPAview *inputView); - -// write stats to output file -bool ppstampStatsOutput (pmConfig *config); -#endif /// Return short version information @@ -90,12 +73,6 @@ ); - -// calculate stats, including MD5 -bool ppstampPixelStats (pmConfig *config, const pmFPAview *inputView); - -// calculate stats from headers and concepts -bool ppstampMetadataStats (pmConfig *config); - -void ppstampFileCheck (pmConfig *config); +psRegion *ppstampCellRegion(const pmCell *cell); +psRegion *ppstampChipRegion(const pmChip *chip); #endif Index: /trunk/ppstamp/src/ppstampArguments.c =================================================================== --- /trunk/ppstamp/src/ppstampArguments.c (revision 15322) +++ /trunk/ppstamp/src/ppstampArguments.c (revision 15323) @@ -4,32 +4,31 @@ #include "ppstamp.h" - -static void usage (void) { - fprintf(stderr, "USAGE: ppstamp -pixcenter x y -pixrange dx dy -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); - - fprintf(stderr, "USAGE: ppstamp -pixcenter x y -celrange dRA dDEC -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); +#include "ohana.h" + +static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothDegrees, double *p1, double *p2); +static bool get2Ints(int argnum, int *pArgc, char **argv, int *p1, int *p2); + +// XXX This usage output is kind of busy + +static void usage (void) +{ + fprintf(stderr, "\n"); + + fprintf(stderr, "USAGE: ppstamp -celcenter RA DEC -arcrange dRA dDEC [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); fprintf(stderr, " ppstamp -celcenter RA DEC -celrange dRA dDEC [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); fprintf(stderr, " ppstamp -celcenter RA DEC -pixrange dx dy [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); - fprintf(stderr, "\n"); - -#ifdef notyet + fprintf(stderr, " ppstamp -pixcenter x y -pixrange dx dy -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); + fprintf(stderr, " ppstamp -pixcenter x y -arcrange dRA dDEC -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); + fprintf(stderr, " ppstamp -pixcenter x y -celrange dRA dDEC -chip chipname [-file INPUT.fits] [-list INPUT.txt] OUTPUT\n"); + fprintf(stderr, "\n"); + + fprintf(stderr, "celcenter & celrange may be specified in HH:MM:SS DD:MM:SS or as 2 values in decimal degrees\n"); + fprintf(stderr, "arcrange may be specfied as 2 values DD:MM:SS DD:MM:SS or decimal degrees\n"); + fprintf(stderr, "\n"); fprintf(stderr, "Optional arguments:\n"); - fprintf(stderr, "\t-stats STATS.mdc: Output statistics into STATS.mdc\n"); - fprintf(stderr, "\n"); - fprintf(stderr, "Input options (single file / file list):\n"); - fprintf(stderr, "\n"); -#endif + fprintf(stderr, " -chip chipName selects chip (only used with -pixcenter)\n"); + fprintf(stderr, "\n"); + exit (2); -} - -void get2Doubles(int argnum, int *pArgc, char **argv, double *p1, double *p2) { - if (*pArgc < 3) { - usage(); - } - psArgumentRemove(argnum, pArgc, argv); - *p1 = atof(argv[argnum]); - psArgumentRemove(argnum, pArgc, argv); - *p2 = atof(argv[argnum]); - psArgumentRemove(argnum, pArgc, argv); } @@ -48,8 +47,4 @@ psString version; version = ppstampVersionLong(); fprintf (stdout, "%s\n", version); psFree (version); -// don't need these -// version = ppStatsVersionLong(); fprintf (stdout, "%s\n", version); psFree (version); -// version = psphotVersionLong(); fprintf (stdout, "%s\n", version); psFree (version); -// version = psastroVersionLong(); fprintf (stdout, "%s\n", version); psFree (version); version = psModulesVersionLong(); fprintf (stdout, "%s\n", version); psFree (version); version = psLibVersionLong(); fprintf (stdout, "%s\n", version); psFree (version); @@ -58,5 +53,5 @@ // load the site-wide configuration information - pmConfig *config = pmConfigRead(&argc, argv, RECIPE_NAME); + pmConfig *config = pmConfigRead(&argc, argv, NULL); if (config == NULL) { psErrorStackPrint(stderr, "Can't find site configuration!\n"); @@ -64,6 +59,4 @@ } - (void) pmConfigRecipeOptions (config, RECIPE_NAME); - options = ppstampOptionsAlloc(); *pOptions = options; @@ -72,6 +65,12 @@ gotCenter = true; options->celestialCenter = false; - get2Doubles(argnum, &argc, argv, &options->centerX, &options->centerY); - } + psArgumentRemove(argnum, &argc, argv); + + if (!get2Ints(argnum, &argc, argv, &options->centerX, &options->centerY)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "invalid pixcenter specification\n"); + usage(); + } + } + if ((argnum = psArgumentGet(argc, argv, "-celcenter"))) { if (gotCenter) { @@ -79,32 +78,65 @@ usage(); } - get2Doubles(argnum, &argc, argv, &options->centerRA, &options->centerDEC); - options->centerRA = DEG_TO_RAD(options->centerRA); - options->centerDEC = DEG_TO_RAD(options->centerDEC); + psArgumentRemove(argnum, &argc, argv); + + double raDeg, decDeg; + if (!get2Angles(argnum, &argc, argv, false, &raDeg, &decDeg)) { + usage(); + } + options->centerRA = DEG_TO_RAD(raDeg); + options->centerDEC = DEG_TO_RAD(decDeg); gotCenter = true; options->celestialCenter = true; } + if (!gotCenter) { - usage(); - } + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify center of region of interest\n"); + usage(); + } + if ((argnum = psArgumentGet(argc, argv, "-pixrange"))) { gotRange = true; options->celestialRange = false; - get2Doubles(argnum, &argc, argv, &options->dX, &options->dY); - } - if ((argnum = psArgumentGet(argc, argv, "-celrange"))) { - fprintf(stderr, "Specifying range in celestial coordinates is not implemented yet\n"); - exit(1); + psArgumentRemove(argnum, &argc, argv); + if (!get2Ints(argnum, &argc, argv, &options->dX, &options->dY)) { + usage(); + } + } + + if ((argnum = psArgumentGet(argc, argv, "-arcrange"))) { if (gotRange) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "can't specify both -pixrange and -celrange\n");; - usage(); - } + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "specify only one of -pixrange or -arcrange or -celrange\n");; + usage(); + } + psArgumentRemove(argnum, &argc, argv); gotRange = true; options->celestialRange = true; - get2Doubles(argnum, &argc, argv, &options->dRA, &options->dDEC); - options->dRA = DEG_TO_RAD(options->dRA); - options->dDEC = DEG_TO_RAD(options->dDEC); - } + + double deg1, deg2; + if (!get2Angles(argnum, &argc, argv, true, °1, °2)) { + usage(); + } + options->dRA = DEG_TO_RAD(deg1); + options->dDEC = DEG_TO_RAD(deg2); + } + if ((argnum = psArgumentGet(argc, argv, "-celrange"))) { + if (gotRange) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "specify one of -pixrange or -arcrange or -celrange\n");; + usage(); + } + psArgumentRemove(argnum, &argc, argv); + gotRange = true; + options->celestialRange = true; + + double deg1, deg2; + if (!get2Angles(argnum, &argc, argv, false, °1, °2)) { + usage(); + } + options->dRA = DEG_TO_RAD(deg1); + options->dDEC = DEG_TO_RAD(deg2); + } + if (!gotRange) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify range\n"); usage(); } @@ -114,9 +146,4 @@ options->chipName = psStringCopy(argv[argnum]); psArgumentRemove(argnum, &argc, argv); - } - - if (!options->celestialCenter && (options->chipName == NULL)) { - fprintf(stderr, "must specify chip name when using -pixcenter\n"); - usage(); } @@ -125,8 +152,13 @@ if (!status) { // TODO: shall we print a specific message? + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify INPUT\n"); usage (); } - if (argc != 2) usage (); + // finally the output file + if (argc != 2) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "must specify OUTPUT\n"); + usage (); + } // Add the input and output images (which remain on the command-line) to the arguments list @@ -137,2 +169,59 @@ return config; } + +static bool validNumber(char *string) +{ + char *p = string; + if ((*p == '+') || (*p == '-')) { + p++; + } + return isdigit(*p); +} + +static bool get2Ints(int argnum, int *pArgc, char **argv, int *p1, int *p2) +{ + if (*pArgc < 2) { + usage(); + } + + if (!validNumber(argv[argnum])) { + return false; + } + + *p1 = atoi(argv[argnum]); + psArgumentRemove(argnum, pArgc, argv); + + if (!validNumber(argv[argnum])) { + return false; + } + + *p2 = atoi(argv[argnum]); + psArgumentRemove(argnum, pArgc, argv); + + return true; +} + +static bool get2Angles(int argnum, int *pArgc, char **argv, bool bothInDegrees, double *p1, double *p2) +{ + bool rval; + + if (*pArgc < 2) { + usage(); + } + + if (bothInDegrees) { + // both values are angles of arc DD:MM:SS or decimal degrees + rval = ohana_dms_to_ddd(p1, argv[argnum]); + if (rval) { + rval = ohana_dms_to_ddd(p1, argv[argnum+1]); + } + } else { + // first value may be in HH:MM:SS + rval = ohana_str_to_radec(p1, p2, argv[argnum], argv[argnum+1]); + } + + psArgumentRemove(argnum, pArgc, argv); + psArgumentRemove(argnum, pArgc, argv); + + return rval; +} Index: /trunk/ppstamp/src/ppstampMakeStamp.c =================================================================== --- /trunk/ppstamp/src/ppstampMakeStamp.c (revision 15322) +++ /trunk/ppstamp/src/ppstampMakeStamp.c (revision 15323) @@ -16,39 +16,62 @@ } ppstampOverlap; -static bool regionContainsPoint(psRegion *region, psPlane *point); -static bool regionContainsRegion(psRegion *outer, psRegion *inner); -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; -} +// 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) + + 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; + } + + return true; +} + +static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell, psRegion *roi) +{ + pmChip *outChip; + pmFPAview *view = pmFPAviewAlloc(0); + + // our output file has a single chip + view->chip = 0; + outChip = pmFPAviewThisChip(view, output->fpa); + psFree(view); + + // XXX we probably should copy some concepts from the chip (skipping those that don't apply) + outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts); + + pmHDU *inHDU = pmHDUGetHighest(input->fpa, inChip, NULL); + pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL); + + if (inHDU->header) { + outHDU->header = psMetadataCopy(outHDU->header, inHDU->header); + } else if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } + + // 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, roi)) { + return false; + } + + return true; +} + +// Build the postage stamp output file static bool makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input, @@ -82,5 +105,4 @@ while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { if (!readout->data_exists) { - // XXX if this happens something is wrong continue; } @@ -92,5 +114,4 @@ } - // XXX Do we have to do somethign to deal with negative parity psImage *subsetImage = psImageSubset(readout->image, options->roi); if (subsetImage) { @@ -99,4 +120,6 @@ psFree(subsetImage); if (outReadout->image) { + // The image has inherited the source's col0 and row0. We don't want this in + // our output so override it. // 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") @@ -126,122 +149,155 @@ -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); + +static bool regionContainsPoint(psRegion *r, psPlane *pt) +{ + if (pt->x < r->x0) + return false; + if (pt->x >= r->x1) + return false; + if (pt->y < r->y0) + return false; + if (pt->y >= r->y1) + return false; + + return true; +} + +// true if the inner region is equal to or completely contained in +// the outer region +static bool regionContainsRegion(psRegion *outer, psRegion *inner) +{ + if ((outer->x0 <= inner->x0) && + (outer->y0 <= inner->y0) && + (outer->x1 >= inner->x1) && + (outer->y1 >= inner->y1)) { + return true; } 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; - pmFPAview *view = pmFPAviewAlloc(0); - - // our output file has a single chip - view->chip = 0; - outChip = pmFPAviewThisChip(view, output->fpa); - psFree(view); - - outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts); - - pmHDU *inHDU = pmHDUGetHighest(input->fpa, inChip, NULL); - pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL); - - if (inHDU->header) { - outHDU->header = psMetadataCopy(outHDU->header, inHDU->header); - } else if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } - - // 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; -} - - -// 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); + return false; + } +} + + +static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip) +{ + // convert from sky to FP + psProject(pt->TP, pt->sky, fpa->toSky); + psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP); + // convert from FP to chip + psPlaneTransformApply(pt->chip, chip->fromFPA, pt->FP); +} + +static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip) +{ + psPlaneTransformApply(pt->FP, chip->toFPA, pt->chip); + psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP); + psDeproject(pt->sky, pt->TP, fpa->toSky); + +} + +static void compareToBox(psRegion *region, psPlane *pt) +{ + if (pt->x < region->x0) + region->x0 = pt->x; + if (pt->x > region->x1) + region->x1 = pt->x; + if (pt->y < region->y0) + region->y0 = pt->y; + if (pt->y > region->y1) + region->y1 = pt->y; +} + +static void findBoundingBox(ppstampOptions *options, pmFPA *fpa, pmChip *chip, pmAstromObj *center) +{ + pmAstromObj *pt = pmAstromObjAlloc(); + + if (!options->celestialCenter) { + // Center was specified in chip coordinates, transform to the sky + chipToSky(center, fpa, chip); + } + + // calculate the four corners of the bounding box in sky coordinates, translate them to + // chip coordinates and build the ROI by comparison. + + options->roi.x0 = INFINITY; + options->roi.x1 = -INFINITY; + options->roi.y0 = INFINITY; + options->roi.y1 = -INFINITY; + + pt->sky->rErr = 0; + pt->sky->dErr = 0; + + pt->sky->r = center->sky->r - options->dRA/2; + pt->sky->d = center->sky->d - options->dDEC/2; + skyToChip(pt, fpa, chip); + compareToBox(&options->roi, pt->chip); + + pt->sky->r = center->sky->r + options->dRA/2; + pt->sky->d = center->sky->d - options->dDEC/2; + skyToChip(pt, fpa, chip); + compareToBox(&options->roi, pt->chip); + + pt->sky->r = center->sky->r + options->dRA/2; + pt->sky->d = center->sky->d + options->dDEC/2; + skyToChip(pt, fpa, chip); + compareToBox(&options->roi, pt->chip); + + pt->sky->r = center->sky->r - options->dRA/2; + pt->sky->d = center->sky->d + options->dDEC/2; + skyToChip(pt, fpa, chip); + compareToBox(&options->roi, pt->chip); + + psFree(pt); +} + + +// find cell on given chip that contains the region of interest +// return the status of the overlap and if the cell completely contains the ROI +// set a pointer to reference the cell + +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 = ppstampCellRegion(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); + fprintf(stderr, "Partial Overlap cell %d\n", view->cell); + fprintf(stderr, " ROI: %s\n", psRegionToString(options->roi)); + fprintf(stderr, " Cell Extent: %s\n", psRegionToString(*cellExtent)); + 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; +} + +// findROI +// calculate the region of interest in chip coordinates and determine whether that region +// whether the region of interest is completely contained in a cell on a given chip. + +static ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view, pmFPAfile *input, pmChip *chip, + pmAstromObj *center, pmCell **ppCell) +{ + psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + psRegion *chipExtent = ppstampChipRegion(chip); bool onChip = false; ppstampOverlap returnval = PPSTAMP_OFF; - psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // set up the astrometry @@ -253,5 +309,4 @@ psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS\n"); psFree(chipExtent); - psFree(center); return PPSTAMP_ERROR; } @@ -259,14 +314,10 @@ if (options->celestialCenter) { - // convert from sky to FP center->sky->r = options->centerRA; center->sky->d = options->centerDEC; - center->sky->rErr = 0; // TODO: is this right? + center->sky->rErr = 0; center->sky->dErr = 0; - psProject(center->TP, center->sky, input->fpa->toSky); - psPlaneTransformApply(center->FP, input->fpa->fromTPA, center->TP); - - // convert from FP to chip + skyToChip(center, input->fpa, chip); psPlaneTransformApply(center->chip, chip->fromFPA, center->FP); @@ -277,16 +328,14 @@ } } else { - // center specified in pixels, user needs to have specified the chip - if (chipName == NULL) { - psError(PS_ERR_UNKNOWN, true, "Failed to find CHIP.NAME\n"); - returnval = PPSTAMP_ERROR; - } - // ppstampArguments insures that options->chipName is not null - if (!strcmp(chipName, options->chipName)) { + // center specified in pixels. + // If the user specified a name of a chip name wait until we get to that one. + // If no chip name was specified, select this one (the first one that had data) + if ((options->chipName == NULL) || !strcmp(chipName, options->chipName)) { psLogMsg("ppstampMakeStamp", 3, "Center on chip: %s\n", chipName); - onChip = true; - onChip = true; center->chip->x = options->centerX; center->chip->y = options->centerY; + center->chip->xErr = 0; + center->chip->yErr = 0; + onChip = true; } } @@ -294,8 +343,5 @@ if (onChip) { if (options->celestialRange) { - // Protect against unimplemented feature - psError(PS_ERR_PROGRAMMING, true, "not ready for Range in celestial coordinates yet.\n"); - exit(PS_EXIT_PROG_ERROR); - // find the bounding box in pixel space that encloses the ROI + findBoundingBox(options, input->fpa, chip, center); } else { int width = options->dX; @@ -303,18 +349,18 @@ // calculate the ROI in chip coordinates - roi->x0 = center->chip->x - width / 2; - roi->x1 = roi->x0 + width; - roi->y0 = center->chip->y - height / 2; - roi->y1 = roi->y0 + height; - - if (regionContainsRegion(chipExtent, roi)) { - returnval = findCell(view, options, input, center, ppCell); - } else { - 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; - } + options->roi.x0 = center->chip->x - width / 2; + options->roi.x1 = options->roi.x0 + width; + options->roi.y0 = center->chip->y - height / 2; + options->roi.y1 = options->roi.y0 + height; + } + + if (regionContainsRegion(chipExtent, &options->roi)) { + returnval = findCell(view, options, input, center, ppCell); + } else { + fprintf(stderr, "Partial Overlap chip %s\n", chipName); + fprintf(stderr, "ROI: %s\n", psRegionToString(options->roi)); + fprintf(stderr, "Chip Extent: %s\n", psRegionToString(*chipExtent)); + + returnval = PPSTAMP_PARTIALLY_ON; } } @@ -363,5 +409,5 @@ pmCell *cell; - ppstampOverlap overlap = findROI(options, view, input, chip, center, &options->roi, &cell); + ppstampOverlap overlap = findROI(options, view, input, chip, center, &cell); if (overlap == PPSTAMP_ON) { @@ -409,32 +455,2 @@ -static bool regionContainsPoint(psRegion *r, psPlane *pt) -{ - // TODO: should this comparison be done with integers since we - // are dealing with pixels? - // Maybe we should round and truncate before we get here. - if (pt->x < r->x0) - return false; - if (pt->x >= r->x1) - return false; - if (pt->y < r->y0) - return false; - if (pt->y >= r->y1) - return false; - - return true; -} - -// true if the inner region is equal to or completely contained in -// the outer region -static bool regionContainsRegion(psRegion *outer, psRegion *inner) -{ - if ((outer->x0 <= inner->x0) && - (outer->y0 <= inner->y0) && - (outer->x1 >= inner->x1) && - (outer->y1 >= inner->y1)) { - return true; - } else { - return false; - } -} Index: /trunk/ppstamp/src/ppstampOptions.c =================================================================== --- /trunk/ppstamp/src/ppstampOptions.c (revision 15322) +++ /trunk/ppstamp/src/ppstampOptions.c (revision 15323) @@ -17,8 +17,8 @@ options->celestialCenter = false; - options->centerX = NAN; - options->centerY = NAN; - options->centerRA = NAN; - options->centerDEC = NAN; + options->centerX = 0; + options->centerY = 0; + options->centerRA = 0; + options->centerDEC = 0; options->celestialRange = false; options->dX = 0; Index: /trunk/ppstamp/src/ppstampRegion.c =================================================================== --- /trunk/ppstamp/src/ppstampRegion.c (revision 15323) +++ /trunk/ppstamp/src/ppstampRegion.c (revision 15323) @@ -0,0 +1,70 @@ +#include + +#include "pmHDU.h" +#include "pmFPA.h" + +// Functions to calculate the image space boundaries of a given Chip or Cell +// These are calculated in coordinates that match the wcs transformation +// (adapted from pmFPAExtent.c) + +// return cell pixels bounding the readout +static psRegion *ppstampReadoutRegion(const pmReadout *readout) +{ + PS_ASSERT_PTR_NON_NULL(readout, NULL); + + psImage *image = readout->image; // Image from which to get dimensions + if (!image) { + return NULL; + } + + // This is the difference between this function and pmReadoutExtent. + // pmReadoutExtent ignores the col0, row0 of the readout + + int col0 = image->col0; + int row0 = image->row0; + + return psRegionAlloc(col0, col0 + image->numCols, + row0, row0 + image->numRows); +} + +// return chip pixels bounding the cell (all readouts) +psRegion *ppstampCellRegion(const pmCell *cell) +{ + PS_ASSERT_PTR_NON_NULL(cell, NULL); + + psArray *readouts = cell->readouts; // Array of component readouts + psRegion *cellExtent = psRegionAlloc(INFINITY, 0, INFINITY, 0); // Extent of cell + for (long i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // Readout of interest + psRegion *roExtent = ppstampReadoutRegion(readout); // Extent of readout + cellExtent->x0 = PS_MIN(cellExtent->x0, roExtent->x0); + cellExtent->x1 = PS_MAX(cellExtent->x1, roExtent->x1); + cellExtent->y0 = PS_MIN(cellExtent->y0, roExtent->y0); + cellExtent->y1 = PS_MAX(cellExtent->y1, roExtent->y1); + psFree(roExtent); + } + // Unlike pmCellExtent we don't add in cell.x0, cell.y0. The appropriate offset + // was included in the roExtent by ppstampReadoutRegion + + return cellExtent; +} + +// return chip pixels included in all cells +psRegion *ppstampChipRegion(const pmChip *chip) +{ + PS_ASSERT_PTR_NON_NULL(chip, NULL); + + psArray *cells = chip->cells; // Array of component cells + psRegion *chipExtent = psRegionAlloc(INFINITY, 0, INFINITY, 0); // Extent of chip + for (long i = 0; i < cells->n; i++) { + pmCell *cell = cells->data[i]; // Cell of interest + psRegion *cellExtent = ppstampCellRegion(cell); // Extent of cell + chipExtent->x0 = PS_MIN(chipExtent->x0, cellExtent->x0); + chipExtent->x1 = PS_MAX(chipExtent->x1, cellExtent->x1); + chipExtent->y0 = PS_MIN(chipExtent->y0, cellExtent->y0); + chipExtent->y1 = PS_MAX(chipExtent->y1, cellExtent->y1); + psFree(cellExtent); + } + + return chipExtent; +}