Changeset 16382

Timestamp:

Feb 8, 2008, 11:43:52 AM (18 years ago)

Author:

Paul Price

Message:

Code compiles and is ready for preliminary testing. Some features (e.g., rejection based on initial combination) are disabled for now, marked with XXX in the code.

Location:

branches/pap_branch_080207/ppStack/src

Files:

• Makefile.am (modified) (1 diff)
• ppStack.h (modified) (2 diffs)
• ppStackLoop.c (modified) (15 diffs)
• ppStackPhotometry.c (modified) (1 diff)
• ppStackReadout.c (modified) (9 diffs)

branches/pap_branch_080207/ppStack/src/Makefile.am

@@ -9 +9 @@
         ppStackCamera.c         \
         ppStackLoop.c           \
+        ppStackPSF.c            \
         ppStackReadout.c        \
+        ppStackPhotometry.c     \
         ppStackVersion.c        \
         ppStackMatch.c

branches/pap_branch_080207/ppStack/src/ppStack.h

@@ -19 +19 @@
     );
 
+// Determine target PSF for input images
+pmPSF *ppStackPSF(const pmConfig *config, // Configuration
+                  int numCols, int numRows, // Size of image
+                  const psList *list    // List of input PSFs
+    );
+
 // Perform stacking on a readout
-bool ppStackReadout(pmConfig *config,   // Configuration
-                    const pmFPAview *view // View for readout
+bool ppStackReadout(const pmConfig *config,   // Configuration
+                    pmReadout *outRO,   // Output readout
+                    const psArray *readouts // Input readouts
+    );
+
+// Perform photometry on stack
+bool ppStackPhotometry(pmConfig *config, // Configuration
+                       const pmReadout *readout, // Readout to be photometered
+                       const pmFPAview *view // View to readout
     );
 
@@ -35 +48 @@
 
 /// Convolve image to match specified seeing
-bool ppStackMatch(pmReadout *output,    ///< Convolved readout
-                  const pmReadout *input, ///< Readout to be convolved
+bool ppStackMatch(pmReadout *readout, ///< Readout to be convolved; replaced with output
                   const pmReadout *sourcesRO, ///< Readout with sources
                   const pmPSF *psf,     ///< Target PSF

branches/pap_branch_080207/ppStack/src/ppStackLoop.c

@@ -13 +13 @@
 // Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
 
+#if 0
 // All files in the system
 static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
@@ -23 +24 @@
                             "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", "PSPHOT.INPUT.CMF",
                             0 };
+#endif
 
 // Files required in preparation for convolution
 static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 };
 
-// Files for convolution; these should be turned on one by one
-static char *convolutionFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
-                                    "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT",
-                                    0 };
-
 // Input files for the combination
-static char *inCombineFiles = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 };
+static char *inCombineFiles[] = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 };
 
 // Output files for the combination and photometry
-static char *outCombineFiles = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
-                                 "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
-                                 "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
-                                 "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
-                                 "PSPHOT.INPUT.CMF", 0 };
+static char *outCombineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
+                                   "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
+                                   "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
+                                   "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
+                                   "PSPHOT.INPUT.CMF", 0 };
 
 
@@ -103 +100 @@
 
 
-
 // Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
 static pmFPAview *filesIterateDown(pmConfig *config // Configuration
@@ -158 +154 @@
 bool ppStackLoop(pmConfig *config)
 {
+    assert(config);
+
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
 
@@ -182 +180 @@
         return false;
     }
+    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    int numScans = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of scans to read at once
+
+
+
+    // XXX Will need to adjust this when we're doing the convolution
+    int overlap = 0;                    // Overlap between consecutive scans
+
+
+
 
     // Preparation iteration: Load the sources, and get a target PSF model
-    psArray *sources = NULL;            // Sources to use for PSF matching
+    pmReadout *sources = NULL;          // Readout with sources to use for PSF matching
     pmPSF *targetPSF = NULL;            // Target PSF
     {
@@ -194 +202 @@
         }
 
-        pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Readout
         // We want to hang on to the 'sources' even when its host FPA is blown away
-        sources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"));
+        sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"));
         if (!sources) {
             psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
@@ -207 +214 @@
                                                                "^PPSTACK.INPUT$"); // Iterator
         psMetadataItem *fileItem; // Item from iteration
-        int fileNum = 0;        // Number of file
         psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
+        int numCols = 0, numRows = 0;   // Size of image
         while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
             assert(fileItem->type == PS_DATA_UNKNOWN);
             pmFPAfile *inputFile = fileItem->data.V; // An input file
             pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
-            pmPSF *psf = psMetadataLookupPtr(NULL, ro->parent->parent->analysis, "PSPHOT.PSF"); // PSF
+            pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
             if (!psf) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
                 psFree(view);
                 psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
                 return false;
             }
             psListAdd(psfList, PS_LIST_TAIL, psf);
