Changeset 15318

Timestamp:

Oct 15, 2007, 4:49:27 PM (19 years ago)

Author:

bills

Message:

Using mosaic was problematic because coordinate system changes,
so process readouts directly. Various other changes.

Location:

trunk/ppstamp/src

Files:

• ppstampCleanup.c (modified) (1 diff)
• ppstampMakeStamp.c (modified) (14 diffs)

trunk/ppstamp/src/ppstampCleanup.c

@@ -20 +20 @@
 
     // 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;

trunk/ppstamp/src/ppstampMakeStamp.c

@@ -5 +5 @@
 #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 +18 @@
 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 +73 @@
     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 +92 @@
         }
 
-        // 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 +99 @@
             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 +119 @@
     }
 
-    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 +209 @@
     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 +217 @@
         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 +272 @@
 
         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 +309 @@
 
             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 +341 @@
     // 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 +357 @@
 
         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-            // TODO log error
             psError(PS_ERR_UNKNOWN, false, "failed to load chip");
             status = false;
@@ -278 +362 @@
         }
 
-        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 +386 @@
             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;
         }
     }