Changeset 17227

Timestamp:

Mar 28, 2008, 5:08:31 PM (18 years ago)

Author:

Paul Price

Message:

Merging in pap_branch_080320 --- modernised version of ppMerge.

Location:

trunk/ppMerge/src

Files:

• Makefile.am (modified) (1 diff)
• ppMerge.c (modified) (3 diffs)
• ppMerge.h (modified) (1 diff)
• ppMergeArguments.c (added)
• ppMergeCamera.c (added)
• ppMergeFiles.c (added)
• ppMergeLoop.c (added)
• ppMergeMask.c (modified) (2 diffs)
• ppMergeScaleZero.c (modified) (1 diff)

trunk/ppMerge/src/Makefile.am

@@ -6 +6 @@
 ppMerge_SOURCES =               \
         ppMerge.c               \
-        ppMergeCheckInputs.c    \
-        ppMergeCombine.c        \
-        ppMergeConfig.c         \
-        ppMergeData.c           \
-        ppMergeMask.c           \
-        ppMergeMaskSuspect.c            \
-        ppMergeMaskBad.c                \
-        ppMergeMaskWrite.c              \
-        ppMergeMaskGrow.c               \
-        ppMergeMaskAverageConcepts.c    \
-        ppMergeOptions.c        \
+        ppMergeArguments.c      \
+        ppMergeCamera.c         \
+        ppMergeFiles.c          \
         ppMergeScaleZero.c      \
-        ppMergeVersion.c
+        ppMergeLoop.c           \
+        ppMergeMask.c
 
 
 noinst_HEADERS =                \
-        ppMerge.h               \
-        ppMergeCheckInputs.h    \
-        ppMergeCombine.h        \
-        ppMergeConfig.h         \
-        ppMergeData.h           \
-        ppMergeMask.h           \
-        ppMergeOptions.h        \
-        ppMergeScaleZero.h      \
-        ppMergeVersion.h
+        ppMerge.h
 
 CLEANFILES = *~

trunk/ppMerge/src/ppMerge.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <string.h>
 #include <pslib.h>
 #include <psmodules.h>
 
 #include "ppMerge.h"
-#include "ppMergeConfig.h"
-#include "ppMergeData.h"
-#include "ppMergeOptions.h"
-#include "ppMergeCheckInputs.h"
-#include "ppMergeCombine.h"
-#include "ppMergeScaleZero.h"
-#include "ppMergeMask.h"
 
 //#include "ppMem.h"
@@ -27 +21 @@
 {
     psLibInit(NULL);
-    psMemSetThreadSafety(false);           // Turn off thread safety, for more
+    psMemSetThreadSafety(false);
     psTimerStart(TIMERNAME);
 
-    // Parse the configuration and arguments
-    // Open the input image(s)
-    // Determine camera, format from header if not already defined
-    // Construct camera in preparation for reading
-    pmConfig *config = ppMergeConfig(argc, argv);
+    psExit exitValue = PS_EXIT_SUCCESS; // Exit value for program
+
+    pmConfig *config = pmConfigRead(&argc, argv, PPMERGE_RECIPE); // Configuration
     if (!config) {
-        psErrorStackPrint(stderr, "Unable to get configuration.");
-        exit(EXIT_FAILURE);
+        psErrorStackPrint(stderr, "Error reading configuration.");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
     }
 
+    if (!ppMergeArguments(argc, argv, config)) {
+        psErrorStackPrint(stderr, "Error reading arguments.");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame
+    switch (type) {
+      case PPMERGE_TYPE_MASK:
+        if (!ppMergeMask(config)) {
+            psErrorStackPrint(stderr, "Error generating mask.");
+            exitValue = PS_EXIT_DATA_ERROR;
+            goto die;
+        }
+        break;
+      case PPMERGE_TYPE_BIAS:
+      case PPMERGE_TYPE_DARK:
+      case PPMERGE_TYPE_SHUTTER:
+      case PPMERGE_TYPE_FLAT:
+      case PPMERGE_TYPE_FRINGE:
+        if (!ppMergeScaleZero(config)) {
+            psErrorStackPrint(stderr, "Error getting scale and zero-points.");
+            exitValue = PS_EXIT_DATA_ERROR;
+            goto die;
+        }
+        if (!ppMergeLoop(config)) {
+            psErrorStackPrint(stderr, "Error performing merge.");
+            exitValue = PS_EXIT_PROG_ERROR;
+            goto die;
+        }
+        break;
+      default:
+        psAbort("Invalid frame type: %x", type);
+    }
+
+
+    // Output the statistics
+    bool mdok;                          // Status of MD lookup
+    psString statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS.NAME"); // Statistics file name
+    if (mdok && statsName && strlen(statsName) > 0) {
+        psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename
+        FILE *statsFile = fopen(resolved, "w"); // Output statistics file
+        if (!statsFile) {
+            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.", resolved);
+            psFree(resolved);
+            exitValue = PS_EXIT_CONFIG_ERROR;
+            goto die;
+        }
+        psFree(resolved);
+        psMetadata *stats = psMetadataLookupMetadata(&mdok, config->arguments, "STATS.DATA"); // Statistics
+        if (!stats) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find statistics");
+            exitValue = PS_EXIT_PROG_ERROR;
+            goto die;
+        }
+        psString statsOut = psMetadataConfigFormat(stats); // String to write out
+        fprintf(statsFile, "%s", statsOut);
+        psFree(statsOut);
+        fclose(statsFile);
+    }
+
+
+
+
+#if 0
     // Set various tasks (define optional operations)
     ppMergeOptions *options = ppMergeOptionsParse(config);
@@ -81 +139 @@
     }
 
-#if 0
-    pmFPAPrint(stdout, data->out, true, true);
-#endif
-
     // Cleaning up
     psFree(data);
     psFree(options);
+#endif
+
+ die:
+    psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub"));
+    psTimerStop();
+
     psFree(config);
-
-    pmConceptsDone();
     pmConfigDone();
     psLibFinalize();
 
-//    ppMemCheck();
-
-    return EXIT_SUCCESS;
+    exit(exitValue);
 }

