Index: trunk/ppMerge/src/ppMergeMask.c
===================================================================
--- trunk/ppMerge/src/ppMergeMask.c	(revision 15862)
+++ trunk/ppMerge/src/ppMergeMask.c	(revision 15913)
@@ -43,36 +43,55 @@
         pmChip *chipIn;                 // Input chip of interest
         while ((chipIn = pmFPAviewNextChip(view, fpaIn, 1))) {
+
+	    // handle chip vs cell statistics & avoid reading the data twice
+
+	    // load the data of all cells 
             pmCell *cellIn;             // Input cell of interest
             while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
-                if (!pmCellRead(cellIn, fits, config->database)) {
-                    continue;
-                }
-                if (cellIn->readouts->n == 0) {
-                    continue;
-                }
-
-                pmCell *cellOut = pmFPAviewThisCell(view, fpaOut); // Output cell
-                // Suspect pixels image
-                psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT");
-                bool first = suspect ? false : true;
-
-                pmReadout *roIn;        // Input readout of interest
-                while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
-                    if (!roIn->image) {
-                        continue;
-                    }
-                    psTrace("ppMerge", 4, "Flagging suspect pixels in chip %d, cell %d, ro %d\n",
-                            view->chip, view->cell, view->readout);
-                    float frac = options->sample / (float)(roIn->image->numCols * roIn->image->numRows);
-                    suspect = pmMaskFlagSuspectPixels(suspect, roIn, options->maskSuspect,
-                                                      options->combine->maskVal, frac, rng);
-                }
-
-                if (first) {
-                    psMetadataAddImage(cellOut->analysis, PS_LIST_TAIL, "MASK.SUSPECT", 0,
-                                       "Suspect pixels", suspect);
-                    psFree(suspect);
-                }
-
+                if (!pmCellRead(cellIn, fits, config->database)) continue;
+                if (cellIn->readouts->n == 0) continue;
+            }
+
+	    // calculate the readout statistics either for each readout, or across the entire chip
+	    if (options->statsByChip) {
+		ppMergeMaskChipStats (chipIn, options, rng);
+	    } else {
+		// calculate the stats for each cell independently
+		while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
+		    if (cellIn->readouts->n == 0) continue;
+		    pmReadout *roIn;        // Input readout of interest
+		    while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
+			if (!roIn->image) continue;
+			psTrace("ppMerge", 4, "Measure statistics for chip %d, cell %d, ro %d\n",
+				view->chip, view->cell, view->readout);
+			ppMergeMaskReadoutStats (roIn, options, rng);
+		    }
+		}
+	    }
+
+	    // apply the measured statistics to determine the outliers to be masked
+	    while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
+		if (cellIn->readouts->n == 0) continue;
+
+		pmCell *cellOut = pmFPAviewThisCell(view, fpaOut); // Output cell
+		// Suspect pixels image
+		psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT");
+		bool first = suspect ? false : true;
+
+		pmReadout *roIn;        // Input readout of interest
+		while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
+		    if (!roIn->image) {
+			continue;
+		    }
+		    psTrace("ppMerge", 4, "Flagging suspect pixels in chip %d, cell %d, ro %d\n",
+			    view->chip, view->cell, view->readout);
+		    suspect = pmMaskFlagSuspectPixels(suspect, roIn, options->maskSuspect, options->combine->maskVal);
+		}
+
+		if (first) {
+		    psMetadataAddImage(cellOut->analysis, PS_LIST_TAIL, "MASK.SUSPECT", 0,
+				       "Suspect pixels", suspect);
+		    psFree(suspect);
+		}
                 pmCellFreeData(cellIn);
             }
@@ -195,2 +214,141 @@
     return true;
 }
+
+bool ppMergeMaskReadoutStats(const pmReadout *readout, 
+			     ppMergeOptions *options, // Options
+			     psRandom *rng)
+{
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
+    PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false);
+    PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
+    if (readout->mask) {
+        PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false);
+        PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false);
+    }
+
+    psImage *image = readout->image;    // Image of interest
+    psImage *mask = readout->mask;      // Corresponding mask
+
+    if (rng) {
+        psMemIncrRefCounter(rng);
+    } else {
+        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
+
+    // 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 = options->sample;
+
+    if (!psImageBackground(stats, NULL, image, mask, options->combine->maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Failure to measure image statistics.\n");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+    if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
+        psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
+
+    psFree(rng);
+    return true;
+}
+
+bool ppMergeMaskChipStats (const pmChip *chip,
+			   ppMergeOptions *options,
+			   psRandom *rng)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+
+    // XXX note that this now will accept any of several stats options
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+
+    if (rng) {
+	psMemIncrRefCounter(rng);
+    } else {
+	rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
+
+    // accumulate a vector of data values using options->sample per readout
+    psVector *values = psVectorAllocEmpty(options->sample, PS_TYPE_F32); // Vector containing subsample
+
+    for (int nCell = 0; nCell < chip->cells->n; nCell++) {
+	pmCell *cell = chip->cells->data[nCell];
+	if (cell->readouts->n == 0) continue;
+
+	for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
+	    pmReadout *readout = cell->readouts->data[nReadout];
+	    if (!readout->image) continue;
+
+	    psTrace("ppMerge", 4, "Measure statistics for cell %d, readout %d\n", nCell, nReadout);
+
+	    psImage *image = readout->image;    // Image of interest
+	    psImage *mask = readout->mask;      // Corresponding mask
+
+	    // Size of image
+	    long nx = image->numCols;
+	    long ny = image->numRows;
+	    const int Npixels = nx*ny;	// Total number of pixels
+	    const int Nsubset = PS_MIN(options->sample, Npixels); // Number of pixels in subset
+	    
+	    values = psVectorRealloc (values, values->n + Nsubset); // make sure we have enough space
+
+	    // select a subset of the image pixels to measure the stats
+	    for (long i = 0; i < Nsubset; i++) {
+		double frnd = psRandomUniform(rng);
+		int pixel = Npixels * frnd;
+		int ix = pixel % nx;
+		int iy = pixel / nx;
+
+		if (!isfinite(image->data.F32[iy][ix])) continue;
+		if (mask && (mask->data.U8[iy][ix] & options->combine->maskVal)) continue;
+
+		float value = image->data.F32[iy][ix];
+		values->data.F32[values->n] = value;
+		values->n ++;
+	    }
+	}
+    }
+
+    // calculate the statistics
+    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for chip");
+	psFree(values);
+	psFree(stats);
+	psFree(rng);
+	return false;
+    }
+    if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
+	psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
+	psFree(values);
+	psFree(stats);
+	psFree(rng);
+	return false;
+    }
+
+    // supply the stats to the readout analysis metadata
+    for (int nCell = 0; nCell < chip->cells->n; nCell++) {
+	pmCell *cell = chip->cells->data[nCell];
+	if (cell->readouts->n == 0) continue;
+
+	for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
+	    pmReadout *readout = cell->readouts->data[nReadout];
+	    if (!readout->image) continue;
+
+	    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
+	    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
+	}
+    }
+
+    psFree(values);
+    psFree(stats);
+    psFree(rng);
+    return true;
+}