Index: trunk/ppStack/src/ppStackReadout.c
===================================================================
--- trunk/ppStack/src/ppStackReadout.c	(revision 13464)
+++ trunk/ppStack/src/ppStackReadout.c	(revision 13489)
@@ -17,5 +17,11 @@
 bool ppStackReadout(pmConfig *config, const pmFPAview *view)
 {
-    // XXX Somehow need to add an HDU in to the output so that we can actually write it out!
+    // Get the recipe values
+    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
+    float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold
+    float convolveRej = psMetadataLookupF32(NULL, config->arguments, "CONVOLVE.REJ"); // Convolution threshold
+    float extent = psMetadataLookupF32(NULL, config->arguments, "EXTENT"); // Extent of convolution
+    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
+    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
 
     // Get the output target
@@ -31,21 +37,44 @@
                                                             "^PPSTACK.INPUT$"); // Iterator over input files
     psMetadataItem *item;               // Item from iteration
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
+    psRandom *rng = psRandomAlloc(0, PS_RANDOM_TAUS); // Random number generator
     while ((item = psMetadataGetAndIncrement(inputIter))) {
         assert(item->type == PS_DATA_UNKNOWN);
         pmFPAfile *inputFile = item->data.V; // An input file
         pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
+
+        // Brain-dead background subtraction
+        if (!psImageBackground(stats, ro->image, ro->mask, maskBad, rng)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to get image background on image %ld", stack->n);
+            psFree(stats);
+            psFree(rng);
+            psFree(inputIter);
+            psFree(stack);
+            psFree(outRO);
+            return false;
+        }
+        (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32));
+
+        // Normalise each input by the exposure time
+        float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
+        if (!isfinite(exposure)) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "CELL.EXPOSURE is not set for input file %ld", stack->n);
+            psFree(stats);
+            psFree(rng);
+            psFree(inputIter);
+            psFree(outRO);
+            psFree(stack);
+            return false;
+        }
+        (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32));
         pmStackData *data = pmStackDataAlloc(ro, SEEING, WEIGHT); // Data to stack
         psArrayAdd(stack, ARRAY_BUFFER, data);
     }
     psFree(inputIter);
+    psFree(stats);
+    psFree(rng);
 
-    // Get the recipe values
-    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
-    float rej = psMetadataLookupF32(NULL, config->arguments, "REJ"); // Rejection threshold
-    float extent = psMetadataLookupF32(NULL, config->arguments, "EXTENT"); // Extent of convolution
-    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
-    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
-
-    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, iter, rej)) {
+    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, iter, combineRej)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
         psFree(stack);
@@ -54,5 +83,27 @@
     }
 
-    if (!pmStackReject(stack, maskBad, extent, rej)) {
+#if 1
+    {
+        psFits *fits = psFitsOpen("combined.fits", "w");
+        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
+        psFitsClose(fits);
+    }
+
+    for (int i = 0; i < num; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        psImage *inspected = psPixelsToMask(NULL, data->pixels,
+                                            psRegionSet(0, data->sky->image->numCols,
+                                                        0, data->sky->image->numRows),
+                                            maskBlank);
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "inspect%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, inspected, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
+    if (!pmStackReject(stack, maskBad, extent, convolveRej)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to reject input pixels.");
         psFree(stack);
@@ -61,5 +112,21 @@
     }
 
-    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, 0, 0.0)) {
+#if 1
+    for (int i = 0; i < num; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        psImage *rejected = psPixelsToMask(NULL, data->pixels,
+                                           psRegionSet(0, data->sky->image->numCols,
+                                                       0, data->sky->image->numRows),
+                                           maskBlank);
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "reject%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, rejected, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
+    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, 0, combineRej)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
         psFree(stack);
@@ -68,4 +135,8 @@
     }
 
+    outRO->data_exists = true;
+    outRO->parent->data_exists = true;
+    outRO->parent->parent->data_exists = true;
+
     psFree(stack);
     psFree(outRO);                      // Drop reference