trunk/ppMerge/src/ppMerge.h

@@ -2 +2 @@
 #define PP_MERGE_H
 
-#define TIMERNAME "ppMerge"
-#define PPMERGE_RECIPE "PPMERGE"
+#define TIMERNAME "ppMerge"             // Name for timer
+#define PPMERGE_RECIPE "PPMERGE"        // Recipe name
+
+// Type of frame to merge
+typedef enum {
+    PPMERGE_TYPE_UNKNOWN,               // Unknown type
+    PPMERGE_TYPE_BIAS,                  // Bias frame
+    PPMERGE_TYPE_DARK,                  // (Multi-)Dark frame
+    PPMERGE_TYPE_MASK,                  // Mask frame
+    PPMERGE_TYPE_SHUTTER,               // Shutter frame
+    PPMERGE_TYPE_FLAT,                  // Flat-field frame (dome or sky)
+    PPMERGE_TYPE_FRINGE,                // Fringe frame
+} ppMergeType;
+
+// Files, for activation
+typedef enum {
+    PPMERGE_FILES_ALL,                  // All files
+    PPMERGE_FILES_INPUT,                // Input files
+    PPMERGE_FILES_OUTPUT                // Output files
+} ppMergeFiles;
+
+// Parse command-line arguments and recipe
+bool ppMergeArguments(int argc, char *argv[], // Command-line arguments
+                      pmConfig *config  // Configuration
+    );
+
+// Set up camera files
+bool ppMergeCamera(pmConfig *config     // Configuration
+    );
+
+// Measure scale and zero-points
+bool ppMergeScaleZero(pmConfig *config  // Configuration
+    );
+
+// Main loop to do the merging
+bool ppMergeLoop(pmConfig *config       // Configuration
+    );
+
+// Main loop for masks
+bool ppMergeMask(pmConfig *config       // Configuration
+    );
+
+// Read nominated input file
+bool ppMergeFileReadInput(const pmConfig *config, // Configuration
+                          pmReadout *readout, // Readout into which to read
+                          int num,      // Number of file in sequence
+                          int rows      // Number of rows to read at once
+    );
+
+// Open nominated input file
+bool ppMergeFileOpenInput(pmConfig *config, // Configuration
+                          const pmFPAview *view, // View to open
+                          int num       // Number of file in sequence
+    );
+
+// Set the data level for files specified by name; return an array of the files
+psArray *ppMergeFileDataLevel(const pmConfig *config, // Configuration
+                              const char *name, // Name of files
+                              pmFPALevel level // Level for file data level
+    );
+
+// Activate/deactivate a list of files
+bool ppMergeFileActivate(const pmConfig *config, // Configuration
+                         ppMergeFiles files, // Files to turn on/off
+                         bool state     // Activation state
+    );
+
+// Activate/deactivate a single element for a list; return array of files
+psArray *ppMergeFileActivateSingle(const pmConfig *config, // Configuration
+                                   ppMergeFiles files, // Files to turn on/off
+                                   bool state,   // Activation state
+                                   int num // Number of file in sequence
+    );
+
+// Return name of output pmFPAfile
+psString ppMergeOutputFile(const pmConfig *config // Configuration
+    );
 
 #endif

trunk/ppMerge/src/ppMergeMask.c

@@ -5 +5 @@
 #include <stdio.h>
 #include <string.h>
-#include <unistd.h>
 #include <assert.h>
+
 #include <pslib.h>
 #include <psmodules.h>
@@ -12 +12 @@
 
 #include "ppMerge.h"
