Changeset 15323 for trunk/ppstamp/src/ppstampMakeStamp.c

Timestamp:

Oct 16, 2007, 4:00:33 PM (19 years ago)

Author:

bills

Message:

Solved problem with chip and cell boundary calculations by using different
functions than pmChipPixels and pmCellExtents. They don't do their calculations
in image coordinates.

Implemented specifying width/height of the ROI in celestial coordinates.
Removed unneeded WCS calculation methods. (The bug I was chasing was actually
not in ppstamp. It was a bug in ds9

File:

• trunk/ppstamp/src/ppstampMakeStamp.c (modified) (12 diffs)

trunk/ppstamp/src/ppstampMakeStamp.c

@@ -16 +16 @@
 } 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 +105 @@
     while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
         if (!readout->data_exists) {
-            // XXX if this happens something is wrong
             continue;
         }
@@ -92 +114 @@
         }
 
-        // XXX Do we have to do somethign to deal with negative parity
         psImage *subsetImage = psImageSubset(readout->image, options->roi);
         if (subsetImage) {
@@ -99 +120 @@
             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 +149 @@
 
 
-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 +309 @@
         psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS\n");
         psFree(chipExtent);
-        psFree(center);
         return PPSTAMP_ERROR;
     }
@@ -259 +314 @@
     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 +328 @@
         }
     } 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 +343 @@
     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 +349 @@
 
             // 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 +409 @@
 
         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 +455 @@
 
 
-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;
-    }
-}