Changeset 27840 for branches/simtest_nebulous_branches/ppStack/src/ppStackPrepare.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• ppStack/src (modified) (1 prop)
• ppStack/src/ppStackPrepare.c (modified) (12 diffs)

branches/simtest_nebulous_branches/ppStack/src

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

branches/simtest_nebulous_branches/ppStack/src/ppStackPrepare.c

@@ -26 +26 @@
     float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
     if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.RADIUS in recipe");
         return false;
     }
     float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
     if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.SIGMA in recipe");
         return false;
     }
     float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
     if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.FRAC in recipe");
         return false;
     }
@@ -42 +42 @@
     psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
     if (!mdok || !maskValStr) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find MASK.VAL in recipe");
         return false;
     }
@@ -53 +53 @@
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
     if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
+        psError(PPSTACK_ERR_DATA, false, "Unable to measure background for image");
         psFree(stats);
         psFree(rng);
@@ -126 +126 @@
     options->inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
     psVectorInit(options->inputMask, 0);
+    options->exposures = psVectorAlloc(options->num, PS_TYPE_F32);
+    psVectorInit(options->exposures, NAN);
 
     pmFPAfileActivate(config->files, false, NULL);
@@ -134 +136 @@
     }
 
-    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, "^PPSTACK.INPUT$");
-    psMetadataItem *fileItem; // Item from iteration
     psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
     int numCols = 0, numRows = 0;   // Size of image
-    int index = 0;                  // Index for file
-    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
-        assert(fileItem->type == PS_DATA_UNKNOWN);
-        pmFPAfile *inputFile = fileItem->data.V; // An input file
+    options->sumExposure = 0.0;
+    for (int i = 0; i < num; i++) {
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
+        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
+
+        options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE");
+        options->sumExposure += options->exposures->data.F32[i];
 
         // Get list of PSFs, to determine target PSF
         if (options->convolve) {
-            pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
+            pmChip *chip = pmFPAviewThisChip(view, file->fpa); // The chip: holds the PSF
             pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
             if (!psf) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
-                psFree(view);
-                psFree(fileIter);
+                psError(PPSTACK_ERR_PROG, false, "Unable to find PSF.");
+                psFree(view);
                 psFree(psfs);
                 return false;
             }
-            psfs->data[index] = psMemIncrRefCounter(psf);
+            psfs->data[i] = psMemIncrRefCounter(psf);
             psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF");
 
-            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
+            pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
             pmHDU *hdu = pmHDUFromCell(cell);
             assert(hdu && hdu->header);
@@ -163 +165 @@
             int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
             if (naxis1 <= 0 || naxis2 <= 0) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
-                psFree(view);
-                psFree(fileIter);
+                psError(PPSTACK_ERR_PROG, false, "Unable to determine size of image from PSF.");
+                psFree(view);
                 psFree(psfs);
                 return false;
@@ -176 +177 @@
 
 
-        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
-        psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
-        if (!sources) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
-            return NULL;
-        }
-        options->sourceLists->data[index] = psMemIncrRefCounter(sources);
+        bool redoPhot = psMetadataLookupBool(NULL, recipe, "PHOT");
+
+        pmDetections *detections = NULL;
+        if (options->convolve || options->matchZPs || options->photometry || redoPhot) {
+            pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout with sources
+            detections = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.DETECTIONS"); // Sources
+            if (!detections || !detections->allSources) {
+                psWarning("No detections found for image %d --- rejecting.", i);
+                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_CAL;
+                continue;
+            }
+            psAssert (detections->allSources, "missing sources?");
+
+            options->sourceLists->data[i] = psMemIncrRefCounter(detections->allSources);
+        }
 
         // Re-do photometry if we don't trust the source lists
-        if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
-            psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
+        if (redoPhot) {
+            psTrace("ppStack", 2, "Photometering input %d of %d....\n", i, num);
             pmFPAfileActivate(config->files, false, NULL);
-            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
-            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
+            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i);
+            if (options->convolve) {
+                pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
+            }
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File
             pmFPAview *photView = ppStackFilesIterateDown(config);
             if (!photView) {
@@ -198 +210 @@
             pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
 
-            if (!ppStackInputPhotometer(ro, sources, config)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
+            if (!ppStackInputPhotometer(ro, detections->allSources, config)) {
+                psError(psErrorCodeLast(), false, "Unable to do photometry on input sources");
                 psFree(view);
                 psFree(photView);
@@ -213 +225 @@
             ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
         }
-
-        index++;
-    }
-    psFree(fileIter);
+    }
+
+    psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
+    bool havePSFs = false;                                // Do we have any PSFs?
+    options->inputSeeing = psVectorAlloc(num, PS_TYPE_F32);
+    psVectorInit(options->inputSeeing, NAN);
+    for (int i = 0; i < num; i++) {
+        pmPSF *psf = psfs->data[i];     // PSF for image
+        if (!psf) {
+            continue;
+        }
+        havePSFs = true;
+
+        int xNum = PS_MAX(psf->trendNx, 1), yNum = PS_MAX(psf->trendNy, 1); // Number of realisations
+        double sumFWHM = 0.0;           // FWHM for image
+        int numFWHM = 0;                // Number of FWHM measurements
+        for (int y = 0; y < yNum; y++) {
+            float yPos = ((float)y + 0.5) / (float)yNum * numRows;
+            for (int x = 0; x < xNum; x++) {
+                float xPos = ((float)x + 0.5) / (float)xNum * numCols;
+                float fwhm = pmPSFtoFWHM(psf, xPos, yPos); // FWHM for image
+                if (isfinite(fwhm)) {
+                    sumFWHM += fwhm;
+                    numFWHM++;
+                }
+            }
+        }
+        if (numFWHM == 0) {
+            options->inputSeeing->data.F32[i] = NAN;
+            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+            psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", i);
+        } else {
+            options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM;
+        }
+
+        psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
+    }
+    if (havePSFs) {
+        psLogMsg("ppStack", PS_LOG_INFO, "%s", log);
+    }
+    psFree(log);
 
     // Generate target PSF
@@ -223 +272 @@
         psFree(psfs);
         if (!options->psf) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psError(psErrorCodeLast(), false, "Unable to determine output PSF.");
             psFree(view);
             return false;
@@ -229 +278 @@
         psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
                          "Target PSF for stack", options->psf);
-
-        pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.OUTPUT"); // Output chip
+        options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * numCols, 0.5 * numRows); // FWHM for target
+        psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
+
+        pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
         psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
                          "Target PSF", options->psf);
@@ -238 +289 @@
     // Zero point calibration
     if (!ppStackSourcesTransparency(options, view, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
+        psError(PPSTACK_ERR_DATA, false, "Unable to calculate transparency differences");
         psFree(view);
         return false;

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h