-#include "ppMergeData.h"
-#include "ppMergeMask.h"
-
-
-// Generate a mask
-bool ppMergeMask(ppMergeData *data,  // Data
-                 ppMergeOptions *options, // Options
-                 pmConfig *config    // Configuration
-    )
+
+static bool mergeMask(pmConfig *config, // Configuration
+                      const pmFPAview *view, // View to chip
+                      bool writeOut,     // Write output?
+                      psRandom *rng,    // Random number generator
+                      psMetadata *stats // Statistics output
+                      )
 {
-    for (int i = 0; i < 2; i++) {
-        if (!ppMergeMaskSuspect (data, options, config)) {
-            psError(PS_ERR_UNKNOWN, true, "failed on mask suspect.\n");
-            return false;
-        }
-
-        if (!ppMergeMaskBad (data, options, config)) {
-            psError(PS_ERR_UNKNOWN, true, "failed on mask bad.\n");
-            return false;
-        }
-    }
-
-    if (!ppMergeMaskAverageConcepts (data, options, config)) {
-        psError(PS_ERR_UNKNOWN, true, "failed on average concepts.\n");
-        return false;
-    }
-
-    if (!ppMergeMaskGrow (data, options, config)) {
-        psError(PS_ERR_UNKNOWN, true, "failed on mask grow.\n");
-        return false;
-    }
-
-    if (!ppMergeMaskWrite (data, options, config)) {
-        psError(PS_ERR_UNKNOWN, true, "failed on mask write.\n");
-        return false;
-    }
+    assert(config);
+    assert(view);
+    assert(view->chip != -1 && view->cell == -1 && view->readout == -1);
+
+    bool mdok;                          // Status of MD lookup
+    int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of input files
+    psStatsOptions meanStat = psMetadataLookupS32(NULL, config->arguments, "MEAN"); // Statistic for mean
+    psStatsOptions stdevStat = psMetadataLookupS32(NULL, config->arguments, "STDEV"); // Statistic for stdev
+    psMaskType maskVal = psMetadataLookupU8(NULL, config->arguments, "MASKVAL"); // Value to mask
+    int sample = psMetadataLookupS32(NULL, config->arguments, "SAMPLE"); // Size of sample for statistics
+    bool chipStats = psMetadataLookupBool(&mdok, config->arguments, "MASK.CHIPSTATS"); // Statistics on chip?
+    float maskSuspect = psMetadataLookupF32(NULL, config->arguments, "MASK.SUSPECT"); // Threshold for suspect pixels
+    float maskBad = psMetadataLookupF32(NULL, config->arguments, "MASK.BAD"); // Threshold for bad pixels
+    pmMaskIdentifyMode maskMode = psMetadataLookupS32(NULL, config->arguments, "MASK.MODE"); // Mode for identifying bad pixels
+    int maskGrow = psMetadataLookupS32(NULL, config->arguments, "MASK.GROW"); // Radius to grow mask
+    psMaskType maskGrowVal = psMetadataLookupU8(NULL, config->arguments, "MASK.GROWVAL"); // Value for grown mask
+
+    psStats *statistics = psStatsAlloc(meanStat | stdevStat); // Statistics for background
+
+    psString outName = ppMergeOutputFile(config); // Name of output file
+    pmChip *outChip = pmFPAfileThisChip(config->files, view, outName); // Output chip
+    psFree(outName);
+
+    // For each input file, get the statistics, which can be calculated at the chip or cell levels
+    psVector *values = psVectorAlloc(sample, PS_TYPE_F32); // Pixel values for statistics
+    pmFPAview *inView = pmFPAviewAlloc(0); // View for input
+    for (int i = 0; i < numFiles; i++) {
+        pmFPAfileActivate(config->files, false, NULL);
+        psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Input files
+        pmFPAfile *input = files->data[0]; // Input file
+        psFree(files);
+        pmFPA *inFPA = input->fpa;  // Input FPA
+        *inView = *view;
+
+        int valueIndex = 0;             // Index for vector of pixel values
+
+        pmCell *inCell;                 // Input cell
+        while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
+            pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell
+            if (!hdu || hdu->blankPHU) {
+                // No data here
+                continue;
+            }
+            psTrace("ppMerge", 1, "Getting suspect pixels for file %d chip %d cell %d",
+                    i, inView->chip, inView->cell);
+
+            if (!pmFPAfileIOChecks(config, inView, PM_FPA_BEFORE)) {
+                psFree(inView);
+                goto MERGE_MASK_ERROR;
+            }
+
+            if (!ppMergeFileOpenInput(config, view, i)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i);
+                psFree(inView);
+                goto MERGE_MASK_ERROR;
+            }
+
+            if (inCell->readouts->n > 1) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "File %d chip %d cell %d contains more than one readout (%ld)",
+                        i, inView->chip, inView->cell, inCell->readouts->n);
+                psFree(inView);
+                goto MERGE_MASK_ERROR;
+            }
+
+            pmReadout *readout;
+            if (inCell->readouts && inCell->readouts->n == 1) {
+                readout = psMemIncrRefCounter(inCell->readouts->data[0]); // Input readout
+            } else {
+                readout = pmReadoutAlloc(inCell);
+            }
+
+            if (!ppMergeFileReadInput(config, readout, i, 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to read readout %d", i);
+                psFree(inView);
+                psFree(readout);
+                goto MERGE_MASK_ERROR;
+            }
+
+            pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell
+            pmReadout *outRO = NULL;    // Output readout
+            if (outCell->readouts && outCell->readouts->n == 1) {
+                outRO = psMemIncrRefCounter(outCell->readouts->data[0]);
+            } else {
+                outRO = pmReadoutAlloc(outCell);
+            }
+            psImage *outMask = outRO->mask;    // Output mask image (for iterative generation of mask)
+
+            psImage *image = readout->image, *mask = readout->mask; // Image and mask
+            int numCols = readout->image->numCols, numRows = readout->image->numRows; // Image size
+            int numPix = numCols * numRows; // Number of pixels
+            int numCells = chipStats ? readout->parent->parent->cells->n : 1; // Number of cells
+            int num = PS_MIN(numPix, sample / numCells); // Number of values to add
+            if (!chipStats) {
+                valueIndex = 0;
+            }
+            for (int i = 0; i < num; i++) {
+                int pixel = numPix * psRandomUniform(rng);
+                int x = pixel % numCols;
+                int y = pixel / numCols;
+                if ((mask && (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal)) ||
+                    !isfinite(image->data.F32[y][x]) ||
+                    (outMask && (outMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) {
+                    continue;
+                }
+
+                values->data.F32[valueIndex++] = image->data.F32[y][x];
+            }
+
+            if (!chipStats) {
+                values->n = valueIndex;
+                if (!psVectorStats(statistics, values, NULL, NULL, 0)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on readout.");
+                    psFree(inView);
+                    psFree(readout);
+                    goto MERGE_MASK_ERROR;
+                }
+
+                if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat),
+                                             psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i);
+                    psFree(inView);
+                    psFree(readout);
+                    goto MERGE_MASK_ERROR;
+                }
+                pmCellFreeData(inCell);
+
+                if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) {
+                    psFree(inView);
+                    psFree(readout);
+                    goto MERGE_MASK_ERROR;
+                }
+            }
+            psFree(readout);
+            psFree(outRO);
+        }
+
+        // Additional run through cells if we want chip-level statistics
+        if (chipStats && valueIndex > 0) {
+            values->n = valueIndex;
+            if (!psVectorStats(statistics, values, NULL, NULL, 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on chip.");
+                goto MERGE_MASK_ERROR;
+            }
+            inView->cell = -1;
+            while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
+                pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell
+                if (!hdu || hdu->blankPHU) {
+                    // No data here
+                    continue;
+                }
+                pmReadout *readout = inCell->readouts->data[0]; // Readout of interest
+
+                inView->readout = 0;
+                pmReadout *outRO = pmFPAfileThisReadout(config->files, inView, "PPMERGE.OUTPUT.MASK");
+
+                if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat),
+                                             psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i);
+                    goto MERGE_MASK_ERROR;
+                }
+                pmCellFreeData(inCell);
+
+                inView->readout = -1;
+                if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) {
+                    psFree(inView);
+                    goto MERGE_MASK_ERROR;
+                }
+            }
+        }
+    }
+    psFree(inView);
+    psFree(statistics); statistics = NULL;
+    psFree(values); values = NULL;
+
+
+    // Another run through the chip to threshold on the suspects
+    pmFPAfileActivate(config->files, false, NULL);
+    if (writeOut) {
+        ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true);
+    }
+    pmFPA *outFPA = outChip->parent;    // Output FPA
+    pmCell *outCell;                    // Output cell
+    pmFPAview *outView = pmFPAviewAlloc(0); // View into output FPA
+    *outView = *view;
+    while ((outCell = pmFPAviewNextCell(outView, outFPA, 1))) {
+        pmHDU *hdu = pmHDUFromCell(outCell); // HDU for cell
+        if (!hdu || hdu->blankPHU) {
+            // No data here
+            continue;
+        }
+
+        psTrace("ppMerge", 1, "Getting bad pixels for chip %d cell %d", outView->chip, outView->cell);
+
+        assert(outCell->readouts && outCell->readouts->n == 1);
+        pmReadout *outRO = outCell->readouts->data[0]; // Output readout
+        if (!pmMaskIdentifyBadPixels(outRO, maskVal, maskBad, maskMode)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to mask bad pixels");
+            goto MERGE_MASK_ERROR;
+        }
+
+        // Supplementary outputs
+        {
+            // The counts image is fairly useless, but it preserves the model
+            pmCell *countsCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.COUNT");
+            pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel
+            countsRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32);
+            psImageInit(countsRO->image, numFiles);
+            countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true;
+            psFree(countsRO);
+
+            pmCell *sigmaCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.SIGMA");
+            pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with suspect image
+            psImage *suspect = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT);
+            sigmaRO->image = psImageCopy(sigmaRO->image, suspect, PS_TYPE_F32);
+            psMetadataRemoveKey(outRO->analysis, PM_MASK_ANALYSIS_SUSPECT);
+            sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true;
+            psFree(sigmaRO);
+        }
+
+        if (maskGrowVal > 0) {
+            psImage *grown = psImageGrowMask(NULL, outRO->mask, maskVal, maskGrow, maskGrowVal); // Grown mask
+            psFree(outRO->mask);
+            outRO->mask = grown;
+        }
+
+        if (writeOut) {
+            if (!pmFPAfileIOChecks(config, outView, PM_FPA_BEFORE)) {
+                psFree(outView);
+                goto MERGE_MASK_ERROR;
+            }
+
+            outRO->data_exists = outCell->data_exists = outChip->data_exists = true;
+
+            // Average concepts
+            psList *cells = psListAlloc(NULL); // List of cells, for concept averaging
+            for (int i = 0; i < numFiles; i++) {
+                pmFPAfile *inFile = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file
+                pmCell *inCell = pmFPAviewThisCell(outView, inFile->fpa); // Input cell
+                psListAdd(cells, PS_LIST_TAIL, inCell);
+            }
+            if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, true)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
+                psFree(cells);
+                psFree(outRO);
+                psFree(outView);
+                return false;
+            }
+            psFree(cells);
+
+            // Statistics on the merged cell using a fake image
+            outRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32);
+            psImageInit(outRO->image, 1.0);
+            if (!ppStatsFPA(stats, outRO->parent->parent->parent, outView,
+                            maskVal | pmConfigMask("BLANK", config), config)) {
+                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.");
+                psFree(outRO);
+                psFree(outView);
+                return false;
+            }
+            psFree(outRO->image);
+            outRO->image = NULL;
+
+            // Write
+            if (!pmFPAfileIOChecks(config, outView, PM_FPA_AFTER)) {
+                psFree(outView);
+                goto MERGE_MASK_ERROR;
+            }
+        }
+    }
+    psFree(outView);
+
+    ppMergeFileActivate(config, PPMERGE_FILES_ALL, true);
 
     return true;
