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

Timestamp:

Oct 23, 2007, 11:50:52 AM (19 years ago)

Author:

bills

Message:

Extract postage stamp from chip mosaic image. Changed -celcenter to -skycenter
Range in celestial coordinates not specified with -arcrange (arcseconds)
Better error handling.
Allow files without WCS if we don't need it

Location:

trunk/ppstamp/src

Files:

• . (modified) (1 prop)
• ppstampMakeStamp.c (modified) (23 diffs)

trunk/ppstamp/src

@@ -6 +6 @@
 config.h.in
 stamp-h1
+ppstampErrorCodes.c
+ppstampErrorCodes.h

@@ -6 +6 @@
 config.h.in
 stamp-h1
+ppstampErrorCodes.c
+ppstampErrorCodes.h

trunk/ppstamp/src/ppstampMakeStamp.c

@@ -16 +16 @@
 } ppstampOverlap;
 
+static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
+
 // convert the input chip's transforms to the output
 static bool convertWCS(pmHDU *outHDU, pmFPA *outFPA, pmChip *outChip, pmChip *inChip, psRegion *roi)
@@ -24 +26 @@
 
     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");
+        psError(PS_ERR_UNKNOWN, false, "Failed to write WCS to output\n");
         return false;
     }
@@ -34 +37 @@
 }
 
-static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, pmCell *inCell, psRegion *roi)
+static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options)
 {
     pmChip    *outChip;
@@ -44 +47 @@
     psFree(view);
 
-    // XXX we probably should copy some concepts from the chip (skipping those that don't apply)
+    // copy data from the input chip header to the output.
+    // since some of the keywords might be duplicated we may not want to copy both
+
+    // copy the fpa concepts
     outChip->parent->concepts = psMetadataCopy(outChip->parent->concepts, inChip->parent->concepts);
 
-    pmHDU *inHDU  = pmHDUGetHighest(input->fpa, inChip, NULL);
+    pmHDU *inHDU  = pmHDUFromChip(inChip);
     pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
 
@@ -56 +62 @@
     }
 
-    // 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;
-    }
+    // 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)) {
+            return false;
+        }
+    } else {
+        psWarning("No WCS present in input\n");
+    }
+
+    ppstampVersionMetadata(outHDU->header, options);
 
     return true;
 }
 
