Index: trunk/ppStack/src/ppStackReadout.c
===================================================================
--- trunk/ppStack/src/ppStackReadout.c	(revision 21183)
+++ trunk/ppStack/src/ppStackReadout.c	(revision 21477)
@@ -10,6 +10,4 @@
 #include "ppStack.h"
 
-//#define TESTING                  // Write debugging output?
-
 bool ppStackReadoutInitialThread(psThreadJob *job)
 {
@@ -20,10 +18,9 @@
     pmConfig *config = args->data[1];   // Configuration
     pmReadout *outRO = args->data[2];   // Output readout
-    psArray *subRegions = args->data[3]; // Regions for PSF-matching
-    psArray *subKernels = args->data[4]; // Kernels for PSF-matching
-    psVector *weightings = args->data[5]; // Weightings (1/noise^2) for each image
-    psVector *addVariance = args->data[6]; // Additional variance when rejecting
-
-    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels,
+    psVector *mask = args->data[3];     // Mask for inputs
+    psVector *weightings = args->data[4]; // Weightings (1/noise^2) for each image
+    psVector *addVariance = args->data[5]; // Additional variance when rejecting
+
+    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
                                              weightings, addVariance);
 
@@ -42,9 +39,11 @@
     pmConfig *config = args->data[1];   // Configuration
     pmReadout *outRO = args->data[2];   // Output readout
-    psArray *rejected = args->data[3];  // Rejected pixels
-    psVector *weightings = args->data[4];  // Weighting (1/noise^2) for each image
-
-    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected,
-                                      weightings); // Status of operation
+    psVector *mask = args->data[3];     // Mask for inputs
+    psArray *rejected = args->data[4];  // Rejected pixels
+    psVector *weightings = args->data[5];  // Weighting (1/noise^2) for each image
+    psVector *addVariance = args->data[6];  // Weighting (1/noise^2) for each image
+
+    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
+                                      weightings, addVariance); // Status of operation
 
     thread->busy = false;
@@ -85,14 +84,10 @@
 
 psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
-                               const psArray *regions, const psArray *kernels, const psVector *weightings,
-                               const psVector *addVariance)
+                               const psVector *mask, const psVector *weightings, const psVector *addVariance)
 {
     assert(config);
     assert(outRO);
     assert(readouts);
-    assert(regions);
-    assert(kernels);
-    assert(readouts->n == regions->n);
-    assert(regions->n == kernels->n);
+    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
     assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
     assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
@@ -107,4 +102,5 @@
     float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
     float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
+    float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction
     bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
     bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
@@ -123,5 +119,5 @@
     for (int i = 0; i < num; i++) {
         pmReadout *ro = readouts->data[i];
-        if (!ro) {
+        if (!ro || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
             // Bad image
             continue;
@@ -148,21 +144,10 @@
     }
 
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, combineRej, combineSys,
-                        true, useVariance, safe)) {
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
+                        combineRej, combineSys, combineDiscard, true, useVariance, safe)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
         psFree(stack);
         return false;
     }
-
-#ifdef TESTING
-    {
-        psString name = NULL;           // Name of image
-        psStringAppend(&name, "combined_initial_%03d.fits", sectionNum);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
-        psFitsClose(fits);
-    }
-#endif
 
     // Save list of pixels to inspect
@@ -189,5 +174,6 @@
 
 bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
-                         const psArray *rejected, const psVector *weightings)
+                         const psVector *mask, const psArray *rejected, const psVector *weightings,
+                         const psVector *addVariance)
 {
     assert(config);
@@ -196,4 +182,5 @@
     assert(rejected);
     assert(readouts->n == rejected->n);
+    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
     assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
 
@@ -217,12 +204,10 @@
     int numGood = num;                  // Number of good inputs: images that haven't been completely rejected
     for (int i = 0; i < num; i++) {
-        if (!rejected->data[i]) {
+        pmReadout *ro = readouts->data[i];
+        if (!ro || !rejected->data[i] || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
             // Image completely rejected
             numGood--;
             continue;
         }
-
-        pmReadout *ro = readouts->data[i];
-        assert(ro);
 
 #if 0
@@ -242,10 +227,10 @@
         }
 
-        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], NAN);
+        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
         data->reject = psMemIncrRefCounter(rejected->data[i]);
         stack->data[i] = data;
     }
 
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN,
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN, NAN,
                         numGood != num, useVariance, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
@@ -254,15 +239,4 @@
     }
 
-#ifdef TESTING
-    {
-        psString name = NULL;           // Name of image
-        psStringAppend(&name, "combined_final_%03d.fits", sectionNum);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
-        psFitsClose(fits);
-    }
-#endif
-
     pmCell *outCell = outRO->parent;    // Output cell
     pmChip *outChip = outCell->parent;  // Output chip