+
+
+MERGE_MASK_ERROR:
+    psFree(statistics);
+    psFree(values);
+    return false;
 }
 
-bool ppMergeMaskReadoutStats(const pmReadout *readout, 
-                             const pmReadout *output, 
-                             ppMergeOptions *options, // Options
-                             psRandom *rng)
+
+
+
+
+bool ppMergeMask(pmConfig *config)
 {
-    PS_ASSERT_PTR_NON_NULL(readout, false);
-    PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
-    PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false);
-    PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
-    if (readout->mask) {
-        PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false);
-        PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false);
-    }
-
-    psImage *mask = NULL;
-    psImage *image = readout->image;    // Image of interest
-
-    if (output) {
-        mask = output->mask;      // Corresponding mask
-    }
-
-    if (rng) {
-        psMemIncrRefCounter(rng);
-    } else {
-        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
-    }
-
-    // XXX note that this now will accept any of several stats options
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    stats->nSubsample = options->sample;
-
-    if (!psImageBackground(stats, NULL, image, mask, options->combine->maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Failure to measure image statistics.\n");
-        psFree(stats);
-        psFree(rng);
-        return false;
-    }
-    if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
-        psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
-        psFree(stats);
-        psFree(rng);
-        return false;
-    }
-
-    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
-    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
-
+    assert(config);
+
+    bool mdok;                          // Status of MD lookup
+    int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
+    bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights?
+    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Number of rejection iterations
+
+    PS_ASSERT_INT_POSITIVE(iter, false);
+
+    pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+
+    psMetadata *stats = NULL;           // Statistics for output
+    if (psMetadataLookup(config->arguments, "STATS.NAME")) {
+        stats = psMetadataAlloc();
+        psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "STATS.DATA", 0, "Statistics output", stats);
+    }
+
+    psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images
+    psFree(inputs);
+    if (haveMasks) {
+        psArray *masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT);
+        psFree(masks);
+    }
+    if (haveWeights) {
+        psArray *weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT);
+        psFree(weights);
+    }
+
+    if (!ppMergeFileActivate(config, PPMERGE_FILES_INPUT, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to activate files.");
+        goto PPMERGE_MASK_ERROR;
+    }
+    if (!ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to activate files.");
+        goto PPMERGE_MASK_ERROR;
+    }
+
+    psString outName = ppMergeOutputFile(config); // Name of output file
+    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, outName); // Output file
+    psFree(outName);
+    assert(output && output->fpa);
+    pmFPA *outFPA = output->fpa;        // Output FPA
+
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        goto PPMERGE_MASK_ERROR;
+    }
+    pmChip *outChip;                    // Chip of interest
+    while ((outChip = pmFPAviewNextChip(view, outFPA, 1))) {
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            goto PPMERGE_MASK_ERROR;
+        }
+
+        for (int i = 0; i < iter; i++) {
+            if (!mergeMask(config, view, (i == iter - 1), rng, stats)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to merge chip %d", view->chip);
+                goto PPMERGE_MASK_ERROR;
+            }
+        }
+
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            goto PPMERGE_MASK_ERROR;
+        }
+    }
+
+    psList *fpaList = psListAlloc(NULL);// List of FPAs for concept averaging
+    for (int i = 0; i < numFiles; i++) {
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file
+        psListAdd(fpaList, PS_LIST_TAIL, file->fpa);
+    }
+    if (!pmConceptsAverageFPAs(outFPA, fpaList)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts.");
+        psFree(fpaList);
+        goto PPMERGE_MASK_ERROR;
+    }
+    psFree(fpaList);
+
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        goto PPMERGE_MASK_ERROR;
+    }
+
+    psFree(stats);
+    psFree(view);
     psFree(rng);
