Changeset 34735 for branches/eam_branches/ipp-20120905/pstamp/src/ppstampMakeStamp.c

Timestamp:

Nov 29, 2012, 2:15:01 PM (14 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20120905

Files:

• . (modified) (1 prop)
• pstamp/src/ppstampMakeStamp.c (modified) (15 diffs)

branches/eam_branches/ipp-20120905/pstamp/src/ppstampMakeStamp.c

@@ -19 +19 @@
 static void chipToSky(pmAstromObj *pt, pmFPA *fpa, pmChip *chip);
 static bool setMaskedToNAN(pmConfig *config, psImage *image, psImage *mask, psImage *variance);
+static bool copySources(pmReadout *outputReadout, pmReadout *inputReadout, pmReadout *astromReadout, psRegion extractRegion);
 
 // convert the input chip's transforms to the output
@@ -80 +81 @@
 }
 
-static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options, pmAstromObj *center)
+static bool copyMetadata(pmFPAfile *output, pmFPAfile *input, pmChip *inChip, ppstampOptions *options, pmAstromObj *center, pmFPAfile *astrom, pmFPAview *inView)
 {
     pmChip    *outChip;
-    pmFPAview *view = pmFPAviewAlloc(0);
+    pmFPAview *outView = pmFPAviewAlloc(0);
 
     // our output file has a single chip
-    view->chip = 0;
-    outChip = pmFPAviewThisChip(view, output->fpa);
-    psFree(view);
-
-    // copy data from the input chip header to the output.
+    outView->chip = 0;
+    outChip = pmFPAviewThisChip(outView, output->fpa);
+    psFree(outView);
+    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
+
+    // copy data from the input's chip header to the output.
     // since some of the keywords might be duplicated we may not want to copy both
 
-    pmHDU *inHDU  = pmHDUFromChip(inChip);
-    pmHDU *outHDU = pmHDUGetHighest(output->fpa, outChip, NULL);
+    pmHDU *inHDU  = NULL;
+    if (strcmp(options->stage, "raw") && input != astrom) {
+        // Copy the header from the astrometry file since it contains the psphot and psastro results
+        pmChip *astromChip = pmFPAviewThisChip(inView, astrom->fpa);
+        inHDU  = pmHDUFromChip(astromChip);
+    } else {
+        inHDU  = pmHDUFromChip(inChip);
+    }
 
     if (inHDU->header) {
@@ -102 +110 @@
     }
 
-    // copy the fpa concepts
+    // copy the fpa and chip concepts
     pmConceptsCopyFPA(output->fpa, input->fpa, false, false);
+    // Needed to preserve FPA.ZP
     pmConceptsCopyChip(outChip, inChip, false);
+
+    // copy the cell concepts (needed to get CELL.EXPOSURE which gets used to set EXPTIME)
     pmCell *outCell = outChip->cells->data[0]; // The only output cell
 
@@ -122 +133 @@
     }
 
-    // If input had WCS convert it for stamp
+    // If input had WCS convert it for the new image coordinate system of our stamp
     if (input->fpa->toSky) {
         // copy the input fpa's transforms
@@ -135 +146 @@
         }
     } else {
-        psWarning("No WCS present in input\n");
-    }
-
+        psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "No WCS present in input.");
+    }
+
+    // Save our specific metadata
     psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "RA_DEG", PS_META_REPLACE, "Right Ascension of stamp center", RAD_TO_DEG(center->sky->r));
     psMetadataAddF64(outHDU->header, PS_LIST_TAIL, "DEC_DEG", PS_META_REPLACE, "Declination of stamp center", RAD_TO_DEG(center->sky->d));
+
+    // XXX: TODO: save the region of interest boundary and the offset to be applied. 
+    // XXX: Gene thinks that there may be a way to handle this in the WCS
+
+    // Now a few sanity checks
+    // Make sure zero point in the fpa concepts matches the observed value in the input header
+    bool status;
+    float zpt_obs = psMetadataLookupF32(&status, outHDU->header, "ZPT_OBS");
+    if (status) {
+        psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.ZP", PS_META_REPLACE, "Magnitude zero point", zpt_obs);
+    }
 
     ppstampVersionMetadata(outHDU->header, options);
@@ -324 +347 @@
 
             if (options->censorMasked && !setMaskedToNAN(config, outReadout->image, outReadout->mask, outReadout->variance)) {
-                 psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
-                 status = false;
-                 break;
+                psError(PS_ERR_UNKNOWN, false, "failed to create postage stamp mask image\n");
+                status = false;
+                break;
             }
         }
@@ -335 +358 @@
         status = true;
 
+        if (options->writeCMF) {
+            pmReadout *astromReadout = NULL;
+            if (astrom->fpa != input->fpa) {
+                astromReadout = pmFPAviewThisReadout(srcView, astrom->fpa);
+            }
+            if (!copySources(outReadout, astromReadout, astromReadout, extractRegion)) {
+                psError(PS_ERR_UNKNOWN, false, "failed to extract sources from region of interest.\n");
+                status = false;
+                break;
+            }
+        }
+
         psFree(outReadout); // drop reference
     }
@@ -341 +376 @@
 
     if (status) {
-        status = copyMetadata(output, input, inChip, options, center);
+        // For raw stage we need to do more work to build a proper header
+        status = copyMetadata(output, input, inChip, options, center, astrom, view);
     }
     return status ? PS_EXIT_SUCCESS : PS_EXIT_UNKNOWN_ERROR;
