Changeset 10169 for trunk/psLib/src/imageops/psImageBackground.c

Timestamp:

Nov 22, 2006, 7:26:30 PM (19 years ago)

Author:

magnier

Message:

added psVectorStats for sky measurement, with multiple stats options; old code selected with ROBUST_QUARTILE

File:

• trunk/psLib/src/imageops/psImageBackground.c (modified) (4 diffs)

trunk/psLib/src/imageops/psImageBackground.c

@@ -11 +11 @@
 #include "psError.h"
 
-
+// XXX allow the user to choose the stats method?
+// (SAMPLE_MEAN, CLIPPED_MEAN, ROBUST_MEDIAN, FITTED_MEAN)
 psStats *psImageBackground(const psImage *image, const psImage *mask, psMaskType maskValue,
-                           double fmin, double fmax, long nMax, psRandom *rng)
+                           double fmin, double fmax, long nMax, psStatsOptions option, psRandom *rng)
 {
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
@@ -32 +33 @@
     long ny = image->numRows;
 
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE); // Statistics, for return
-
     int Nsubset = PS_MIN(nMax, nx*ny);  // Number of pixels in nubset
     int Npixels = nx*ny;                // Total number of pixels
@@ -43 +42 @@
     float max = values->data.F32[0];
 
+    // select a subset of the image pixels to measure the stats
     long n = 0;                         // Number of actual pixels in subset
     for (long i = 0; i < Nsubset; i++) {
@@ -62 +62 @@
     values->n = n;
 
-    int imin = fmin * n;
-    int imax = fmax * n;
-    int npts = imax - imin + 1;
+    psStats *stats = psStatsAlloc(option); // Statistics, for return
 
-    if (!psVectorSort(values, values)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
-        psFree(stats);
-        psFree(values);
-        return NULL;
+    if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+        // use simple stats code (old verions)
+        // XXX this hack is just for testing, drop when I am happy with the psStats version of the values
+
+        int imin = fmin * n;
+        int imax = fmax * n;
+        int npts = imax - imin + 1;
+
+        if (!psVectorSort(values, values)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
+            psFree(stats);
+            psFree(values);
+            return NULL;
+        }
+
+        // Subtract the median when we add the numbers, so we don't get numerical problems
+        float median = npts % 2 ? 0.5 * (values->data.F32[npts/2 - 1] + values->data.F32[npts/2]) :
+                       values->data.F32[npts/2];
+        double value = 0;
+        for (long i = imin; (i <= imax) && (i < n); i++) {
+            value += values->data.F32[i] - median;
+        }
+        value = value / npts + median;
+        stats->robustMedian = value;
+        stats->robustUQ = values->data.F32[imax];
+        stats->robustLQ = values->data.F32[imin];
+    } else {
+        // XXX leave this as a psphot user option (passed in as part of stats?)
+        if (stats->options & PS_STAT_FITTED_MEAN) {
+            stats->clipSigma = 1.0;
+        }
+        psVectorStats (stats, values, NULL, NULL, 0);
     }
-
-    // Subtract the median when we add the numbers, so we don't get numerical problems
-    float median = npts % 2 ? 0.5 * (values->data.F32[npts/2 - 1] + values->data.F32[npts/2]) :
-                   values->data.F32[npts/2];
-    double value = 0;
-    for (long i = imin; (i <= imax) && (i < n); i++) {
-        value += values->data.F32[i] - median;
-    }
-    value = value / npts + median;
-
-    stats->robustMedian = value;
-    stats->robustUQ = values->data.F32[imax];
-    stats->robustLQ = values->data.F32[imin];
 
     psFree(values);
     return stats;
-    // XXX correct for selection bias??
 }