+
     return true;
+
+PPMERGE_MASK_ERROR:
+    psFree(stats);
+    psFree(view);
+    psFree(rng);
+    return false;
 }
 
-bool ppMergeMaskChipStats (const pmChip *chip,
-                           const pmChip *output, 
-                           ppMergeOptions *options,
-                           psRandom *rng)
-{
-    PS_ASSERT_PTR_NON_NULL(chip, false);
-
-    // XXX note that this now will accept any of several stats options
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-
-    if (rng) {
-        psMemIncrRefCounter(rng);
-    } else {
-        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
-    }
-
-    // accumulate a vector of data values using options->sample per readout
-    psVector *values = psVectorAllocEmpty(options->sample, PS_TYPE_F32); // Vector containing subsample
-
-    for (int nCell = 0; nCell < chip->cells->n; nCell++) {
-        pmCell *cell = chip->cells->data[nCell];
-        if (cell->readouts->n == 0) continue;
-
-        for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
-            pmReadout *readout = cell->readouts->data[nReadout];
-            if (!readout->image) continue;
-
-            psTrace("ppMerge", 4, "Measure statistics for cell %d, readout %d\n", nCell, nReadout);
-
-            psImage *image = readout->image;    // Image of interest
-
-            pmCell *cellOutput = output->cells->data[nCell];
-            pmReadout *readoutOutput = NULL;
-            psImage *mask = NULL;
-            if (cellOutput->readouts->n > 0) {
-                readoutOutput = cellOutput->readouts->data[0];
-                mask = readoutOutput->mask;      // Corresponding mask
-            }
-
-            // Size of image
-            long nx = image->numCols;
-            long ny = image->numRows;
-            const int Npixels = nx*ny;  // Total number of pixels
-            const int Nsubset = PS_MIN(options->sample, Npixels); // Number of pixels in subset
-            
-            values = psVectorRealloc (values, values->n + Nsubset); // make sure we have enough space
-
-            // select a subset of the image pixels to measure the stats
-            for (long i = 0; i < Nsubset; i++) {
-                double frnd = psRandomUniform(rng);
-                int pixel = Npixels * frnd;
-                int ix = pixel % nx;
-                int iy = pixel / nx;
-
-                if (!isfinite(image->data.F32[iy][ix])) continue;
-                if (mask && (mask->data.U8[iy][ix] & options->combine->maskVal)) continue;
-
-                float value = image->data.F32[iy][ix];
-                values->data.F32[values->n] = value;
-                values->n ++;
-            }
-        }
-    }
-
-    // no valid data, skip the chip
-    if (!values->n) {
-        psFree(values);
-        psFree(stats);
-        psFree(rng);
-        return true;
-    }
-
-    // calculate the statistics
-    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for chip");
-        psFree(values);
-        psFree(stats);
-        psFree(rng);
-        return false;
-    }
-    if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
-        psError(PS_ERR_UNKNOWN, false, "invalid image statistics (nan).\n");
-        psFree(values);
-        psFree(stats);
-        psFree(rng);
-        return false;
-    }
-
-    // supply the stats to the readout analysis metadata
-    for (int nCell = 0; nCell < chip->cells->n; nCell++) {
-        pmCell *cell = chip->cells->data[nCell];
-        if (cell->readouts->n == 0) continue;
-
-        for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
-            pmReadout *readout = cell->readouts->data[nReadout];
-            if (!readout->image) continue;
-
-            psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
-            psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
-        }
-    }
-
-    psLogMsg ("ppMerge", PS_LOG_INFO, "statistics for chip: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
-
-    psFree(values);
-    psFree(stats);
-    psFree(rng);
-    return true;
-}

trunk/ppMerge/src/ppMergeScaleZero.c

@@ -10 +10 @@
 
 #include "ppMerge.h"
-#include "ppMergeScaleZero.h"
 
 // Get the scale and zero for each chip of each input
