Index: trunk/ppStack/src/ppStackCombineFinal.c
===================================================================
--- trunk/ppStack/src/ppStackCombineFinal.c	(revision 25096)
+++ trunk/ppStack/src/ppStackCombineFinal.c	(revision 26076)
@@ -10,9 +10,16 @@
 #include "ppStackLoop.h"
 
-bool ppStackCombineFinal(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
+//#define TESTING                         // Enable test output
+
+bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances,
+                         ppStackOptions *options, pmConfig *config, bool safe, bool normalise)
 {
     psAssert(stack, "Require stack");
     psAssert(options, "Require options");
     psAssert(config, "Require configuration");
+
+    if (!target->mask) {
+        target->mask = psImageAlloc(target->image->numCols, target->image->numRows, PS_TYPE_IMAGE_MASK);
+    }
 
     stack->lastScan = 0;            // Reset read
@@ -30,9 +37,12 @@
         }
 
-        // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
+        // call: ppStackReadoutFinal(config, target, readouts, rejected)
         psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
+        psArrayAdd(job->args, 1, target);
         psArrayAdd(job->args, 1, thread);
         psArrayAdd(job->args, 1, options);
         psArrayAdd(job->args, 1, config);
+        PS_ARRAY_ADD_SCALAR(job->args, safe, PS_TYPE_U8);
+        PS_ARRAY_ADD_SCALAR(job->args, normalise, PS_TYPE_U8);
         if (!psThreadJobAddPending(job)) {
             psFree(job);
@@ -48,31 +58,43 @@
 
     // Sum covariance matrices
-    double sumWeights = 0.0;            // Sum of weights
-    for (int i = 0; i < options->num; i++) {
-        if (options->inputMask->data.U8[i]) {
-            psFree(options->covariances->data[i]);
-            options->covariances->data[i] = NULL;
-            continue;
+    if (covariances) {
+        double sumWeights = 0.0;            // Sum of weights
+        for (int i = 0; i < options->num; i++) {
+            if (options->inputMask->data.U8[i]) {
+                psFree(covariances->data[i]);
+                covariances->data[i] = NULL;
+                continue;
+            }
+            psKernel *covar = covariances->data[i]; // Covariance matrix
+            if (!covar) {
+                continue;
+            }
+            float weight = options->weightings->data.F32[i]; // Weight to apply
+            psBinaryOp(covar->image, covar->image, "*", psScalarAlloc(weight, PS_TYPE_F32));
+            sumWeights += weight;
         }
-        psKernel *covar = options->covariances->data[i]; // Covariance matrix
-        if (!covar) {
-            continue;
+        if (sumWeights > 0.0) {
+            target->covariance = psImageCovarianceSum(covariances);
+            psBinaryOp(target->covariance->image, target->covariance->image, "/",
+                       psScalarAlloc(sumWeights, PS_TYPE_F32));
+            psImageCovarianceTransfer(target->variance, target->covariance);
         }
-        float weight = options->weightings->data.F32[i]; // Weight to apply
-        psBinaryOp(covar->image, covar->image, "*", psScalarAlloc(weight, PS_TYPE_F32));
-        sumWeights += weight;
-    }
-    if (sumWeights > 0.0) {
-        pmReadout *outRO = options->outRO;  // Output readout
-        outRO->covariance = psImageCovarianceSum(options->covariances);
-        psBinaryOp(outRO->covariance->image, outRO->covariance->image, "/",
-                   psScalarAlloc(sumWeights, PS_TYPE_F32));
-        psFree(options->covariances); options->covariances = NULL;
-        psImageCovarianceTransfer(outRO->variance, outRO->covariance);
+    } else {
+        target->covariance = psImageCovarianceNone();
     }
 
 #ifdef TESTING
-    pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits");
-    ppStackWriteImage("combined_final.fits", NULL, outRO->image, config);
+    static int pass = 0;                // Pass through
+    psString name = NULL;               // Name of file
+    psStringAppend(&name, "combined_image_final_%d.fits", pass);
+    pass++;
+    ppStackWriteImage(name, NULL, target->image, config);
+    psStringSubstitute(&name, "mask", "image");
+    ppStackWriteImage(name, NULL, target->mask, config);
+    psStringSubstitute(&name, "variance", "mask");
+    ppStackWriteImage(name, NULL, target->variance, config);
+    psFree(name);
+
+    pmStackVisualPlotTestImage(target->image, "combined_image_final.fits");
 #endif