Changeset 31158 for trunk/ppStack/src/ppStackPrepare.c

Timestamp:

Apr 4, 2011, 1:16:41 PM (15 years ago)

Author:

eugene

Message:

improvements to the output headers; some code organization; plug leaks; user-specified limits to input seeing ranges

Location:

trunk/ppStack

Files:

• . (modified) (2 props)
• src/ppStackPrepare.c (modified) (9 diffs)

trunk/ppStack

@@ -17 +17 @@
 autom4te.cache
 Doxyfile
+a.out.dSYM

@@ -17 +17 @@
 autom4te.cache
 Doxyfile
+a.out.dSYM

trunk/ppStack/src/ppStackPrepare.c

@@ -1 @@
-#include "ppStack.h"
-
-#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
-                          PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
-                          PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
+# include "ppStack.h"
+
+# define RE_PHOTOMETER 0
+
+# define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
+                           PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
+                           PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
 
 bool ppStackInputPhotometer(const pmReadout *ro, const psArray *sources, const pmConfig *config)
@@ -13 +15 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
+# if (RE_PHOTOMETER)
     bool mdok;                          // Status of MD lookup
 
@@ -38 +41 @@
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
 
+    // Measure background
     psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
-
-    // Measure background
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
@@ -52 +54 @@
     psFree(stats);
     psFree(rng);
-
-    // Photometer sources
-    int numCols = image->numCols, numRows = image->numRows; // Size of image
+# endif
+
+    // Examine the sources (photometer if RE_PHOTOMETER is set)
     int numSources = sources->n; // Number of sources
     int numGood = 0;            // Number of good sources
@@ -64 +66 @@
             continue;
         }
-        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
+
+// This code re-measures the magnitude in the assumption that the warp analysis failed It would
+// be better to catch and correct such cases.  In any case, we should check on the validity of
+// this assumption (eg, compare the mean psf-ap offset) and raise an error on that?
+
+# if (RE_PHOTOMETER)
+        
+        int numCols = image->numCols, numRows = image->numRows; // Size of image
+        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
         int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
         // Range of integration
@@ -96 +106 @@
             }
         }
+# else
+        numGood++;
+        maxMag = PS_MAX(maxMag, source->psfMag);
+# endif
     }
     psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);
@@ -217 +231 @@
         }
     }
+
+    bool mdok = false;
+    float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
+    float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
 
     psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
@@ -251 +269 @@
         }
 
+        // reject any input images which exceed the specified max FWHM
+        if (isfinite(maxFWHM) && (options->inputSeeing->data.F32[i] > maxFWHM)) {
+            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+            psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f maxFWHM) --- rejected.", i, options->inputSeeing->data.F32[i], maxFWHM);
+        }
+
         psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
     }
@@ -257 +281 @@
     }
     psFree(log);
+
+    // We should have the ability to filter the input list based on the seeing:
+    // * reject above some max value and/or min value
+    // * reject images out-of-line with the rest of the images
+
+    // measure stats
+    psStats *fwhmStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+    psVectorStats (fwhmStats, options->inputSeeing, NULL, options->inputMask, 0xff);
+    psLogMsg("ppStack", PS_LOG_INFO, "Input FWHMs : %f +/- %f", fwhmStats->clippedMean, fwhmStats->clippedStdev);
+
+    if (isfinite(clipFWHMnSig)) {
+        float fwhmLimit = fwhmStats->clippedMean + clipFWHMnSig * fwhmStats->clippedStdev;
+        fwhmLimit = isfinite(maxFWHM) ? PS_MIN (maxFWHM, fwhmLimit) : fwhmLimit;
+        for (int i = 0; i < num; i++) {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+            if (options->inputSeeing->data.F32[i] > fwhmLimit) {
+                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+                psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f fwhmLimit) --- rejected.", i, options->inputSeeing->data.F32[i], fwhmLimit);
+            }
+        }
+    }
 
     // Generate target PSF

@@ -17 +17 @@
 autom4te.cache
 Doxyfile
+a.out.dSYM