-bool ppMergeScaleZero(psImage **scales, // The scales for each integration/cell
-                      psImage **zeros, // The zeroes for each integration/cell
-                      psArray **shutters, // The shutter correction data for each cell
-                      ppMergeData *data,// The data
-                      const ppMergeOptions *options, // The options
-                      const pmConfig *config // The configuration
-    )
+bool ppMergeScaleZero(pmConfig *config)
 {
-    assert(data);
-    assert(options);
     assert(config);
 
-    if (!options->scale && !options->zero && !options->shutter) {
-        return true;                    // We did everything we were asked for
+    ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame
+    int numInputs = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    int numCells = psMetadataLookupS32(NULL, config->arguments, "INPUTS.CELLS"); // Number of cells
+    psStatsOptions meanStat = psMetadataLookupS32(NULL, config->arguments, "MEAN"); // Statistic for mean
+    psStatsOptions stdevStat = psMetadataLookupS32(NULL, config->arguments, "STDEV"); // Statistic for stdev
+    psMaskType maskVal = psMetadataLookupU8(NULL, config->arguments, "MASKVAL"); // Value to mask
+    int shutterSize = psMetadataLookupS32(NULL, config->arguments, "SHUTTER.SIZE"); // Size of shutter region
+
+    psVector *gains = NULL;             // Gains for each cell
+    psArray *shutters = NULL;           // Shutter data for each cell
+    psStats *stats = NULL;              // Statistics for background
+    psImage *background = NULL;         // Background measurements per cell per file
+
+    switch (type) {
+      case PPMERGE_TYPE_BIAS:
+      case PPMERGE_TYPE_DARK:
+        // Nothing to measure
+        return true;
+      case PPMERGE_TYPE_FLAT:
+      case PPMERGE_TYPE_FRINGE:
+        gains = psVectorAlloc(numCells, PS_TYPE_F32);
+        background = psImageAlloc(numCells, numInputs, PS_TYPE_F32);
+        psImageInit(background, NAN);
+        stats = psStatsAlloc(meanStat);
+        break;
+      case PPMERGE_TYPE_SHUTTER:
+        shutters = psArrayAlloc(numCells);
+        break;
+      case PPMERGE_TYPE_MASK:
+      default:
+        break;
     }
-
-    assert(config->camera);             // Need the camera configuration
-    assert(config->arguments);          // Need the list of files
-
-    psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
-    assert(filenames);                  // It should be here --- it's put here in ppMergeConfig
-
-    // Sanity checks
-    assert(!options->scale || scales);
-    assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F32 &&
-                                   (*scales)->numCols == data->numCells &&
-                                   (*scales)->numRows == filenames->n));
-    assert(!options->zero || zeros);
-    assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F32 &&
-                                 (*zeros)->numCols == data->numCells &&
-                                 (*zeros)->numRows == filenames->n));
-    assert(!options->shutter || shutters);
-    assert(!shutters || !*shutters || (*shutters)->n == data->numCells);
-
-    // Allocate the outputs
-    if (options->scale) {
-        if (*scales) {
-            psFree(*scales);
-        }
-        *scales = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32);
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+    pmFPAview *view = NULL;             // View into FPA
+
+    for (int i = 0; i < numInputs; i++) {
+        pmFPAfileActivate(config->files, false, NULL);
+        psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Activated files
+        pmFPAfile *input = files->data[0]; // Representative file; should be the image (not mask or weight)
+        pmFPA *fpa = input->fpa;        // FPA of interest
+        view = pmFPAviewAlloc(0);       // View to component of interest
+        int cellNum = 0;                // Index for cell
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            goto ERROR;
+        }
+        pmChip *chip;                   // Chip of interest
+        while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
+            if (!chip->file_exists) {
+                continue;
+            }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                goto ERROR;
+            }
+
+            pmCell *cell;               // Cell of interest
+            while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
+                if (!cell->file_exists) {
+                    continue;
+                }
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    goto ERROR;
+                }
+
+                if (!cell->data_exists) {
+                    continue;
+                }
+
+                if (cell->readouts->n > 1) {
+                    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                            "File %d chip %d cell %d contains more than one readout (%ld)",
+                            i, view->chip, view->cell, cell->readouts->n);
+                    goto ERROR;
+                }
+                pmReadout *readout = cell->readouts->data[0]; // Readout of interest
+
+                switch (type) {
+                  case PPMERGE_TYPE_FLAT:
+                  case PPMERGE_TYPE_FRINGE: {
+                      // Extract the gain
+                      float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
+                      if (!isfinite(gain)) {
+                          psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                                  "CELL.GAIN for file %d chip %d cell %d is not set.",
+                                  i, view->chip, view->cell);
+                          goto ERROR;
+                      }
+                      gains->data.F32[cellNum] = gain;
+
+                      // Measure the background
+                      if (!psImageBackground(stats, NULL, readout->image, readout->mask, maskVal, rng)) {
+                          psError(PS_ERR_UNKNOWN, false,
+                                  "Unable to get statistics for file %d chip %d cell %d",
+                                  i, view->chip, view->cell);
+                          goto ERROR;
+                      }
+                      background->data.F32[i][cellNum] = psStatsGetValue(stats, meanStat);
+                      break;
+                  }
+                  case PPMERGE_TYPE_SHUTTER: {
+                      pmShutterCorrectionData *shutter = shutters->data[cellNum]; // Shutter correction data
+                      if (!shutter) {
+                          shutter = pmShutterCorrectionDataAlloc(readout->image->numCols,
+                                                                 readout->image->numRows,
+                                                                 shutterSize);
+                          shutters->data[cellNum] = shutter;
+                      }
+                      if (!pmShutterCorrectionAddReadout(shutter, readout, meanStat, stdevStat,
+                                                         maskVal, rng)) {
+                          psError(PS_ERR_UNKNOWN, false,
+                                  "Can't add file %d chip %d cell %d to shutter correction.",
+                                  i, view->chip, view->cell);
+                          goto ERROR;
+                      }
+                      break;
+                  }
+                  case PPMERGE_TYPE_MASK:
+                  default:
+                    psAbort("Should never get here.");
+                }
+
+                cellNum++;
+
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                    goto ERROR;
+                }
+            }
+
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                goto ERROR;
+            }
+        }
+
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            goto ERROR;
+        }
+
+        psFree(view);
+
+#if 0
+        // Reset files for reading again
+        for (int i = 0; i < files->n; i++) {
+            pmFPAfile *file = files->data[i]; // File of interest
+        }
+#endif
+        psFree(files);
     }
