Changeset 14282 for trunk/psModules/src/imcombine/pmStack.c

Timestamp:

Jul 17, 2007, 5:54:01 PM (19 years ago)

Author:

Paul Price

Message:

Use median as an estimator of the mean for rejection.

File:

• trunk/psModules/src/imcombine/pmStack.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-18 02:16:46 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-18 03:54:01 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -114 +114 @@
 // Generate a standard deviation for the combination
 // Not using psVectorStats because it has additional allocations which slow things down
-static bool combinationStdev(float *stdev, // Mean value, to return
+static bool combinationStdev(float *median, // Median value, to return
+                             float *stdev, // Standard deviation value, to return
                              const psVector *values, // Values to combine
                              const psVector *masks, // Mask to apply
@@ -137 +138 @@
         return false;
     }
+    *median = num % 2 ? (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 :
+        sortBuffer->data.F32[num / 2];
     *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
 
@@ -269 +272 @@
     psMaskType ignore = maskVal | bad;  // Ignore values with this mask value
     for (int i = 0; i < numIter && numClipped > 0; i++) {
-        float stdev;                    // Standard deviation of the combination, for rejection
+        float median, stdev;                    // Median and stdev of the combination, for rejection
         // Only get the mean if it's not the first iteration (because we got the mean earlier)
-        if ((i != 0 && !combinationMean(&mean, pixelData, weights, pixelMasks, maskVal)) ||
-            !combinationStdev(&stdev, pixelData, pixelMasks, maskVal, sort)) {
+        if (!combinationStdev(&median, &stdev, pixelData, pixelMasks, maskVal, sort)) {
             image->data.F32[y][x] = NAN;
             mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
@@ -285 +287 @@
                 continue;
             }
-            float diff = pixelData->data.F32[j] - mean;
+            float diff = pixelData->data.F32[j] - median;
             if (fabsf(diff) > limit) {
                 // Add the pixel as one to inspect