Index: /trunk/ppStack/src/ppStackArguments.c
===================================================================
--- /trunk/ppStack/src/ppStackArguments.c	(revision 13488)
+++ /trunk/ppStack/src/ppStackArguments.c	(revision 13489)
@@ -80,5 +80,6 @@
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stat", 0, "Statistics file", NULL);
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
-    psMetadataAddF32(arguments, PS_LIST_TAIL, "-rej", 0, "Rejection thresold (sigma)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0, "Combination rejection thresold (sigma)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-convolve-rej", 0, "Convolution rejection thresold (sigma)", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-extent", 0, "Extent of convolution (sigma)", NAN);
     psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-bad", 0, "Mask value for bad pixels", 0);
@@ -89,5 +90,5 @@
     }
 
-    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "IMAGE.LIST", 0,
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "IMAGES.LIST", 0,
                      "Name of the input image list", argv[1]);
     psMetadataAddStr(config->arguments, PS_LIST_TAIL, "MASKS.LIST", 0,
@@ -104,9 +105,10 @@
     }
 
-    VALUE_ARG_RECIPE_INT("-iter",       "ITER",       S32, 0);
-    VALUE_ARG_RECIPE_FLOAT("-rej",      "REJ",        F32);
-    VALUE_ARG_RECIPE_FLOAT("-extent",   "EXTENT",     F32);
-    VALUE_ARG_RECIPE_INT("-mask-bad",   "MASK.BAD",   U8, 0);
-    VALUE_ARG_RECIPE_INT("-mask-blank", "MASK.BLANK", U8, 0);
+    VALUE_ARG_RECIPE_INT("-iter",           "ITER",         S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-combine-rej",  "COMBINE.REJ",  F32);
+    VALUE_ARG_RECIPE_FLOAT("-convolve-rej", "CONVOLVE.REJ", F32);
+    VALUE_ARG_RECIPE_FLOAT("-extent",       "EXTENT",       F32);
+    VALUE_ARG_RECIPE_INT("-mask-bad",       "MASK.BAD",     U8, 0);
+    VALUE_ARG_RECIPE_INT("-mask-blank",     "MASK.BLANK",   U8, 0);
 
     psTrace("ppStack", 1, "Done reading command-line arguments\n");
Index: /trunk/ppStack/src/ppStackCamera.c
===================================================================
--- /trunk/ppStack/src/ppStackCamera.c	(revision 13488)
+++ /trunk/ppStack/src/ppStackCamera.c	(revision 13489)
@@ -104,5 +104,11 @@
 
     // Output image
-    pmFPAfile *output = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT");
+    pmFPA *fpa = pmFPAConstruct(config->camera); // FPA to contain the output
+    if (!fpa) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
+        return false;
+    }
+    pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSTACK.OUTPUT");
+    psFree(fpa);                        // Drop reference
     if (!output) {
         psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT"));
@@ -113,4 +119,41 @@
         return false;
     }
+
+    // Generate the correct structure
+    pmFPALevel phuLevel = pmFPAPHULevel(output->format); // Level at which PHU goes
+    pmFPAview *view = pmFPAviewAlloc(0);// View for current level
+    if (phuLevel == PM_FPA_LEVEL_FPA) {
+        if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
+            psFree(fpa);
+            psFree(view);
+            return NULL;
+        }
+    } else {
+        pmChip *chip;                       // Chip from FPA
+        while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
+            if (phuLevel == PM_FPA_LEVEL_CHIP) {
+                if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
+                    psFree(fpa);
+                    psFree(view);
+                    return NULL;
+                }
+            } else {
+                pmCell *cell;                   // Cell from chip
+                while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
+                    if (phuLevel == PM_FPA_LEVEL_CELL) {
+                        if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
+                            psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
+                            psFree(fpa);
+                            psFree(view);
+                            return NULL;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    psFree(view);
 
     // Output mask
Index: /trunk/ppStack/src/ppStackLoop.c
===================================================================
--- /trunk/ppStack/src/ppStackLoop.c	(revision 13488)
+++ /trunk/ppStack/src/ppStackLoop.c	(revision 13489)
@@ -30,7 +30,7 @@
     }
 
-    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
-    if (!output) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n");
+    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Token input file
+    if (!input) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n");
         return false;
     }
@@ -45,5 +45,5 @@
 
     pmChip *chip;                       // Chip of interest
-    while ((chip = pmFPAviewNextChip(view, output->fpa, 1)) != NULL) {
+    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
             return false;
@@ -51,5 +51,5 @@
 
         pmCell *cell;                // Cell of interest
-        while ((cell = pmFPAviewNextCell(view, output->fpa, 1)) != NULL) {
+        while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
             if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                 return false;
@@ -67,5 +67,5 @@
 
             pmReadout *readout;         // Readout of interest
-            while ((readout = pmFPAviewNextReadout(view, output->fpa, 1))) {
+            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1))) {
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                     return false;
@@ -85,5 +85,5 @@
             // Perform statistics on the cell
             if (stats) {
-                ppStats(stats, output->fpa, view, config);
+                ppStats(stats, input->fpa, view, config);
             }
 
Index: /trunk/ppStack/src/ppStackReadout.c
===================================================================
--- /trunk/ppStack/src/ppStackReadout.c	(revision 13488)
+++ /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