-    if (options->zero) {
-        if (*zeros) {
-            psFree(*zeros);
-        }
-        *zeros = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32);
+
+    psFree(rng); rng = NULL;
+    psFree(stats); stats = NULL;
+
+    // Store results
+    switch (type) {
+      case PPMERGE_TYPE_FRINGE:
+        psMetadataAddImage(config->arguments, PS_LIST_TAIL, "ZEROS", 0,
+                           "Zero to subtract from each input cell", background);
+        // Flow through
+      case PPMERGE_TYPE_FLAT: {
+          // Need to normalize over the focal plane
+          if (psTraceGetLevel("ppMerge") > 9) {
+              for (int i = 0; i < gains->n; i++) {
+                  psTrace("ppMerge", 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]);
+              }
+          }
+          psVector *fluxes = NULL;        // Solution to fluxes
+          if (!pmFlatNormalize(&fluxes, &gains, background)) {
+              psError(PS_ERR_UNKNOWN, false, "Normalisation failed to converge --- continuing anyway.");
+              psFree(fluxes);
+              goto ERROR;
+          }
+
+          psMetadataAddVector(config->arguments, PS_LIST_TAIL, "SCALES", 0,
+                              "Scale to divide into each input file", fluxes);
+          psFree(fluxes);               // Drop reference
+          break;
+      }
+      case PPMERGE_TYPE_SHUTTER:
+        psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SHUTTER", 0,
+                           "Shutter data", shutters);
+        break;
+      case PPMERGE_TYPE_MASK:
+      default:
+        psAbort("Should never get here.");
     }
-    psRandom *rng = NULL;               // Random number generator
-    if (options->shutter) {
-        if (*shutters) {
-            psFree(*shutters);
-        }
-        *shutters = psArrayAlloc(data->numCells);
-        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
-    }
-
-    bool fromConcepts = false;          // Do we get the scale and zero points from the concepts?
-    bool first = true;                  // Are we on the first cell (that sets the standard for the rest)?
-    bool done = false;                  // Are we done going through the list?
-    bool mdok = true;                   // Status of MD lookup
-    for (long i = 0; i < data->in->n && !done; i++) {
-        pmFPA *fpa = data->in->data[i]; // The FPA
-        if (!fpa) {
-            continue;
-        }
-        long cellNum = -1;              // Number of the cell
-        psArray *chips = fpa->chips;    // The array of chips
-        for (long j = 0; j < chips->n && !done; j++) {
-            pmChip *chip = chips->data[j]; // The chip
-            if (!chip) {
-                continue;
-            }
-            psArray *cells = chip->cells; // The array of cells
-            for (long k = 0; k < cells->n && !done; k++) {
-                pmCell *cell = cells->data[k]; // The cell
-                if (!cell) {
-                    continue;
-                }
-                cellNum++;
-
-                if (options->scale) {
-                    float scale = psMetadataLookupF32(&mdok, cell->concepts, "PPMERGE.SCALE"); // The scale
-                    if (mdok && !isnan(scale)) {
-                        if (!first && !fromConcepts) {
-                            psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
-                                     "for some, but not all cells --- we will re-measure it for all cells.");
-                            done = true;
-                            continue;
-                        }
-                        fromConcepts = true;
-                        (*scales)->data.F32[i][cellNum] = scale;
-                        psTrace("ppMerge", 9, "Scale for input %ld, chip %ld, cell %ld: %f\n",
-                                i, j, k, scale);
-                    } else if (!first && fromConcepts) {
-                        psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
-                                 "for some, but not all cells --- we will re-measure it for all cells.");
-                        fromConcepts = false;
-                        done = true;
-                        continue;
-                    }
-                }
-
-                if (options->zero) {
-                    float zero = psMetadataLookupF32(&mdok, cell->concepts, "PPMERGE.ZERO"); // The zero
-                    if (mdok && !isnan(zero)) {
-                        if (!first && !fromConcepts) {
-                            psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
-                                     "for some, but not all cells --- we will re-measure it for all cells.");
-                            done = true;
-                            continue;
-                        }
-                        fromConcepts = true;
-                        (*zeros)->data.F32[i][cellNum] = zero;
-                        psTrace("ppMerge", 9, "Zero for input %ld, chip %ld, cell %ld: %f\n", i, j, k, zero);
-                    } else if (!first && fromConcepts) {
-                        psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set "
-                                 "for some, but not all cells --- we will re-measure it for all cells.");
-                        fromConcepts = false;
-                        done = true;
-                        continue;
-                    }
-                }
-
-                first = false;
-            }
-        }
-    }
-
-    if (fromConcepts) {
-        // We've already done everything we need to
-        psFree(rng);
-        return true;
-    }
-
-    psImage *background = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); // Background measurements
-    psImageInit(background, NAN);
-    psStats *bgStats = psStatsAlloc(options->mean); // Statistic to measure the background
-    psVector *gains = NULL;             // The gains for each cell
-    if (options->scale) {
-        gains = psVectorAlloc(data->numCells, PS_TYPE_F32);
-    }
-
-    pmFPAview *view = pmFPAviewAlloc(0);
-
-    bool status = true;                 // Status of getting the scale and zero --- did everything go right?
-    for (int i = 0; i < filenames->n; i++) {
-        psString name = filenames->data[i]; // The name of the file
-        if (!name || strlen(name) == 0) {
-            continue;
-        }
-        psTrace("ppMerge", 9, "Opening %s to get background...\n", name);
-        psFits *inFile = data->files->data[i]; // The FITS file to read
-        pmFPA *fpa = data->in->data[i]; // The FPA for this input
-        int cellNum = -1;               // Number of the cell
-        psArray *chips = fpa->chips;    // Array of chips
-        for (int j = 0; j < chips->n; j++) {
-            pmChip *chip = chips->data[j]; // The chip of interest
-            if (!chip) {
-                continue;
-            }
-            psArray *cells = chip->cells; // Array of cells
-            for (int k = 0; k < cells->n; k++) {
-                pmCell *cell = cells->data[k]; // The cell of interest
-                if (!cell) {
-                    continue;
-                }
-                cellNum++;
-
-                if (!pmCellReadHeader(cell, inFile)) {
-                    continue;
-                }
-
-                // Scaling by the background
-                if (options->scale || options->zero) {
-                    if (!pmCellRead(cell, inFile, config->database)) {
-                        // Nothing here
-                        pmCellFreeData(cell);
-                        continue;
-                    }
-
-                    if (cell->readouts->n > 1) {
-                        psWarning("File %s chip %d cell %d contains more than one "
-                                 "readout --- ignoring all but the first.\n", name, j, k);
-                        status = false;
-                    }
-
-                    pmReadout *readout = cell->readouts->data[0]; // The readout of interest
-                    psImage *image = readout->image; // The pixels of interest
-                    if (!image) {
-                        pmCellFreeData(cell);
-                        continue;
-                    }
-
-                    // Get the gain
-                    if (options->scale) {
-                        bool mdok = true;   // Status of MD lookup
-                        gains->data.F32[cellNum] = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN");
-                        if (!mdok || isnan(gains->data.F32[cellNum])) {
-                            psWarning("CELL.GAIN for file %s chip %d cell %d is not "
-                                     "set.\n", name, j, k);
-                            gains->data.F32[cellNum] = NAN;
-                            status = false;
-                        }
-                    }
-
-                    // Get the background
-                    int sampleSize = (image->numCols * image->numRows) / options->sample; // Size of sample
-                    psVector *sample = psVectorAlloc(sampleSize, PS_TYPE_F32); // Sample of the image
-                    psVector *sampleMask = NULL; // Mask for sample
-                    if (readout->mask) {
-                        sampleMask = psVectorAlloc(sampleSize, PS_TYPE_U8);
-                    }
-                    psImage *mask = readout->mask; // The mask image
-                    for (long i = 0; i < sampleSize; i++) {
-                        int j = i * options->sample; // Index into image
-                        int x = j % image->numCols; // x index
-                        int y = j / image->numCols; // y index
-                        sample->data.F32[i] = image->data.F32[y][x];
-                        if (readout->mask) {
-                            sampleMask->data.U8[i] = mask->data.U8[y][x];
-                        }
-                    }
-                    status = psVectorStats(bgStats, sample, sampleMask, NULL, options->combine->maskVal);
-                    if (!status) {
-                      psTrace("ppMerge", 3, "failed to get stats for for %s, cell %d is %f\n", name, cellNum,
-                              background->data.F32[i][cellNum]);
-                      psErrorClear();
-                    }
-                    psFree(sample);
-                    psFree(sampleMask);
-                    background->data.F32[i][cellNum] = psStatsGetValue(bgStats, options->mean);
-                    psTrace("ppMerge", 3, "Background for %s, cell %d is %f\n", name, cellNum,
-                            background->data.F32[i][cellNum]);
-                }
-
-                // Shutter correction
-                if (options->shutter) {
-                    if (!pmCellRead(cell, inFile, config->database)) {
-                        // Nothing here
-                        pmCellFreeData(cell);
-                        continue;
-                    }
-
-                    if (cell->readouts->n > 1) {
-                        psWarning("File %s chip %d cell %d contains more than one "
-                                  "readout --- ignoring all but the first.\n", name, j, k);
-                        status = false;
-                    }
-                    pmReadout *readout = cell->readouts->data[0]; // The readout of interest
-
-                    pmShutterCorrectionData *shutter = (*shutters)->data[cellNum]; // Shutter correction data
-                    if (!shutter) {
-                        shutter = pmShutterCorrectionDataAlloc(readout->image->numCols,
-                                                               readout->image->numRows,
-                                                               options->shutterSize);
-                        (*shutters)->data[cellNum] = shutter;
-                    }
-                    if (!pmShutterCorrectionAddReadout(shutter, readout, options->mean, options->stdev,
-                                                       options->combine->maskVal, rng)) {
-                        psWarning("Can't add file %s chip %d cell %d to shutter correction --- ignored.",
-                                  name, j, k);
-                        status = false;
-                    }
-                }
-
-
-                pmCellFreeData(cell);
-            }
-            pmChipFreeData(chip);
-        }
-        pmFPAFreeData(fpa);
-    }
+
+    psFree(background);
+    psFree(shutters);
+    return true;
+
+ERROR:
+    // Common path for errors
+    psFree(gains);
+    psFree(background);
+    psFree(shutters);
     psFree(rng);
