Index: trunk/ppStack/src/ppStackReadout.c
===================================================================
--- trunk/ppStack/src/ppStackReadout.c	(revision 23577)
+++ trunk/ppStack/src/ppStackReadout.c	(revision 26076)
@@ -25,11 +25,9 @@
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
 
-    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
-                                             weightings, addVariance);
-
-    job->results = inspect;
+    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
+                                         weightings, addVariance);
     thread->busy = false;
 
-    return inspect ? true : false;
+    return job->results ? true : false;
 }
 
@@ -39,16 +37,19 @@
 
     psArray *args = job->args;          // Arguments
-    ppStackThread *thread = args->data[0]; // Thread
-    ppStackOptions *options = args->data[1]; // Options
-    pmConfig *config = args->data[2];   // Configuration
-
-    pmReadout *outRO = options->outRO;  // Output readout
+    pmReadout *target = args->data[0];  // Output readout
+    ppStackThread *thread = args->data[1]; // Thread
+    ppStackOptions *options = args->data[2]; // Options
+    pmConfig *config = args->data[3];   // Configuration
+    bool safety = PS_SCALAR_VALUE(args->data[4], U8);    // Safety switch on?
+    bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images?
+
     psVector *mask = options->inputMask; // Mask for inputs
     psArray *rejected = options->rejected; // Rejected pixels
     psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
-
-    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
-                                      weightings, addVariance); // Status of operation
+    psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
+
+    bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, rejected,
+                                      weightings, addVariance, safety, norm); // Status of operation
 
     thread->busy = false;
@@ -63,24 +64,35 @@
 
     psArray *args = job->args;  // Input arguments
-    psArray *inspect = args->data[0]; // Array of pixel arrays
-    int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest
-
-    psArray *inputs = inspect->data[index]; // Array of interest
-    psPixels *output = NULL;    // Output pixel list
-    for (int i = 0; i < inputs->n; i++) {
-        psPixels *input = inputs->data[i]; // Input pixel list
-        if (!input || input->n == 0) {
-            continue;
-        }
-        output = psPixelsConcatenate(output, input);
-    }
-
-    if (!output) {
-        // If there are no pixels to inspect, then just fake it
-        output = psPixelsAllocEmpty(0);
-    }
-
-    psFree(inputs);
-    inspect->data[index] = output;
+    psArray *inspects = args->data[0]; // Array of pixel arrays
+    psArray *rejects = args->data[1];  // Array of pixel arrays
+    int index = PS_SCALAR_VALUE(args->data[2], S32); // Index of interest
+
+    psArray *inInspects = inspects->data[index]; // Array of interest
+    psArray *inRejects = rejects->data[index]; // Array of interest
+    psAssert(inInspects->n == inRejects->n, "Size should be the same");
+    psPixels *outInspect = NULL, *outReject = NULL; // Output pixel lists
+    for (int i = 0; i < inInspects->n; i++) {
+        psPixels *inInspect = inInspects->data[i]; // Input pixel list
+        if (inInspect && inInspect->n > 0) {
+            outInspect = psPixelsConcatenate(outInspect, inInspect);
+        }
+        psPixels *inReject = inRejects->data[i]; // Input pixel list
+        if (inReject && inReject->n > 0) {
+            outReject = psPixelsConcatenate(outReject, inReject);
+        }
+    }
+
+    // If there are no pixels to inspect, then just fake it
+    if (!outInspect) {
+        outInspect = psPixelsAllocEmpty(0);
+    }
+    if (!outReject) {
+        outReject = psPixelsAllocEmpty(0);
+    }
+
+    psFree(inspects->data[index]);
+    inspects->data[index] = outInspect;
+    psFree(rejects->data[index]);
+    rejects->data[index] = outReject;
 
     return true;
@@ -104,5 +116,5 @@
 
     bool mdok;                          // Status of MD lookup
-    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
+    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
     float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
     float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
@@ -114,6 +126,8 @@
     int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
@@ -149,6 +163,6 @@
     }
 
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
-                        combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) {
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
+                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
         psFree(stack);
@@ -156,6 +170,7 @@
     }
 
-    // Save list of pixels to inspect
+    // Save lists of pixels
     psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
+    psArray *reject = psArrayAlloc(num);  // List of pixels rejected
     for (int i = 0; i < num; i++) {
         pmStackData *data = stack->data[i]; // Data for this image
@@ -168,4 +183,5 @@
         }
         inspect->data[i] = psMemIncrRefCounter(data->inspect);
+        reject->data[i] = psMemIncrRefCounter(data->reject);
     }
     psFree(stack);
@@ -173,5 +189,9 @@
     sectionNum++;
 
-    return inspect;
+    psArray *results = psArrayAlloc(2); // Array of results
+    results->data[0] = inspect;
+    results->data[1] = reject;
+
+    return results;
 }
 
@@ -180,5 +200,5 @@
 bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
                          const psVector *mask, const psArray *rejected, const psVector *weightings,
-                         const psVector *addVariance)
+                         const psVector *addVariance, bool safety, const psVector *norm)
 {
     assert(config);
@@ -196,13 +216,11 @@
 
     bool mdok;                          // Status of MD lookup
-    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
-    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
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
@@ -211,22 +229,15 @@
     psArray *stack = psArrayAlloc(num); // Array for stacking
 
-    bool entire = true;                 // Combine entire image?
-    if (rejected) {
-        // We have rejection from a previous combination: combine without flagging pixels to inspect
-        entire = false;
-        safe = false;
-        iter = 0;
-        combineRej = NAN;
-        combineSys = NAN;
-    }
+    // We have rejection from a previous combination: combine without flagging pixels to inspect
+    safe &= safety;
+    int iter = 0;
+    float combineRej = NAN;
+    float combineSys = NAN;
+    float combineDiscard = NAN;
 
     for (int i = 0; i < num; i++) {
         pmReadout *ro = readouts->data[i];
-        if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) {
-            // Image completely rejected since previous combination
-            entire = true;
-            continue;
-        } else if (mask->data.U8[i]) {
-            // Image completely rejected before original combination
+        if (mask->data.U8[i]) {
+            // Image completely rejected
             continue;
         }
@@ -252,9 +263,14 @@
         data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
         stack->data[i] = data;
-    }
-
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0,
-                        iter, combineRej, combineSys, combineDiscard,
-                        entire, useVariance, safe, !rejected)) {
+
+        if (norm) {
+            float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
+            psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(normalise, PS_TYPE_F32));
+            psBinaryOp(ro->variance, ro->variance, "*", psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32));
+        }
+    }
+
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
+                        combineSys, combineDiscard, useVariance, safe, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
         psFree(stack);