Index: /branches/eam_branch_20080719/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /branches/eam_branch_20080719/psModules/src/imcombine/pmSubtraction.c	(revision 18729)
+++ /branches/eam_branch_20080719/psModules/src/imcombine/pmSubtraction.c	(revision 18730)
@@ -27,5 +27,5 @@
 
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
-
+#define MIN_SAMPLE_STATS    7           // Minimum number to use sample statistics; otherwise use quartiles
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -513,25 +513,44 @@
     // with large N.  Therefore, let's just treat this as a Gaussian distribution.
 
-    double mean = 0.0;                  // Mean deviation
     int numStamps = 0;                  // Number of used stamps
+    psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask, for statistics
+    psVectorInit(mask, 0);
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
         if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            mask->data.PS_TYPE_MASK_DATA[i] = 0xff;
             continue;
         }
-        mean += deviations->data.F32[i];
         numStamps++;
     }
-    mean /= numStamps;
-
-    double rms = 0.0;                   // Standard deviation
-    for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
-            continue;
-        }
-        rms += PS_SQR(deviations->data.F32[i] - mean);
-    }
-    rms = sqrt(rms / (numStamps - 1));
+    psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", numStamps);
+
+    if (numStamps == 0) {
+        psError(PS_ERR_UNKNOWN, true, "No good stamps found.");
+        psFree(mask);
+        return -1;
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV |
+                                  PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics for deviatns
+    if (!psVectorStats(stats, deviations, NULL, mask, 0xff)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for deviations.");
+        psFree(stats);
+        psFree(mask);
+        return -1;
+    }
+    psFree(mask);
+
+    double mean, rms;                 // Mean and RMS of deviations
+    if (numStamps < MIN_SAMPLE_STATS) {
+        mean = stats->sampleMean;
+        rms = stats->sampleStdev;
+    } else {
+        mean = stats->sampleMedian;
+        rms = 0.74 * (stats->sampleUQ - stats->sampleLQ);
+    }
+
+    psTrace("psModules.imcombine", 1, "Mean: %f\n", mean);
+    psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms);
 
     if (rmsPtr) {
@@ -540,9 +559,4 @@
     if (numPtr) {
         *numPtr = numStamps;
-    }
-
-    if (numStamps == 0) {
-        psError(PS_ERR_UNKNOWN, true, "No good stamps found.");
-        return -1;
     }