+#ifdef notdef
+// update the ROI to account for any change in the coordinate system in the
+// mosaic process
+// THERE isn't any change. This is something I was trying to make the megacam
+// transforms come out correctly.
+static bool updateROI(ppstampOptions *options, pmFPAfile *mosaic, pmChip *mosaicChip)
+{
+    // set up the astrometry
+    pmHDU *hdu = pmHDUFromChip(mosaicChip);
+    PS_ASSERT_PTR_NON_NULL(hdu, 1)
+    PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
+
+
+    if (!pmAstromReadWCS(mosaic->fpa, mosaicChip, hdu->header, 1.0)) {
+        psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS from mosaic, cannot proceed\n");
+        return false;
+    }
+    pmAstromObj *center = pmAstromObjAlloc();
+
+    center->sky->r = options->centerRA;
+    center->sky->d = options->centerDEC;
+    center->sky->rErr = 0;
+    center->sky->dErr = 0;
+
+    skyToChip(center, mosaic->fpa, mosaicChip);
+
+    int width = options->dX;
+    int height = options->dY;
+
+    fprintf(stderr, "ROI before: %s\n", psRegionToString(options->roi));
+
+    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;
+
+    fprintf(stderr, "ROI after:  %s\n", psRegionToString(options->roi));
+
+
+    return true;
+}
+#endif
+
 // Build the postage stamp output file
 
 static bool makeStamp(pmConfig *config, ppstampOptions *options, pmFPAfile *input, 
-                pmChip *inChip, pmCell *inCell, pmFPAview *view)
+                pmChip *inChip, pmFPAview *view)
 {
     int status = false;
@@ -91 +147 @@
     outview->chip = 0;
     outview->cell = 0;
-    outview->readout = 0;
     pmCell *target =  pmFPAviewThisCell(outview, output->fpa); // Target cell
     psFree(outview);
     outview = NULL;
 
-    //    psMetadataPrint(stderr, inChip->concepts, 0);
-    //   psMetadataPrint(stderr, inCell->concepts, 0);
-    // these default to zero would that be ok?
+    //   psMetadataPrint(stderr, inChip->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 the tricky bits of parity, binning, etc.
+    pmFPAfile *mosaic = ppstampBuildMosaic(config, input, view);
+    if (mosaic == NULL) {
+        return false;
+    }
+
+    // one cell one chip
+    pmFPAview *mosaicView = pmFPAviewAlloc(0);
+    mosaicView->chip = 0;
+    mosaicView->cell = 0;
+
     pmReadout *readout;
-    while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
+    while ((readout = pmFPAviewNextReadout (mosaicView, mosaic->fpa, 1)) != NULL) {
         if (!readout->data_exists) {
+            psError(PS_ERR_UNKNOWN, false, "no data in mosaic readout!\n");
             continue;
         }
@@ -113 +182 @@
             break;
         }
-
+    
         psImage *subsetImage = psImageSubset(readout->image, options->roi);
         if (subsetImage) {
@@ -120 +189 @@
             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.
+                // The image has inherited the subset's col0 and row0. We don't want this in
+                // our output so override the values.
                 // 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")
+                // (What I really needed was a function that does
+                //     "create a copy of this portion of an image and return it as a new image with
+                //      origin 0,0 ")
                 P_PSIMAGE_SET_COL0(outReadout->image, 0);
                 P_PSIMAGE_SET_ROW0(outReadout->image, 0);
 
                 outReadout->data_exists = true;
-                outReadout->parent->data_exists = true;         // cell
-                outReadout->parent->parent->data_exists = true; // chip
+                outReadout->parent->data_exists = true;
+                outReadout->parent->parent->data_exists = true;
                 status = true;
             } else {
@@ -142 +213 @@
     }
 
+    // XXXX remove this or make writing the mosaic to a file a valid option
+    static bool writeMosaic = false;
+    if (writeMosaic) {
+        mosaic->save = true;
+        mosaicView->cell = -1;
+        pmFPAfileIOChecks(config, mosaicView, PM_FPA_AFTER);
+    }
+
+    psFree(mosaicView);
+
     if (status) {
-        status = copyMetadata(output, input, inChip, inCell, &options->roi);
+        status = copyMetadata(output, input, inChip, options);
     }
     return status;
@@ -181 +262 @@
 static void skyToChip(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
 {
-    // convert from sky to FP
+    // convert from sky to TP to FP
     psProject(pt->TP, pt->sky, fpa->toSky);
     psPlaneTransformApply(pt->FP, fpa->fromTPA, pt->TP);
@@ -190 +271 @@
 static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip)
 {
+    // chip to FP
     psPlaneTransformApply(pt->FP, chip->toFPA, pt->chip);
+    // FP to TP to sky
     psPlaneTransformApply(pt->TP, fpa->toTPA, pt->FP);
     psDeproject(pt->sky, pt->TP, fpa->toSky);
@@ -208 +291 @@
 }
 
+// For roi width and height given in sky coordinates find a bounding box in chip coordinates
+// that encloses the requested range
 static void findBoundingBox(ppstampOptions *options, pmFPA *fpa, pmChip *chip, pmAstromObj *center)
 {
@@ -249 +334 @@
 
     psFree(pt);
-}
-
-
+
+    // save the width and height in case we need them later
+    options->dX = options->roi.x1 - options->roi.x0;
+    options->dY = options->roi.y1 - options->roi.y0;
+}
+
+#ifdef notdef
+// I'm not using this any more. We mosaic the chip before extracting the pixels
+// so we don't need to deal with cells at this level.
+
+// findCell
 // 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 
@@ -263 +356 @@
     view->cell = -1;
     while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
-        psRegion *cellExtent = ppstampCellRegion(cell);
+        psRegion *cellExtent = pmCellExtent(cell);
         if (regionContainsPoint(cellExtent, center->chip)) {
             if (regionContainsRegion(cellExtent, &options->roi)) {
@@ -288 +381 @@
     return PPSTAMP_OFF;
 }
+#endif
+
+
 
 // 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.
+// is completely contained on a single chip
 
 static ppstampOverlap findROI(ppstampOptions *options, pmFPAview *view, pmFPAfile *input, pmChip *chip,
-                       pmAstromObj *center, pmCell **ppCell)
+                       pmAstromObj *center)
 {
     psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
-    psRegion    *chipExtent = ppstampChipRegion(chip);
+    psRegion    *chipBounds = ppstampChipRegion(chip);
     bool        onChip = false;
     ppstampOverlap   returnval = PPSTAMP_OFF;
+
+//    fprintf(stderr, "ppstampChipBounds: %s\n", psRegionToString(*chipBounds));
 
     // set up the astrometry
@@ -307 +405 @@
 
     if (!pmAstromReadWCS(input->fpa, chip, hdu->header, 1.0)) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS\n");
-        psFree(chipExtent);
-        return PPSTAMP_ERROR;
+        // we can live without WCS if the ROI is specified in pixels
+
+        if (options->celestialCenter || options->celestialRange) {
+            psError(PPSTAMP_ERR_DATA, false, "Failed to Read WCS, cannot proceed\n");
+            psFree(chipBounds);
+            return PPSTAMP_ERROR;
+        }
     }
 
@@ -320 +422 @@
 
         skyToChip(center, input->fpa, chip);
-        psPlaneTransformApply(center->chip, chip->fromFPA, center->FP);
-
-        if (regionContainsPoint(chipExtent, center->chip)) {
-            psLogMsg("ppstampMakeStamp", 3, "Found center (%f %f) on chip: %s\n", 
+
+        if (regionContainsPoint(chipBounds, center->chip)) {
+            psLogMsg("ppstampMakeStamp", 2, "Found center (%f %f) on chip: %s\n", 
                 center->chip->x, center->chip->y, chipName);
             onChip = true;
@@ -332 +433 @@
         // 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);
+            psLogMsg("ppstampMakeStamp", 2, "Center on chip: %s\n", chipName);
             center->chip->x = options->centerX;
             center->chip->y = options->centerY;
@@ -355 +456 @@
         }
 
-        if (regionContainsRegion(chipExtent, &options->roi)) {
-            returnval = findCell(view, options, input, center, ppCell);
+
+        if (regionContainsRegion(chipBounds, &options->roi)) {
+            // returnval = findCell(view, options, input, center, ppCell);
+            psLogMsg("ppstampMakeStamp", 2, "ROI contained on: %s\n", chipName);
+            returnval = PPSTAMP_ON;
         } 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));
+            fprintf(stderr, "Chip Extent: %s\n", psRegionToString(*chipBounds));
 
             returnval = PPSTAMP_PARTIALLY_ON;
-        }
-    }
-    psFree(chipExtent);
+
+        }
+    }
+    psFree(chipBounds);
 
     return returnval;