+
+            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
+            pmHDU *hdu = pmHDUFromCell(cell);
+            assert(hdu && hdu->header);
+            int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns
+            int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
+            if (naxis1 <= 0 || naxis2 <= 0) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
+                psFree(view);
+                psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
+                return false;
+            }
+            if (numCols == 0 && numRows == 0) {
+                numCols = naxis1;
+                numRows = naxis2;
+            }
         }
         psFree(fileIter);
         psFree(view);
 
-        psArray *psfArray = psListToArray(psfList); // Array of PSFs
+        targetPSF = ppStackPSF(config, numCols, numRows, psfList);
         psFree(psfList);
-        targetPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel,
-                                  psfOrder, psfOrder);
-        psFree(psfArray);
         if (!targetPSF) {
             psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
@@ -241 +265 @@
 
     // Generate convolutions and write them to disk
-    psArray *cells = psArrayAlloc(numFiles);    // Cells for convolved images --- a handle for reading again
-    for (int i = 0; i < numFiles; i++) {
+    psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
+    for (int i = 0; i < num; i++) {
         pmFPAfileActivate(config->files, false, NULL);
         psArray *files = fileActivationSingle(config, prepareFiles, true, i);
         pmFPAview *view = filesIterateDown(config);
         if (!view) {
+            psFree(sources);
+            psFree(targetPSF);
             return false;
         }
@@ -254 +280 @@
 
         // Background subtraction, scaling and normalisation is performed automatically by the image matching
-        if (!ppStackMatch(ro, sources, outPSF, config)) {
+        if (!ppStackMatch(readout, sources, targetPSF, config)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i);
-            psFree(convolved);
+            psFree(sources);
+            psFree(targetPSF);
             return false;
         }
@@ -290 +317 @@
         filesIterateUp(config);
     }
+    psFree(sources);
+    psFree(targetPSF);
 
     // Set files to read at the readout level --- then when we iterate down, the pmFPAfileIO stuff should take
@@ -295 +324 @@
     fileSetDataLevel(config, inCombineFiles, PM_FPA_LEVEL_CHUNK);
     pmFPAfileActivate(config->files, false, NULL);
-    fileActivation(inCombineFiles, true);
-    fileActivation(outCombineFiles, true);
+    fileActivation(config, inCombineFiles, true);
+    fileActivation(config, outCombineFiles, true);
 
     {
@@ -304 +333 @@
             return false;
         }
-        pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.OUTPUT"); // Output readout
-        psFree(view);
-
-        psArray *readouts = psArrayAlloc(numFiles); // Readouts for convolved images
-        for (int i = 0; i < numFiles; i++) {
+        pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
+        pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+
+        psArray *readouts = psArrayAlloc(num); // Readouts for convolved images
+        for (int i = 0; i < num; i++) {
             readouts->data[i] = pmReadoutAlloc(cells->data[i]); // Readout into which to read
         }
+        psFree(cells);
 
         // Read convolutions by chunks
-        int numRead = 0;                    // Number of inputs read
-        int numChunk = 0;                   // Chunk number
         bool more = true;               // More to read?
         for (int numChunk = 0; more; numChunk++) {
-            for (int i = 0; i < numFiles; i++) {
+            for (int i = 0; i < num; i++) {
                 pmReadout *readout = readouts->data[i];
                 assert(readout);
@@ -336 +364 @@
                     psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
                     psFree(readouts);
-                    psFree(cells);
+                    psFree(outRO);
+                    psFree(view);
+                    return false;
                 }
             }
 
+#if 0
             if (!ppStackReadout(config, outRO, readouts)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
+                psFree(readouts);
+                psFree(outRO);
+                psFree(view);
                 return false;
             }
-
-            for (int i = 0; i < numFiles && more; i++) {
+#else
+            printf("Stack...\n");
+#endif
+
+            for (int i = 0; i < num && more; i++) {
                 pmReadout *readout = readouts->data[i];
                 assert(readout);
@@ -364 +401 @@
             }
         }
+        psFree(readouts);
 
         // Get rid of the input files
-        fileActivation(outCombineFiles, false);
+        fileActivation(config, outCombineFiles, false);
         filesIterateUp(config);
 
-        ppStackPhotometry(config, outRO);
+        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY") &&
+            !ppStackPhotometry(config, outRO, view)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
+            psFree(outRO);
+            psFree(view);
+            return false;
+        }
+
+        // Statistics on output
+        if (stats) {
+            ppStatsFPA(stats, outCell->parent->parent, view, maskBlank, config);
+        }
 
         // Put version information into the header
         pmHDU *hdu = pmHDUFromCell(outRO->parent);
         if (!hdu) {
-            psError(PS_ERR_UNKNOWN, true, "Unable to find HDU for output.");
+            psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
             return false;
         }
@@ -382 +431 @@
         ppStackVersionMetadata(hdu->header);
 
-        // Statistics on output
-        if (stats) {
-            ppStatsFPA(stats, output->fpa, view, maskBlank, config);
-        }
-
         // Write out the output files
-        filesActivation(outCombineFiles, true);
+        fileActivation(config, outCombineFiles, true);
         filesIterateUp(config);
-    }
-
-    psFree(view);
+
+        psFree(outRO);
+        psFree(view);
+    }
+
 
     // Write out summary statistics

branches/pap_branch_080207/ppStack/src/ppStackPhotometry.c

@@ -10 +10 @@
 #include "ppStack.h"
 
-bool ppStackPhotometry(const pmConfig *config, const pmReadout *readout, const pmFPAview *view)
+bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view)
 {
     pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-    pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
+    pmFPACopy(photFile->fpa, readout->parent->parent->parent);
 
     if (!psphotReadout(config, view)) {

branches/pap_branch_080207/ppStack/src/ppStackReadout.c

@@ -10 +10 @@
 #include "ppStack.h"
 
-#define ARRAY_BUFFER 16                 // Number to add to array at a time
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
-//#define NO_CONVOLUTION                  // Don't perform convolutions?
-//#define CONVOLUTION_FILES               // Write convolutions?
 //#define REJECTION_FILES                 // Write rejection mask?
 //#define INSPECTION_FILES                // Write inspection mask?
 
-bool ppStackReadout(pmConfig *config, const pmFPAview *view)
+bool ppStackReadout(const pmConfig *config, pmReadout *outRO, const psArray *readouts)
 {
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+
     // Get the recipe values
-    bool mdok;                          // Status of MD lookup
     int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
     float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold
     psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
-    float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
-    bool photometry = psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY"); // Perform photometry?
-    int psfInstances = psMetadataLookupS32(&mdok, config->arguments, "PSF.INSTANCES"); // Number of instances for PSF
-    float psfRadius = psMetadataLookupF32(NULL, config->arguments, "PSF.RADIUS"); // Radius for PSF
-    const char *psfModel = psMetadataLookupStr(NULL, config->arguments, "PSF.MODEL"); // Model for PSF
-    int psfOrder = psMetadataLookupS32(&mdok, config->arguments, "PSF.ORDER"); // Spatial order for PSF
-
-    // Get the output target
-    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
-    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
-    pmFPA *outFPA = outCell->parent->parent; // Output FPA
-
-    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-    assert(num > 0);
-
-    // Get the input sources
-    psArray *stack = psArrayAllocEmpty(ARRAY_BUFFER); // The stack of inputs
-    pmReadout *sources = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Sources for stamps
-    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
-                                                           "^PPSTACK.INPUT$"); // Iterator over input files
-    psMetadataItem *fileItem;           // Item from iteration
-    int fileNum = 0;                    // Number of file
-    psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
-    pmPSF *outPSF = NULL;               // Ouptut PSF
-    int numCols = 0, numRows = 0;       // Size of image
-    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
-        assert(fileItem->type == PS_DATA_UNKNOWN);
-        pmFPAfile *inputFile = fileItem->data.V; // An input file
-        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
-        pmCell *cell = ro->parent;      // The cell
-        pmChip *chip = cell->parent;    // The chip: holds the PSF
-
-        bool mdok;                      // Status of MD lookup
-        pmPSF *psf = psMetadataLookupPtr(&mdok, chip->analysis, "PSPHOT.PSF");
-        if (mdok && psf) {
-            psListAdd(psfList, PS_LIST_TAIL, psf);
-            if (numCols == 0 && numRows == 0) {
-                numCols = ro->image->numCols;
-                numRows = ro->image->numRows;
-            }
-        }
-    }
-
-    bool convolve = false;              // Convolve the input images?
-    if (psfList->n > 0) {
-        psArray *psfArray = psListToArray(psfList); // Array of PSFs
-        outPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel,
-                               psfOrder, psfOrder);
-        psFree(psfArray);
-        if (!outPSF) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
-            // XXX Free stuff
-            return false;
-        }
-        convolve = true;
-        PS_ASSERT_PTR_NON_NULL(sources, false);
-    }
-
-    // Iterate through again to get the convolved images (or not)
-
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
-    psRandom *rng = psRandomAlloc(0, PS_RANDOM_TAUS); // Random number generator
+//    float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    pmCell *outCell = outRO->parent;    // Output cell
+    pmChip *outChip = outCell->parent;  // Output chip
+    pmFPA *outFPA = outChip->parent;    // Output FPA
+
     float totExposure = 0.0;            // Total exposure time
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
     psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
