Index: trunk/ppMerge/src/ppMergeMask.c
===================================================================
--- trunk/ppMerge/src/ppMergeMask.c	(revision 15913)
+++ trunk/ppMerge/src/ppMergeMask.c	(revision 15937)
@@ -22,193 +22,30 @@
     )
 {
-    psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
-    assert(filenames);                  // It should be here --- it's put here in ppMergeConfig
-
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
-    pmFPA *fpaOut = data->out;          // Output FPA
-
-    // Iterate over each file
-    for (int i = 0; i < filenames->n; i++) {
-        if (! filenames->data[i] || strlen(filenames->data[i]) == 0) {
-            continue;
-        }
-        psFits *fits = data->files->data[i]; // FITS file handle
-        if (!fits) {
-            continue;
-        }
-        psTrace("ppMerge", 3, "File %d: %s\n", i, (const char*)filenames->data[i]);
-
-        pmFPA *fpaIn = data->in->data[i]; // Input FPA
-        pmFPAview *view = pmFPAviewAlloc(0); // View of FPA, for iteration
-        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;
-            }
-
-	    // 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);
-            }
-            pmChipFreeData(chipIn);
-        }
-        pmFPAFreeData(fpaIn);
-        psFree(view);
-    }
-    psFree(rng);
-
-    pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration
-    if (fpaOut->hdu) {
-        pmFPAUpdateNames(fpaOut, NULL, NULL);
-    }
-    pmFPAWriteMask(fpaOut, data->outFile, config->database, true, false); // Write header only
-    pmChip *chipOut;                    // Output chip of interest
-    while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) {
-        if (chipOut->hdu) {
-            chipOut->data_exists = true;
-            pmFPAUpdateNames(fpaOut, chipOut, NULL);
-        }
-        pmChipWriteMask(chipOut, data->outFile, config->database, true, false); // Write header only
-        pmCell *cellOut;                   // Output cell of interest
-        while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) {
-            if (cellOut->hdu) {
-                chipOut->data_exists = cellOut->data_exists = true;
-                pmFPAUpdateNames(fpaOut, chipOut, cellOut);
-            }
-            pmCellWriteMask(cellOut, data->outFile, config->database, true); // Write header only
-
-            psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT");
-            if (! suspect) {
-                continue;
-            }
-
-            pmReadout *roOut = pmReadoutAlloc(cellOut); // Output readout
-            roOut->mask = pmMaskIdentifyBadPixels(suspect, options->combine->maskVal, filenames->n, options->maskBad, options->maskMode);
-            roOut->data_exists = cellOut->data_exists = chipOut->data_exists = true;
-
-            // Get list of cells for concepts averaging
-            {
-                psList *inCells = psListAlloc(NULL); // List of cells
-                for (int i = 0; i < filenames->n; i++) {
-                    if (! filenames->data[i] || strlen(filenames->data[i]) == 0) {
-                        continue;
-                    }
-                    pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell
-                    psListAdd(inCells, PS_LIST_TAIL, cellIn);
-                }
-                if (!pmConceptsAverageCells(cellOut, inCells, NULL, NULL, true)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
-                    psFree(inCells);
-                    return false;
-                }
-                psFree(inCells);
-            }
-
-            // Statistics on the merged cell
-            if (data->statsFile) {
-                if (!data->stats) {
-                    data->stats = psMetadataAlloc();
-                }
-
-                // Build a fake image to do statistics
-                roOut->image = psImageAlloc(roOut->mask->numCols, roOut->mask->numRows, PS_TYPE_F32);
-                psImageInit(roOut->image, 1.0);
-                if (!ppStatsFPA(data->stats, data->out, view,
-                                options->combine->maskVal | pmConfigMask("BLANK", config),
-                                config)) {
-                    psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.\n");
-                    return false;
-                }
-                psFree(roOut->image);
-                roOut->image = NULL;
-            }
-
-            psFree(roOut);              // Drop reference
-
-            if (cellOut->hdu && !cellOut->hdu->blankPHU) {
-                psTrace("ppMerge", 5, "Writing out cell HDU.\n");
-                pmCellWriteMask(cellOut, data->outFile, config->database, false);
-                pmCellFreeData(cellOut);
-            }
-        }
-
-        if (chipOut->hdu && !chipOut->hdu->blankPHU) {
-            psTrace("ppMerge", 5, "Writing out chip HDU.\n");
-            pmChipWriteMask(chipOut, data->outFile, config->database, false, false);
-            pmChipFreeData(chipOut);
-        }
-    }
-
-    // Get list of FPAs for concepts averaging
-    {
-        psList *inFPAs = psListAlloc(NULL); // List of FPAs
-        for (int i = 0; i < filenames->n; i++) {
-            if (! filenames->data[i] || strlen(filenames->data[i]) == 0) {
-                continue;
-            }
-            pmFPA *fpaIn = data->in->data[i]; // Input FPA
-            psListAdd(inFPAs, PS_LIST_TAIL, fpaIn);
-        }
-
-        if (!pmConceptsAverageFPAs(fpaOut, inFPAs)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts.");
-            psFree(inFPAs);
-            return false;
-        }
-        psFree(inFPAs);
-    }
-
-    if (fpaOut->hdu && !fpaOut->hdu->blankPHU) {
-        psTrace("ppMerge", 5, "Writing out FPA HDU.\n");
-        pmFPAWriteMask(fpaOut, data->outFile, config->database, false, false);
-    }
-    pmFPAFreeData(fpaOut);
-
-    psFree(view);
+    for (int i = 0; i < 2; i++) {
+	if (!ppMergeMaskSuspect (data, options, config)) {
+	    psError(PS_ERR_UNKNOWN, true, "failed on mask suspect.\n");
+	    return false;
+	}
+
+	if (!ppMergeMaskBad (data, options, config)) {
+	    psError(PS_ERR_UNKNOWN, true, "failed on mask bad.\n");
+	    return false;
+	}
+    }
+
+    if (!ppMergeMaskAverageConcepts (data, options, config)) {
+	psError(PS_ERR_UNKNOWN, true, "failed on average concepts.\n");
+	return false;
+    }
+
+    if (!ppMergeMaskGrow (data, options, config)) {
+	psError(PS_ERR_UNKNOWN, true, "failed on mask grow.\n");
+	return false;
+    }
+
+    if (!ppMergeMaskWrite (data, options, config)) {
+	psError(PS_ERR_UNKNOWN, true, "failed on mask write.\n");
+	return false;
+    }
 
     return true;