@@ -374 +479 @@
     bool        status = false;
     bool        returnval = false;;
+    bool        foundOverlap = false;
     pmAstromObj *center = pmAstromObjAlloc();
 
@@ -408 +514 @@
         }
 
-        pmCell *cell;
-        ppstampOverlap overlap = findROI(options, view, input, chip, center, &cell);
-
-        if (overlap == PPSTAMP_ON) {
-
-            returnval = makeStamp(config, options, input, chip, cell, view);
-
+        ppstampOverlap overlap = findROI(options, view, input, chip, center);
+
+        switch (overlap) {
+        case PPSTAMP_OFF:
+            // keep looking
+            break;
+        case PPSTAMP_ON:
+            returnval = makeStamp(config, options, input, chip, view);
             allDone = true;
-        } else if (overlap == PPSTAMP_PARTIALLY_ON) {
-            psError(PS_ERR_UNKNOWN, false, "Region of interest must be confined to a single cell\n");
+            foundOverlap = true;
+            break;
+        case PPSTAMP_PARTIALLY_ON:
+            psError(PS_ERR_UNKNOWN, false, "Region of interest must be confined to a single chip\n");
             // XXX: perhaps print a helpful message about what coordinates would be valid
             returnval = false;
             allDone = true;
-        } else if (overlap == PPSTAMP_ERROR) {
+            foundOverlap = true;
+            break;
+        case PPSTAMP_ERROR:
             returnval = false;
             allDone = true;
+            break;
+        default:
+            psError(PS_ERR_PROGRAMMING, false, "findROI returned unexpected value %d\n", overlap);
+            break;
         }
 
@@ -430 +545 @@
         if (allDone) {
             view->chip = -1;
-            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;
+            break;
         }
     } 
@@ -451 +553 @@
     psFree(view);
 
+    if (!foundOverlap) {
+        fprintf(stderr, "ROI not found in input\n");
+    }
+
     return returnval;
 }

@@ -6 +6 @@
 config.h.in
 stamp-h1
+ppstampErrorCodes.c
+ppstampErrorCodes.h