Index: branches/tap_branches/ppStack/src/ppStackCombineFinal.c
===================================================================
--- branches/tap_branches/ppStack/src/ppStackCombineFinal.c	(revision 25900)
+++ branches/tap_branches/ppStack/src/ppStackCombineFinal.c	(revision 27838)
@@ -10,9 +10,48 @@
 #include "ppStackLoop.h"
 
-bool ppStackCombineFinal(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
+
+//#define TESTING                         // Enable test output
+
+bool ppStackCombineFinal(ppStackThreadData *stack, psArray *covariances, ppStackOptions *options,
+                         pmConfig *config, bool safe, bool normalise, bool grow)
 {
     psAssert(stack, "Require stack");
     psAssert(options, "Require options");
     psAssert(config, "Require configuration");
+
+    pmReadout *outRO = options->outRO;                                      // Output readout
+    pmReadout *expRO = options->expRO;                                      // Exposure readout
+    int numCols = outRO->image->numCols, numRows = outRO->image->numRows; // Size of image
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+    float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor"
+
+    // Grow the list of rejected pixels, if desired
+    psArray *reject = psArrayAlloc(options->num); // Pixels rejected for each image
+    for (int i = 0; i < options->num; i++) {
+        if (options->inputMask->data.U8[i]) {
+            continue;
+        }
+        if (grow) {
+            reject->data[i] = pmStackRejectGrow(options->rejected->data[i], numCols, numRows, poorFrac,
+                                                options->regions->data[i], options->kernels->data[i]);
+            if (!reject->data[i]) {
+                psError(psErrorCodeLast(), false, "Unable to grow rejected pixels for image %d", i);
+                psFree(reject);
+                return false;
+            }
+        } else {
+            reject->data[i] = psMemIncrRefCounter(options->rejected->data[i]);
+        }
+    }
+
+    if (!outRO->mask) {
+        outRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+    }
+    if (!expRO->mask) {
+        expRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+    }
 
     stack->lastScan = 0;            // Reset read
@@ -22,5 +61,6 @@
         if (!status) {
             // Something went wrong
-            psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
+            psError(psErrorCodeLast(), false, "Unable to read chunk %d", numChunk);
+            psFree(reject);
             return false;
         }
@@ -33,8 +73,12 @@
         psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
         psArrayAdd(job->args, 1, thread);
+        psArrayAdd(job->args, 1, reject);
         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);
+            psFree(reject);
             return false;
         }
@@ -43,36 +87,100 @@
 
     if (!psThreadPoolWait(true)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to do final combination.");
+        psError(psErrorCodeLast(), false, "Unable to do final combination.");
+        psFree(reject);
         return false;
     }
 
+    psFree(reject);
+
     // 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;
+    if (covariances) {
+        outRO->covariance = psImageCovarianceAverageWeighted(covariances, options->weightings);
+    } else {
+        outRO->covariance = psImageCovarianceNone();
+    }
+
+    // Propagate WCS
+    bool wcsDone = false;           // Have we done the WCS?
+    for (int i = 0; i < options->num && !wcsDone; i++) {
+        if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
             continue;
         }
-        psKernel *covar = options->covariances->data[i]; // Covariance matrix
-        if (!covar) {
-            continue;
+
+        ppStackThread *thread = stack->threads->data[0]; // Representative stack
+        pmReadout *inRO = thread->readouts->data[i]; // Template readout
+        if (inRO && !wcsDone) {
+            // Copy astrometry over
+            wcsDone = true;
+            pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
+            pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
+            pmChip *outChip = outRO->parent->parent; // Output chip
+            pmFPA *outFPA = outChip->parent; // Output FPA
+            if (!outHDU || !inHDU) {
+                psWarning("Unable to find HDU at FPA level to copy astrometry.");
+            } else {
+                if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
+                    psErrorClear();
+                    psWarning("Unable to read WCS astrometry from input FPA.");
+                    wcsDone = false;
+                } else {
+                    if (!outHDU->header) {
+                        outHDU->header = psMetadataAlloc();
+                    }
+                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
+                        psErrorClear();
+                        psWarning("Unable to write WCS astrometry to output FPA.");
+                        wcsDone = false;
+                    }
+                }
+            }
         }
-        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);
     }
 
+    // Set exposure time correctly
+    {
+        float exptime = 0.0;            // Summed exposure time
+        for (int i = 0; i < options->num; i++) {
+            if (options->inputMask->data.U8[i]) {
+                continue;
+            }
+            exptime += options->exposures->data.F32[i];
+        }
+
+        {
+            psMetadataItem *item = psMetadataLookup(outRO->parent->concepts, "CELL.EXPOSURE");
+            item->data.F32 = exptime;
+        }
+        {
+            psMetadataItem *item = psMetadataLookup(outRO->parent->parent->parent->concepts, "FPA.EXPOSURE");
+            item->data.F32 = exptime;
+        }
+    }
+
+    // Put version information into the header
+    pmHDU *hdu = pmHDUFromCell(outRO->parent);
+    if (!hdu) {
+        psError(PPSTACK_ERR_PROG, false, "Unable to find HDU for output.");
+        return false;
+    }
+    if (!hdu->header) {
+        hdu->header = psMetadataAlloc();
+    }
+    ppStackVersionHeader(hdu->header);
+
+
 #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, outRO->image, config);
+    psStringSubstitute(&name, "mask", "image");
+    ppStackWriteImage(name, NULL, outRO->mask, config);
+    psStringSubstitute(&name, "variance", "mask");
+    ppStackWriteImage(name, NULL, outRO->variance, config);
+    psFree(name);
+
+    pmStackVisualPlotTestImage(outRO->image, "combined_image_final.fits");
 #endif