-
-    psMetadataIteratorSet(fileIter, PS_LIST_HEAD);
-    fileNum = 0;
-    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
-        assert(fileItem->type == PS_DATA_UNKNOWN);
-        pmFPAfile *inputFile = fileItem->data.V; // An input file
-        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        assert(ro);
+        pmFPA *fpa = ro->parent->parent->parent; // Parent FPA
 
         bool mdok;                      // Status of MD lookup
-        float weighting = psMetadataLookupF32(&mdok, inputFile->fpa->analysis,
-                                              "PPSTACK.WEIGHTING"); // Relative weighting
+        float weighting = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weight
         if (!mdok || !isfinite(weighting)) {
-            psWarning("No WEIGHTING supplied for image %d --- set to unity.", fileNum);
+            psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
             weighting = 1.0;
         }
-        float scale = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.SCALE"); // Rel. scale
+
+#if 0
+        // I don't think 'scale' is needed, since the convolution matches the scales
+        float scale = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.SCALE"); // Rel. scale
         if (!mdok || !isfinite(scale)) {
-            psWarning("No SCALE supplied for image %d --- set to unity.", fileNum);
+            psWarning("No SCALE supplied for image %d --- set to unity.", i);
             scale = 1.0;
         }
+#endif
 
         float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
-#if 0
         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(fileIter);
             psFree(fpaList);
             psFree(cellList);