@@ -216,4 +53,5 @@
 
 bool ppMergeMaskReadoutStats(const pmReadout *readout, 
+			     const pmReadout *output, 
 			     ppMergeOptions *options, // Options
 			     psRandom *rng)
@@ -229,6 +67,10 @@
     }
 
+    psImage *mask = NULL;
     psImage *image = readout->image;    // Image of interest
-    psImage *mask = readout->mask;      // Corresponding mask
+
+    if (output) {
+	mask = output->mask;      // Corresponding mask
+    }
 
     if (rng) {
@@ -263,4 +105,5 @@
 
 bool ppMergeMaskChipStats (const pmChip *chip,
+			   const pmChip *output, 
 			   ppMergeOptions *options,
 			   psRandom *rng)
@@ -291,5 +134,12 @@
 
 	    psImage *image = readout->image;    // Image of interest
-	    psImage *mask = readout->mask;      // Corresponding mask
+
+	    pmCell *cellOutput = output->cells->data[nCell];
+	    pmReadout *readoutOutput = NULL;
+	    psImage *mask = NULL;
+	    if (cellOutput->readouts->n > 0) {
+		readoutOutput = cellOutput->readouts->data[0];
+		mask = readoutOutput->mask;      // Corresponding mask
+	    }
 
 	    // Size of image
@@ -318,4 +168,12 @@
     }
 
+    // no valid data, skip the chip
+    if (!values->n) {
+	psFree(values);
+	psFree(stats);
+	psFree(rng);
+	return true;
+    }
+
     // calculate the statistics
     if (!psVectorStats (stats, values, NULL, NULL, 0)) {
@@ -327,5 +185,5 @@
     }
     if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
-	psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
+	psError(PS_ERR_UNKNOWN, false, "invalid image statistics (nan).\n");
 	psFree(values);
 	psFree(stats);
@@ -347,4 +205,6 @@
 	}
     }
+
+    psLogMsg ("ppMerge", PS_LOG_INFO, "statistics for chip: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
 
     psFree(values);