@@ -428 +464 @@
     double dy = 0.5 * options->roip.dDEC / fpa->toSky->Ys;
 
-    if (dx > 5000) {
-        fprintf(stderr, "requested width %f too large reducing to 5000\n", dx);
-        dx = 5000;
-    }
-    if (dy > 5000) {
-        fprintf(stderr, "requested height %f too large reducing to 5000\n", dy);
-        dy = 5000;
+    // XXX: why do we limit this?
+    if (dx > 8000) {
+        psWarning( "requested width %f too large reducing to 8000\n", dx);
+        dx = 8000;
+    }
+    if (dy > 8000) {
+        psWarning( "requested height %f too large reducing to 8000\n", dy);
+        dy = 8000;
     }
 
@@ -480 +517 @@
     ppstampOverlap   returnval = PPSTAMP_OFF;
 
-//    fprintf(stderr, "ppstampChipBounds: %s\n", psRegionToString(*chipBounds));
-
     // set up the astrometry
     pmHDU *hdu = pmFPAviewThisHDU(view, astrom->fpa);
@@ -498 +533 @@
         }
     } else {
-        // we use a default FPA pixel scale of 1.0
         if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0) && mustHaveAstrometry) {
             psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry for input FPA.");
@@ -537 +571 @@
             center->chip->yErr = 0;
             onChip = true;
-            chipToSky(center, input->fpa, chip);
+            if (input->fpa->toSky) {
+                chipToSky(center, input->fpa, chip);
+            }
         }
     }
@@ -550 +586 @@
                 int width  = options->roip.dX;
                 int height = options->roip.dY;
-                if (width > 7000) {
-                    fprintf(stderr, "requested width %d too large reducing to 7000\n", width);
-                    width = 7000;
+                if (width > 8000) {
+                    psWarning( "requested width %d too large reducing to 8000\n", width);
+                    width = 8000;
                 }
-                if (height > 7000) {
-                    fprintf(stderr, "requested height %d too large reducing to 7000\n", height);
-                    height = 7000;
+                if (height > 8000) {
+                    psWarning( "requested height %d too large reducing to 8000\n", height);
+                    height = 8000;
                 }
 
@@ -675 +711 @@
 
     if (!foundOverlap && (returnval == PS_EXIT_SUCCESS)) {
-        fprintf(stderr, "ROI not found in input\n");
+        psWarning( "ROI not found in input\n");
         returnval = PSTAMP_NO_OVERLAP;
     }
@@ -739 +775 @@
     }
 
-    // XXX: this shouldn't be a fprintf
-    fprintf(stderr, "Censored %ld masked pixels\n", numCensored);
+    psLogMsg("ppstamp", PS_LOG_DETAIL, "Censored %ld masked pixels\n", numCensored);
 
     return true;
 }
+
+static bool copySources(pmReadout *outReadout, pmReadout *inReadout, pmReadout *astromReadout, psRegion extractRegion) {
+    bool status;
+
+    // first look for detections in the input readout. Those would have come from the -sources file
+    pmDetections *inDetections = psMetadataLookupPtr (&status, inReadout->analysis, "PSPHOT.DETECTIONS");
+    if (!inDetections) {
+        // No detections on the input readout. Try the astrometry readout if supplied.
+        if (astromReadout) {
+            inDetections = psMetadataLookupPtr (&status, astromReadout->analysis, "PSPHOT.DETECTIONS");
+            if (!inDetections) {
+                psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "no detections found on input or astrometry readout\n");
+                return false;
+            } else {
+                psLogMsg("ppstampMakeStamp", PS_LOG_INFO, "using detections from astrom file\n");
+            }
+        } else {
+            psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "no detections found on input readout\n");
+            return false;
+        }
+    }
+    if (!inDetections->allSources) {
+        psLogMsg("ppstampMakeStamp", PS_LOG_WARN, "no sources array found on detections structure\n");
+        return false;
+    }
+
+    // Ok we have sources, allocate a detections structure and copy them to the output readout
+    psArray *inSources = inDetections->allSources;
+
+    pmDetections *outDetections = pmDetectionsAlloc();
+    psMetadataAddPtr (outReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN,
+        "psphot detections", outDetections);
+
+    psArray *outSources = outDetections->allSources = psArrayAllocEmpty( 100 );
+
+    // include sources that are off the stamp by just a little bit
+    #define EXPAND_SIZE 16
+    psRegion sourceRegion = extractRegion;
+    sourceRegion.x0 -= EXPAND_SIZE;
+    sourceRegion.y0 -= EXPAND_SIZE;
+    sourceRegion.x1 += EXPAND_SIZE;
+    sourceRegion.y1 += EXPAND_SIZE;
+
+    for (int i = 0; i < inSources->n; i++ ) {
+        pmSource *src = inSources->data[i];
+
+        float x = src->peak->xf;
+        float y = src->peak->yf;
+        if (x < sourceRegion.x0 || x > sourceRegion.x1 ||
+            y < sourceRegion.y0 || y > sourceRegion.y1) {
+            continue;
+        }
+        psArrayAdd(outSources, 100, src);
+    }
+
+    return true;
+}