-            psFree(outRO);
             psFree(stack);
             return false;
         }
-#endif
         totExposure += exposure;        // Total exposure time
 
-        // Generate convolved version of input
-        pmReadout *convolved;
-        if (convolve) {
-            convolved = pmReadoutAlloc(NULL); // Convolved version of input readout
-            // Background subtraction, scaling and normalisation is performed automatically by the image
-            // matching
-            if (!ppStackMatch(convolved, ro, sources, outPSF, config)) {
-                psWarning("Unable to match image %d --- ignoring.", fileNum);
-                psErrorClear();
-                psFree(convolved);
-                // XXX Free the bad image so it's not taking up good memory!
-                continue;
-            }
-
-#ifdef CONVOLUTION_FILES
-            if (convolved->image) {
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "convolved%03d_image.fits", fileNum);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, convolved->image, 0, NULL);
-                psFitsClose(fits);
-            }
-            if (convolved->mask) {
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "convolved%03d_mask.fits", fileNum);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, convolved->mask, 0, NULL);
-                psFitsClose(fits);
-            }
-            if (convolved->weight) {
-                psString name = NULL;           // Name of image
-                psStringAppend(&name, "convolved%03d_weight.fits", fileNum);
-                psFits *fits = psFitsOpen(name, "w");
-                psFree(name);
-                psFitsWriteImage(fits, NULL, convolved->weight, 0, NULL);
-                psFitsClose(fits);
-            }
-#endif // CONVOLUTION_FILES
-
-        } else {
-            // No convolution --- just use the unconvolved images as the "convolved" images, with some tweaks
-            convolved = psMemIncrRefCounter(ro);
-            sources = NULL;
-
-            // Brain-dead background subtraction
-            if (!psImageBackground(stats, NULL, 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(fileIter);
-                psFree(fpaList);
-                psFree(cellList);
-                psFree(stack);
-                psFree(outRO);
-                psFree(convolved);
-                return false;
-            }
-            psTrace("ppStack", 3, "Background for image %d is %f\n", fileNum, stats->robustMedian);
-            (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32));
-
-            // Apply scaling
-            (void)psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(1.0 / scale, 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(fileIter);
-                psFree(fpaList);
-                psFree(cellList);
-                psFree(outRO);
-                psFree(convolved);
-                psFree(stack);
-                return false;
-            }
-
-            (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32));
-            if (ro->weight) {
-                (void)psBinaryOp(ro->weight, ro->weight, "/", psScalarAlloc(exposure, PS_TYPE_F32));
-            }
-        }
-
-        if (fileNum == 0) {
+        if (i == 0) {
             // Copy astrometry over
-            pmFPA *fpa = ro->parent->parent->parent; // Template FPA
             pmHDU *hdu = fpa->hdu; // Template HDU
             pmHDU *outHDU = outFPA->hdu; // Output HDU
@@ -224 +77 @@
                 psWarning("Unable to find HDU at FPA level to copy astrometry.");
             } else {
-                if (!pmAstromReadWCS(outFPA, outCell->parent, hdu->header, 1.0)) {
+                if (!pmAstromReadWCS(outFPA, outChip, hdu->header, 1.0)) {
                     psErrorClear();
                     psWarning("Unable to read WCS astrometry from input FPA.");
@@ -231 +84 @@
                         outHDU->header = psMetadataAlloc();
                     }
-                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
+                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
                         psErrorClear();
                         psWarning("Unable to write WCS astrometry to output FPA.");
@@ -239 +92 @@
         }
 
-        // Don't need the original any more!
-        psFree(ro);
-
         // Ensure there is a mask, or pmStackCombine will complain
-        if (!convolved->mask) {
-            convolved->mask = psImageAlloc(convolved->image->numCols, convolved->image->numRows,
-                                           PS_TYPE_MASK);
-            psImageInit(convolved->mask, 0);
-        }
-
-        psListAdd(fpaList, PS_LIST_TAIL, inputFile->fpa);
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        psListAdd(fpaList, PS_LIST_TAIL, fpa);
         psListAdd(cellList, PS_LIST_TAIL, ro->parent);
 
-        pmStackData *data = pmStackDataAlloc(convolved, weighting); // Data to stack
-        psFree(convolved);
-        psArrayAdd(stack, ARRAY_BUFFER, data);
-        psFree(data);                   // Drop reference
-
-        fileNum++;
-    }
-    psFree(fileIter);
-    psFree(stats);
-    psFree(rng);
-
-    if (fileNum == 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not enough good files to combine.");
-        psFree(fpaList);
-        psFree(cellList);
-        psFree(stack);
-        psFree(outRO);
-        return false;
+        stack->data[i] = pmStackDataAlloc(ro, weighting);
     }
 
@@ -277 +109 @@
         psFree(cellList);
         psFree(stack);
-        psFree(outRO);
         return false;
     }
