Index: /branches/eam_branch_20071212/ppMerge/src/ppMergeMaskByImageStats.c
===================================================================
--- /branches/eam_branch_20071212/ppMerge/src/ppMergeMaskByImageStats.c	(revision 15806)
+++ /branches/eam_branch_20071212/ppMerge/src/ppMergeMaskByImageStats.c	(revision 15807)
@@ -15,4 +15,6 @@
 #include "ppMergeMask.h"
 
+// this function measures the statistics of the pixel values for a single readout (median and stdev)
+// and identifies pixels which deviate from the median by more than MASK.SUSPECT * stdev.
 
 // Generate a mask image consisting of pixels which are outliers from the image mean
@@ -65,6 +67,27 @@
                             view->chip, view->cell, view->readout);
                     float frac = options->sample / (float)(roIn->image->numCols * roIn->image->numRows);
+
+		    // measure the global image stats: median and stdev
+		    // XXX note that this now will accept any of several stats options
+		    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+		    stats->nSubsample = frac * image->numCols * image->numRows;
+		    if (!psImageBackground(stats, NULL, image, mask, maskVal, rng) ||
+			!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
+		      psError(PS_ERR_UNKNOWN, false, "Unable to measure image statistics.\n");
+		      psFree(stats);
+		      psFree(rng);
+		      return NULL;
+		    }
+		    psFree(rng);
+
+		    float median = stats->robustMedian; // Median value
+		    float stdev = stats->robustStdev; // Estimate of the standard deviation
+		    psFree(stats);
+
                     suspect = pmMaskFlagSuspectPixels(suspect, roIn, options->maskSuspect,
-                                                      options->combine->maskVal, frac, rng);
+                                                      options->combine->maskVal, median, stdev);
+		    
+                    psMetadataAddF32(roIn->analysis, PS_LIST_TAIL, "MEDIAN", 0, "median", median);
+                    psMetadataAddF32(roIn->analysis, PS_LIST_TAIL, "STDEV", 0, "stdev", stdev);
                 }
 
Index: /branches/eam_branch_20071212/ppMerge/src/ppMergeMaskReadoutByPixelStats.c
===================================================================
--- /branches/eam_branch_20071212/ppMerge/src/ppMergeMaskReadoutByPixelStats.c	(revision 15806)
+++ /branches/eam_branch_20071212/ppMerge/src/ppMergeMaskReadoutByPixelStats.c	(revision 15807)
@@ -9,4 +9,9 @@
 #include <pslib.h>
 #include <psmodules.h>
+
+// this function measures the statistics of the inputs->n representations of a given pixel (median and stdev) and identifies pixels
+// which deviate from the median by more than MASK.PIXEL.SUSPECT * stdev.
+
+// it also identifies pixels for which the stdev is more than MASK.PIXEL.NOISY * the image stdev
 
 // XXX: Maybe add support for S16 and S32 types.  Currently, only F32 supported.
@@ -27,13 +32,5 @@
     }
 
-    // XXX is it possible / desired to use the weight in this analysis?
-    for (int i = 0; i < inputs->n; i++) {
-        pmReadout *readout = inputs->data[i]; // Readout of interest
-        if (params->weights && !readout->weight) {
-            psError(PS_ERR_UNEXPECTED_NULL, true,
-                    "Rejection based on weights requested, but no weights supplied for image %d.\n", i);
-            return false;
-        }
-    }
+    // XXX note : we are ignoring weights for now.
 
     bool first = !output->image;        // First pass through?
@@ -77,43 +74,24 @@
     psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
 
+    // Get the correct statistics option for weights
     psStatsOptions combineStdev = 0; // Statistics option for weights
-    if (params->weights) {
-
-        if (!output->weight) {
-            output->weight = psImageAlloc(xSize, ySize, PS_TYPE_F32);
-        }
-        if (output->weight->numCols < xSize || output->weight->numRows < ySize) {
-            psImage *newWeight = psImageAlloc(xSize, ySize, PS_TYPE_F32);
-            psImageInit(newWeight, 0.0);
-            psImageOverlaySection(newWeight, output->weight, output->col0, output->row0, "=");
-            psFree(output->weight);
-            output->weight = newWeight;
-        }
-
-        if (first) {
-            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
-                             "Using input weights to combine images", "");
-        }
-
-        // Get the correct statistics option for weights
-        switch (params->combine) {
-        case PS_STAT_SAMPLE_MEAN:
-        case PS_STAT_SAMPLE_MEDIAN:
-            combineStdev = PS_STAT_SAMPLE_STDEV;
-            break;
-        case PS_STAT_ROBUST_MEDIAN:
-            combineStdev = PS_STAT_ROBUST_STDEV;
-            break;
-        case PS_STAT_FITTED_MEAN:
-            combineStdev = PS_STAT_FITTED_STDEV;
-            break;
-        case PS_STAT_CLIPPED_MEAN:
-            combineStdev = PS_STAT_CLIPPED_STDEV;
-            break;
-        default:
-            psAbort("Should never get here --- checked params->combine before.\n");
-        }
-        stats->options |= combineStdev;
-    }
+    switch (params->combine) {
+      case PS_STAT_SAMPLE_MEAN:
+      case PS_STAT_SAMPLE_MEDIAN:
+	combineStdev = PS_STAT_SAMPLE_STDEV;
+	break;
+      case PS_STAT_ROBUST_MEDIAN:
+	combineStdev = PS_STAT_ROBUST_STDEV;
+	break;
+      case PS_STAT_FITTED_MEAN:
+	combineStdev = PS_STAT_FITTED_STDEV;
+	break;
+      case PS_STAT_CLIPPED_MEAN:
+	combineStdev = PS_STAT_CLIPPED_STDEV;
+	break;
+      default:
+	psAbort("Should never get here --- checked params->combine before.\n");
+    }
+    stats->options |= combineStdev;
 
     // We loop through each pixel in the output image.  We loop through each input readout.  We determine if
@@ -129,11 +107,4 @@
     psU8 *maskData = mask->data.U8;     // Dereference mask
 
-    psVector *weights = NULL;           // Stack of weights
-    psVector *errors = NULL;            // Stack of errors (sqrt of variance/weights), for psVectorStats
-    psF32 *weightsData = NULL;          // Dereference weights
-    if (params->weights) {
-        weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights
-        weightsData = weights->data.F32;
-    }
     psVector *index = NULL;             // The indices to sort the pixels
 
@@ -172,8 +143,4 @@
     psF32 **outputImage  = output->image->data.F32; // Output image
     psU8  **outputMask   = output->mask->data.U8; // Output mask
-    psF32 **outputWeight = NULL; // Output weight map
-    if (output->weight) {
-        outputWeight = output->weight->data.F32;
-    }
 
     psVector *invScale = NULL;          // Inverse scale; pre-calculated for efficiency