-    psFree(bgStats);
+    psFree(stats);
     psFree(view);
-
-    if (options->scale) {
-        // Need to normalize over the focal plane
-        if (psTraceGetLevel("ppMerge") > 9) {
-            for (int i = 0; i < gains->n; i++) {
-                psTrace("ppMerge", 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]);
-            }
-        }
-        psVector *fluxes = NULL;        // Solution to fluxes
-        if (!pmFlatNormalize(&fluxes, &gains, background)) {
-            psWarning("Normalisation failed to converge --- continuing anyway.\n");
-            status = false;
-        }
-        psFree(gains);
-
-        psImage *scalesDeref = *scales; // Dereference the pointer
-
-        for (int i = 0; i < scalesDeref->numRows; i++) {
-            psF32 bg = fluxes->data.F32[i];
-            for (int j = 0; j < scalesDeref->numCols; j++) {
-                scalesDeref->data.F32[i][j] = bg;
-            }
-        }
-    }
-
-    if (options->zero) {
-        if (!*zeros) {
-            // This is much faster than copying!
-            *zeros = psMemIncrRefCounter(background);
-        } else {
-            *zeros = psImageCopy(*zeros, background, PS_TYPE_F32);
-        }
-    }
-
-    // Diagnostic stuff
-    if (scales && *scales && psTraceGetLevel("ppMerge") > 9) {
-        psImage *scalesDeref = *scales; // Dereference the pointer
-        for (int i = 0; i < scalesDeref->numRows; i++) {
-            for (int j = 0; j < scalesDeref->numCols; j++) {
-                psTrace("ppMerge", 9, "Scale for exposure %d, cell %d is: %f\n", i, j,
-                        scalesDeref->data.F32[i][j]);
-            }
-        }
-    }
-    if (zeros && *zeros && psTraceGetLevel("ppMerge") > 9) {
-        psImage *zerosDeref = *zeros; // Dereference the pointer
-        for (int i = 0; i < zerosDeref->numRows; i++) {
-            for (int j = 0; j < zerosDeref->numCols; j++) {
-                psTrace("ppMerge", 9, "Zero for exposure %d, cell %d is: %f\n", i, j,
-                        zerosDeref->data.F32[i][j]);
-            }
-        }
-    }
-
-    psFree(background);
-
-    return status;
+    return false;
 }