@@ -304 +135 @@
 #endif
 
-    for (int i = 0; i < stack->n; i++) {
+
+
+    // XXX THIS NEEDS WORK
+#if 0
+    for (int i = 0; i < num; i++) {
         pmStackData *data = stack->data[i]; // Data for this image
         pmReadout *readout = data->readout; // Readout for this image
+
+
+        // XXX Need some other mechanism to get the subtraction kernel
+
 
         // Extract the regions and solutions used in the image matching
@@ -344 +183 @@
         psFree(regions);
     }
+#endif
+
 
 #ifdef REJECTION_FILES
@@ -370 +211 @@
         psFree(cellList);
         psFree(stack);
-        psFree(outRO);
         return false;
-    }
-
-    if (!convolve) {
-        // Restore image to counts using the total exposure time
-        (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
-        if (outRO->weight) {
-            (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
-        }
     }
 
@@ -397 +229 @@
     psFree(cellList);
 
-    if (photometry) {
-        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
-
-        if (!psphotReadout(config, view)) {
-            psWarning("Unable to perform photometry on stacked image.");
-            psErrorStackPrint(stderr, "Error stack from photometry:");
-            psErrorClear();
-        }
-
-        pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
-    }
-
     psFree(stack);
-    psFree(outRO);                      // Drop reference
 
     return true;