Changeset 27838 for branches/tap_branches/ppStack/src

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• ppStack/src (modified) (1 prop)
• ppStack/src/Makefile.am (modified) (2 diffs)
• ppStack/src/ppStack.c (modified) (2 diffs)
• ppStack/src/ppStack.h (modified) (6 diffs)
• ppStack/src/ppStackArguments.c (modified) (11 diffs)
• ppStack/src/ppStackCamera.c (modified) (22 diffs)
• ppStack/src/ppStackCleanup.c (modified) (5 diffs)
• ppStack/src/ppStackCombineFinal.c (modified) (4 diffs)
• ppStack/src/ppStackCombineInitial.c (modified) (8 diffs)
• ppStack/src/ppStackCombinePrepare.c (modified) (4 diffs)
• ppStack/src/ppStackConvolve.c (modified) (15 diffs)
• ppStack/src/ppStackErrorCodes.c.in (copied) (copied from trunk/ppStack/src/ppStackErrorCodes.c.in )
• ppStack/src/ppStackErrorCodes.dat (copied) (copied from trunk/ppStack/src/ppStackErrorCodes.dat )
• ppStack/src/ppStackErrorCodes.h.in (copied) (copied from trunk/ppStack/src/ppStackErrorCodes.h.in )
• ppStack/src/ppStackFiles.c (modified) (5 diffs)
• ppStack/src/ppStackFinish.c (modified) (2 diffs)
• ppStack/src/ppStackLoop.c (modified) (7 diffs)
• ppStack/src/ppStackLoop.h (modified) (3 diffs)
• ppStack/src/ppStackMatch.c (modified) (26 diffs)
• ppStack/src/ppStackOptions.c (modified) (4 diffs)
• ppStack/src/ppStackOptions.h (modified) (2 diffs)
• ppStack/src/ppStackPSF.c (modified) (4 diffs)
• ppStack/src/ppStackPhotometry.c (modified) (6 diffs)
• ppStack/src/ppStackPrepare.c (modified) (12 diffs)
• ppStack/src/ppStackReadout.c (modified) (16 diffs)
• ppStack/src/ppStackReject.c (modified) (10 diffs)
• ppStack/src/ppStackSetup.c (modified) (7 diffs)
• ppStack/src/ppStackSources.c (modified) (7 diffs)
• ppStack/src/ppStackThread.c (modified) (10 diffs)
• ppStack/src/ppStackThread.h (modified) (1 diff)

branches/tap_branches/ppStack/src

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

branches/tap_branches/ppStack/src/Makefile.am

@@ -47 +47 @@
         ppStackCleanup.c        \
         ppStackPhotometry.c     \
-        ppStackFinish.c
+        ppStackFinish.c         \
+        ppStackErrorCodes.c
 
 noinst_HEADERS =                \
@@ -53 +54 @@
         ppStackLoop.h           \
         ppStackOptions.h        \
-        ppStackThread.h
+        ppStackThread.h         \
+        ppStackErrorCodes.h
 
-CLEANFILES = *~
+### Error codes.
+BUILT_SOURCES = ppStackErrorCodes.h ppStackErrorCodes.c
+CLEANFILES = ppStackErrorCodes.h ppStackErrorCodes.c
+
+ppStackErrorCodes.h : ppStackErrorCodes.dat ppStackErrorCodes.h.in
+        $(ERRORCODES) --data=ppStackErrorCodes.dat --outdir=. ppStackErrorCodes.h
+
+ppStackErrorCodes.c : ppStackErrorCodes.dat ppStackErrorCodes.c.in ppStackErrorCodes.h
+        $(ERRORCODES) --data=ppStackErrorCodes.dat --outdir=. ppStackErrorCodes.c
+
 
 clean-local:

branches/tap_branches/ppStack/src/ppStack.c

@@ -15 +15 @@
 int main(int argc, char *argv[])
 {
-    psExit exitValue = PS_EXIT_SUCCESS; // Exit value
+    psLibInit(NULL);
     psTimerStart(TIMER_NAME);
     psTimerStart("PPSTACK_STEPS");
-    psLibInit(NULL);
+
+    pmErrorRegister();
+    ppStackErrorRegister();
+    psphotErrorRegister();
 
     pmConfig *config = pmConfigRead(&argc, argv, PPSTACK_RECIPE); // Configuration
+    ppStackOptions *options = NULL;                               // Options for stacking
     if (!config) {
-        psErrorStackPrint(stderr, "Error reading configuration.");
-        exitValue = PS_EXIT_CONFIG_ERROR;
-        goto die;
-    }
-
-    (void) psTraceSetLevel("ppStack", 5);
-
-    if (!ppStackArgumentsSetup(argc, argv, config)) {
-        psErrorStackPrint(stderr, "Error reading arguments.\n");
-        exitValue = PS_EXIT_CONFIG_ERROR;
         goto die;
     }
@@ -38 +32 @@
 
     if (!pmModelClassInit()) {
-        psErrorStackPrint(stderr, "Error initialising model classes.\n");
-        exitValue = PS_EXIT_PROG_ERROR;
+        psError(PPSTACK_ERR_PROG, false, "Unable to initialise model classes.");
         goto die;
     }
 
     if (!psphotInit()) {
-        psErrorStackPrint(stderr, "Error initialising psphot.\n");
-        exitValue = PS_EXIT_PROG_ERROR;
+        psError(PPSTACK_ERR_PROG, false, "Error initialising psphot.");
+        goto die;
+    }
+
+    (void)psTraceSetLevel("ppStack", 5);
+
+    if (!ppStackArgumentsSetup(argc, argv, config)) {
         goto die;
     }
 
     if (!ppStackCamera(config)) {
-        psErrorStackPrint(stderr, "Error setting up input files.\n");
-        exitValue = PS_EXIT_CONFIG_ERROR;
         goto die;
     }
 
     if (!ppStackArgumentsParse(config)) {
-        psErrorStackPrint(stderr, "Error reading arguments.\n");
-        exitValue = PS_EXIT_CONFIG_ERROR;
         goto die;
     }
 
-    if (!ppStackLoop(config)) {
-        psErrorStackPrint(stderr, "Error performing combination.\n");
-        exitValue = PS_EXIT_DATA_ERROR;
+    options = ppStackOptionsAlloc();
+    if (!ppStackLoop(config, options)) {
         goto die;
     }
 
 
-     // Common code for the death.
-die:
-    psTrace("ppStack", 1, "Finished at %f sec\n", psTimerMark(TIMER_NAME));
-    psTimerStop();
+ die:
+    // Common code for the death.
+    {
+        psExit exitValue = ppStackExitCode(PS_EXIT_SUCCESS); // Exit code
 
-    psFree(config);
-    pmModelClassCleanup();
-    pmConfigDone();
-    psLibFinalize();
-    pmVisualClose();
-    exit(exitValue);
+        // Ensure everything closes
+        if (config) {
+            ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
+            ppStackFileActivation(config, PPSTACK_FILES_CONVOLVE, true);
+            ppStackFileActivation(config, PPSTACK_FILES_STACK, true);
+            ppStackFileActivation(config, PPSTACK_FILES_UNCONV, true);
+            ppStackFileActivation(config, PPSTACK_FILES_PHOT, true);
+            if (!ppStackFilesIterateUp(config)) {
+                psError(psErrorCodeLast(), false, "Unable to close files.");
+                exitValue = ppStackExitCode(exitValue);
+                pmFPAfileFreeSetStrict(false);
+            }
+        }
+
+        // Write out summary statistics
+        if (options && options->stats) {
+
+            psMetadataAddS32(options->stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
+                             "Bad data quality flag", options->quality);
+            psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_STACK", 0,
+                             "Total time", psTimerClear("PPSTACK_TOTAL"));
+
+            const char *statsMDC = psMetadataConfigFormat(options->stats);
+            if (!statsMDC || strlen(statsMDC) == 0) {
+                psError(PS_ERR_IO, false, "Unable to get statistics MDC file.");
+                return false;
+            }
+            if (fprintf(options->statsFile, "%s", statsMDC) != strlen(statsMDC)) {
+                psError(PS_ERR_IO, false, "Unable to write statistics MDC file.");
+                return false;
+            }
+            psFree(statsMDC);
+            if (fclose(options->statsFile) == EOF) {
+                psError(PS_ERR_IO, false, "Unable to close statistics MDC file.");
+                return false;
+            }
+            options->statsFile = NULL;
+            pmConfigRunFilenameAddWrite(config, "STATS",
+                                        psMetadataLookupStr(NULL, config->arguments, "STATS"));
+        }
+        psFree(options);
+
+        // Dump configuration
+        bool mdok;                                                                    // Status of MD lookup
+        psString dump = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); // File for config
+        if (dump && !pmConfigDump(config, dump)) {
+            psError(psErrorCodeLast(), false, "Unable to dump configuration.");
+            exitValue = ppStackExitCode(exitValue);
+        }
+
+        psTrace("ppStack", 1, "Finished at %f sec\n", psTimerMark(TIMER_NAME));
+        psTimerStop();
+
+        psFree(config);
+        pmModelClassCleanup();
+        pmConfigDone();
+        psLibFinalize();
+        pmVisualClose();
+
+        exitValue = ppStackExitCode(exitValue);
+        exit(exitValue);
+    }
 }
 

branches/tap_branches/ppStack/src/ppStack.h

@@ -9 +9 @@
 
 #include "ppStackOptions.h"
+#include "ppStackErrorCodes.h"
 
 // Mask values for inputs
 typedef enum {
+    PPSTACK_MASK_NONE   = 0x00,         // Nothing wrong
     PPSTACK_MASK_CAL    = 0x01,         // Photometric calibration failed
-    PPSTACK_MASK_MATCH  = 0x02,         // PSF-matching failed
-    PPSTACK_MASK_CHI2   = 0x04,         // Chi^2 too deviant
-    PPSTACK_MASK_REJECT = 0x08,         // Rejection failed
-    PPSTACK_MASK_BAD    = 0x10,         // Bad image (too many pixels rejected)
+    PPSTACK_MASK_PSF    = 0x02,         // PSF measurement failed
+    PPSTACK_MASK_MATCH  = 0x04,         // PSF-matching failed
+    PPSTACK_MASK_CHI2   = 0x08,         // Chi^2 too deviant
+    PPSTACK_MASK_REJECT = 0x10,         // Rejection failed
+    PPSTACK_MASK_BAD    = 0x20,         // Bad image (too many pixels rejected)
     PPSTACK_MASK_ALL    = 0xff          // All errors
 } ppStackMask;
@@ -24 +27 @@
     PPSTACK_FILES_PREPARE,              // Files for preparation
     PPSTACK_FILES_CONVOLVE,             // Files for convolution
-    PPSTACK_FILES_COMBINE,              // Files for combination
+    PPSTACK_FILES_STACK,                // Stack files
+    PPSTACK_FILES_UNCONV,               // Unconvolved stack files
     PPSTACK_FILES_PHOT                  // Files for photometry
 } ppStackFileList;
@@ -59 +63 @@
 // Perform stacking on a readout
 //
-// Returns an array of pixels to inspect for each input image
+// Returns two arrays: pixels to inspect for each input image, and pixels to reject for each input image.
 psArray *ppStackReadoutInitial(const pmConfig *config,   // Configuration
                                pmReadout *outRO,   // Output readout
@@ -65 +69 @@
                                const psVector *mask, // Mask for input readouts
                                const psVector *weightings, // Weighting factors for each image
+                               const psVector *exposures,  // Exposure time for each image
                                const psVector *addVariance // Additional variance for rejection
     );
@@ -79 +84 @@
 bool ppStackReadoutFinal(const pmConfig *config,   // Configuration
                          pmReadout *outRO,   // Output readout
+                         pmReadout *expRO,   // Exposure readout
                          const psArray *readouts, // Input readouts
                          const psVector *mask, // Mask for input readouts
                          const psArray *rejected, // Array with pixels rejected in each image
                          const psVector *weightings, // Weighting factors for each image
-                         const psVector *addVariance // Additional variance for rejection
+                         const psVector *exposures,  // Exposure times for each image
+                         const psVector *addVariance, // Additional variance for rejection
+                         bool safety,                 // Enable safety switch?
+                         const psVector *norm         // Normalisations to apply
     );
 
@@ -164 +173 @@
     );
 
+/// Return an appropriate exit code based on the error code
+psExit ppStackExitCode(psExit exitValue);
 
 #endif

branches/tap_branches/ppStack/src/ppStackArguments.c

@@ -43 +43 @@
         value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \
         if (!mdok) { \
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
+            psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s", \
                 RECIPENAME, PPSTACK_RECIPE); \
             goto ERROR; \
@@ -58 +58 @@
         value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \
         if (!mdok) { \
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
+            psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s", \
                 RECIPENAME, PPSTACK_RECIPE); \
             goto ERROR; \
@@ -73 +73 @@
         name = psMetadataLookupStr(NULL, recipe, RECIPENAME); \
         if (!name) { \
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
-                RECIPENAME, PPSTACK_RECIPE); \
+            psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s", \
+                    RECIPENAME, PPSTACK_RECIPE);                        \
             goto ERROR; \
         } \
@@ -108 +108 @@
         value = psMetadataLookupStr(NULL, recipe, mdName);
         if (!value) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s",
+            psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s",
                     mdName, PPSTACK_RECIPE);
             return false;
@@ -148 +148 @@
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
-    psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-iter", 0, "Number of rejection iterations per input", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0,
                      "Combination rejection thresold (sigma)", NAN);
@@ -185 +185 @@
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-1", 0, "Estimated systematic error; pass 1", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-2", 0, "Estimated systematic error; pass 2", NAN);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-dir", 0, "Directory for temporary images", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-image", 0, "Suffix for temporary images", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-mask", 0, "Suffix for temporary masks", NULL);
@@ -206 +207 @@
         psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[argNum], false); // Input file info
         if (!inputs || numBad > 0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to cleanly read MDC file with inputs.");
+            psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to cleanly read MDC file with inputs.");
             return false;
         }
@@ -228 +229 @@
     int numThreads = psMetadataLookupS32(NULL, arguments, "-threads"); // Number of threads
     if (numThreads > 0 && !psThreadPoolInit(numThreads)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", numThreads);
+        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to setup %d threads", numThreads);
         return false;
     }
@@ -246 +247 @@
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe
     if (!recipe) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE);
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find recipe %s", PPSTACK_RECIPE);
         goto ERROR;
     }
 
-    VALUE_ARG_RECIPE_INT("-iter",              "ITER",            S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-combine-iter",    "COMBINE.ITER",    F32);
     VALUE_ARG_RECIPE_FLOAT("-combine-rej",     "COMBINE.REJ",     F32);
     VALUE_ARG_RECIPE_FLOAT("-combine-sys",     "COMBINE.SYS",     F32);
@@ -260 +261 @@
     VALUE_ARG_RECIPE_FLOAT("-poor-frac",       "POOR.FRACTION",   F32);
 
-    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.VAL",  recipe);
+    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.IN",   recipe);
     valueArgRecipeStr(arguments, recipe, "-mask-bad",  "MASK.BAD",  recipe);
     valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe);
@@ -316 +317 @@
 
     pmConfigCamerasCull(config, NULL);
-    pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,MASKS,JPEG");
+    pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
 
     return true;

branches/tap_branches/ppStack/src/ppStackCamera.c

@@ -30 +30 @@
         pmFPAfileDefineFromArgs(&found, config, name, "FILENAMES");
     if (!file || !found) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to define file %s from %s", name, filename);
+        psError(psErrorCodeLast(), false, "Unable to define file %s from %s", name, filename);
         return NULL;
     }
     if (file->type != type) {
-        psError(PS_ERR_IO, true, "%s is not of type %s", name, pmFPAfileStringFromType(type));
+        psError(PS_ERR_IO, PPSTACK_ERR_CONFIG, "%s is not of type %s", name, pmFPAfileStringFromType(type));
         return NULL;
     }
@@ -53 +53 @@
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppSim
     if (!recipe) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE);
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find recipe %s", PPSTACK_RECIPE);
         return false;
     }
@@ -67 +67 @@
                                                            "PPSTACK.INPUT.MASK"); // Input masks
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define input masks from RUN metadata.");
+            psError(psErrorCodeLast(), false, "Unable to define input masks from RUN metadata.");
             psFree(runImages);
             return false;
@@ -76 +76 @@
                                                           "PPSTACK.INPUT.VARIANCE"); // Input variances
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define input variances from RUN metadata.");
+            psError(psErrorCodeLast(), false, "Unable to define input variances from RUN metadata.");
             psFree(runImages);
             return false;
@@ -88 +88 @@
                                                          "PPSTACK.INPUT.SOURCES"); // Input sources
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define input sources from RUN metadata.");
+            psError(psErrorCodeLast(), false, "Unable to define input sources from RUN metadata.");
             psFree(runImages);
             return false;
         }
         if (!runSrc) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to define input sources from RUN metadata.");
+            psError(PPSTACK_ERR_CONFIG, true, "Unable to define input sources from RUN metadata.");
             psFree(runImages);
             return false;
@@ -104 +104 @@
                                                          "PPSTACK.INPUT.PSF"); // Input PSFs
                 if (!status) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to define input PSFs from RUN metadata.");
+                    psError(psErrorCodeLast(), false, "Unable to define input PSFs from RUN metadata.");
                     psFree(runImages);
                     return false;
@@ -118 +118 @@
                                                                     "PPSTACK.CONV.KERNEL"); // Conv'n kernels
                 if (!status) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to define convolution kernels from RUN metadata.");
+                    psError(psErrorCodeLast(), false,
+                            "Unable to define convolution kernels from RUN metadata.");
                     psFree(runImages);
                     return false;
                 }
                 if (!runKernel) {
-                    psError(PS_ERR_UNEXPECTED_NULL, true,
+                    psError(PPSTACK_ERR_CONFIG, true,
                             "Unable to define convolution kernels from RUN metadata.");
                     psFree(runImages);
@@ -137 +138 @@
         psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info
         if (!inputs) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find inputs.");
+            psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to find inputs.");
             return false;
         }
@@ -145 +146 @@
         while ((item = psMetadataGetAndIncrement(iter))) {
             if (item->type != PS_DATA_METADATA) {
-                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                psError(PPSTACK_ERR_ARGUMENTS, true,
                         "Component %s of the input metadata is not of type METADATA", item->name);
                 psFree(iter);
@@ -155 +156 @@
             psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image
             if (!image || strlen(image) == 0) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name);
+                psError(PPSTACK_ERR_ARGUMENTS, true, "Component %s lacks IMAGE of type STR", item->name);
                 psFree(iter);
                 return false;
@@ -169 +170 @@
                                               image, PM_FPA_FILE_IMAGE); // File for image
             if (!imageFile) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
+                psError(psErrorCodeLast(), false, "Unable to define file from image %d (%s)", i, image);
                 return false;
             }
@@ -175 +176 @@
             if (mask && strlen(mask) > 0 &&
                 !defineFile(config, imageFile, "PPSTACK.INPUT.MASK", mask, PM_FPA_FILE_MASK)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
+                psError(psErrorCodeLast(), false, "Unable to define file from mask %d (%s)", i, mask);
                 return false;
             }
@@ -183 +184 @@
                 if (!defineFile(config, imageFile, "PPSTACK.INPUT.VARIANCE", variance,
                                 PM_FPA_FILE_VARIANCE)) {
-                    psError(PS_ERR_UNKNOWN, false,
+                    psError(psErrorCodeLast(), false,
                             "Unable to define file from variance %d (%s)", i, variance);
                     return false;
@@ -195 +196 @@
                     havePSFs = true;
                     if (!defineFile(config, imageFile, "PPSTACK.INPUT.PSF", psf, PM_FPA_FILE_PSF)) {
-                        psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
+                        psError(psErrorCodeLast(), false, "Unable to define file from psf %d (%s)", i, psf);
                         return false;
                     }
                 }
             } else if (havePSFs) {
-                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF %d", i);
+                psError(PPSTACK_ERR_CONFIG, true, "Unable to find PSF %d", i);
                 return false;
             }
 
             if (!sources || strlen(sources) == 0) {
-                psError(PS_ERR_UNEXPECTED_NULL, true, "SOURCES not provided for file %d", i);
+                psError(PPSTACK_ERR_CONFIG, true, "SOURCES not provided for file %d", i);
                 return false;
             }
             if (!defineFile(config, imageFile, "PPSTACK.INPUT.SOURCES", sources, PM_FPA_FILE_CMF)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)",
+                psError(psErrorCodeLast(), false, "Unable to define file from sources %d (%s)",
                         i, sources);
                 return false;
@@ -217 +218 @@
                 pmFPAfile *kernel = pmFPAfileDefineOutput(config, imageFile->fpa, "PPSTACK.CONV.KERNEL");
                 if (!kernel) {
-                    psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.CONV.KERNEL"));
+                    psError(psErrorCodeLast(), false,
+                            "Unable to generate output file from PPSTACK.CONV.KERNEL");
                     return false;
                 }
@@ -233 +235 @@
 
     // Output image
-    pmFPA *fpa = pmFPAConstruct(config->camera, config->cameraName); // 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
+    pmFPA *outFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output
+    if (!outFPA) {
+        psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration.");
+        return false;
+    }
+    pmFPAfile *output = pmFPAfileDefineOutput(config, outFPA, "PPSTACK.OUTPUT");
+    psFree(outFPA);                        // Drop reference
     if (!output) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT"));
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT"));
         return false;
     }
     if (output->type != PM_FPA_FILE_IMAGE) {
-        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT is not of type IMAGE");
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT is not of type IMAGE");
         return false;
     }
     output->save = true;
 
-    if (!pmFPAAddSourceFromFormat(fpa, "Stack", output->format)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA.");
-        psFree(fpa);
+    if (!pmFPAAddSourceFromFormat(outFPA, "Stack", output->format)) {
+        psError(psErrorCodeLast(), false, "Unable to generate output FPA.");
         return false;
     }
@@ -259 +260 @@
     pmFPAfile *outMask = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.MASK");
     if (!outMask) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.MASK"));
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.MASK"));
         return false;
     }
     if (outMask->type != PM_FPA_FILE_MASK) {
-        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.MASK is not of type MASK");
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.MASK is not of type MASK");
         return false;
     }
@@ -272 +273 @@
         pmFPAfile *outVariance = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.VARIANCE");
         if (!outVariance) {
-            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.VARIANCE"));
+            psError(psErrorCodeLast(), false, "Unable to generate output file from PPSTACK.OUTPUT.VARIANCE");
             return false;
         }
         if (outVariance->type != PM_FPA_FILE_VARIANCE) {
-            psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.VARIANCE is not of type VARIANCE");
+            psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.VARIANCE is not of type VARIANCE");
             return false;
         }
@@ -282 +283 @@
     }
 
+
+    // Exposure image
+    pmFPA *expFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output
+    if (!expFPA) {
+        psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration.");
+        return false;
+    }
+    pmFPAfile *exp = pmFPAfileDefineOutput(config, expFPA, "PPSTACK.OUTPUT.EXP");
+    psFree(expFPA);                        // Drop reference
+    if (!exp) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.EXP"));
+        return false;
+    }
+    if (exp->type != PM_FPA_FILE_IMAGE) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXP is not of type IMAGE");
+        return false;
+    }
+    exp->save = true;
+
+    if (!pmFPAAddSourceFromFormat(expFPA, "Stack", exp->format)) {
+        psError(psErrorCodeLast(), false, "Unable to generate output FPA.");
+        return false;
+    }
+
+    // Exposure numbers
+    pmFPAfile *expNum = pmFPAfileDefineOutput(config, exp->fpa, "PPSTACK.OUTPUT.EXPNUM");
+    if (!expNum) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.EXPNUM"));
+        return false;
+    }
+    if (expNum->type != PM_FPA_FILE_MASK) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXPNUM is not of type MASK");
+        return false;
+    }
+    expNum->save = true;
+
+    // Weighted exposure
+    pmFPAfile *expWt = pmFPAfileDefineOutput(config, exp->fpa, "PPSTACK.OUTPUT.EXPWT");
+    if (!expWt) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.EXPWT"));
+        return false;
+    }
+    if (expWt->type != PM_FPA_FILE_VARIANCE) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXPWT is not of type VARIANCE");
+        return false;
+    }
+    expWt->save = true;
+
+
     if (havePSFs) {
-        pmFPAfile *targetPSF = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.TARGET.PSF");
+        pmFPA *psfFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain PSF
+        if (!psfFPA) {
+            psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration.");
+            return false;
+        }
+        pmFPAfile *targetPSF = pmFPAfileDefineOutput(config, psfFPA, "PPSTACK.TARGET.PSF");
         if (!targetPSF) {
-            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.TARGET.PSF"));
-            return false;
-        }
+            psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.TARGET.PSF"));
+            return false;
+        }
+        psFree(psfFPA);
         if (targetPSF->type != PM_FPA_FILE_PSF) {
-            psError(PS_ERR_IO, true, "PPSTACK.TARGET.PSF is not of type PSF");
+            psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.TARGET.PSF is not of type PSF");
             return false;
         }
         targetPSF->save = true;
     }
+
+    // Unconvolved stack
+    pmFPA *unconvFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain unconvolved output
+    if (!unconvFPA) {
+        psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration.");
+        return false;
+    }
+    pmFPAfile *unConv = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV");
+    psFree(unconvFPA);                  // Drop reference
+    if (!unConv) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV"));
+        return false;
+    }
+    if (unConv->type != PM_FPA_FILE_IMAGE) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV is not of type IMAGE");
+        return false;
+    }
+    unConv->save = true;
+
+    if (!pmFPAAddSourceFromFormat(unconvFPA, "Stack", unConv->format)) {
+        psError(psErrorCodeLast(), false, "Unable to generate output FPA.");
+        return false;
+    }
+
+    // Unconvolved mask
+    pmFPAfile *unconvMask = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV.MASK");
+    if (!unconvMask) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.MASK"));
+        return false;
+    }
+    if (unconvMask->type != PM_FPA_FILE_MASK) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.MASK is not of type MASK");
+        return false;
+    }
+    unconvMask->save = true;
+
+    // Unconvolved variance
+    if (haveVariances) {
+        pmFPAfile *unconvVariance = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV.VARIANCE");
+        if (!unconvVariance) {
+            psError(psErrorCodeLast(), false, "Unable to generate output file from PPSTACK.UNCONV.VARIANCE");
+            return false;
+        }
+        if (unconvVariance->type != PM_FPA_FILE_VARIANCE) {
+            psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.VARIANCE is not of type VARIANCE");
+            return false;
+        }
+        unconvVariance->save = true;
+    }
+
+
+    // Exposure image
+    pmFPA *unconvExpFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output
+    if (!unconvExpFPA) {
+        psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration.");
+        return false;
+    }
+    pmFPAfile *unconvExp = pmFPAfileDefineOutput(config, unconvExpFPA, "PPSTACK.UNCONV.EXP");
+    psFree(unconvExpFPA);               // Drop reference
+    if (!unconvExp) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.EXP"));
+        return false;
+    }
+    if (unconvExp->type != PM_FPA_FILE_IMAGE) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXP is not of type IMAGE");
+        return false;
+    }
+    unconvExp->save = true;
+
+    if (!pmFPAAddSourceFromFormat(unconvExpFPA, "Stack", unconvExp->format)) {
+        psError(psErrorCodeLast(), false, "Unable to generate output FPA.");
+        return false;
+    }
+
+    // Exposure numbers
+    pmFPAfile *unconvExpNum = pmFPAfileDefineOutput(config, unconvExp->fpa, "PPSTACK.UNCONV.EXPNUM");
+    if (!unconvExpNum) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.MASK"));
+        return false;
+    }
+    if (unconvExpNum->type != PM_FPA_FILE_MASK) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXPNUM is not of type MASK");
+        return false;
+    }
+    unconvExpNum->save = true;
+
+    // Weighted exposure
+    pmFPAfile *unconvExpWt = pmFPAfileDefineOutput(config, unconvExp->fpa, "PPSTACK.UNCONV.EXPWT");
+    if (!unconvExpWt) {
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.EXPWT"));
+        return false;
+    }
+    if (unconvExpWt->type != PM_FPA_FILE_VARIANCE) {
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXPWT is not of type VARIANCE");
+        return false;
+    }
+    unconvExpWt->save = true;
 
     // Output JPEGs
     pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT.JPEG1");
     if (!jpeg1) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG1"));
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG1"));
         return false;
     }
     if (jpeg1->type != PM_FPA_FILE_JPEG) {
-        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.JPEG1 is not of type JPEG");
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.JPEG1 is not of type JPEG");
         return false;
     }
@@ -308 +461 @@
     pmFPAfile *jpeg2 = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT.JPEG2");
     if (!jpeg2) {
-        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG2"));
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG2"));
         return false;
     }
     if (jpeg2->type != PM_FPA_FILE_JPEG) {
-        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.JPEG2 is not of type JPEG");
+        psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.JPEG2 is not of type JPEG");
         return false;
     }
@@ -323 +476 @@
         psMetadataLookupBool(&mdok, config->arguments, "-photometry"); // perform photometry
     if (doPhotom) {
-        // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by
-        // psphotDefineFiles
+        // This pmFPAfile, PSPHOT.INPUT, is just used as a carrier; output files (eg,
+        // PSPHOT.RESID) are defined by psphotDefineFiles
         pmFPAfile *psphotInput = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
         if (!psphotInput) {
-            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.INPUT"));
-            return false;
-        }
+            psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.INPUT"));
+            return false;
+        }
+        // specify the number of psphot input images
+        psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1);
 
         // Define associated psphot input/output files
         if (!psphotDefineFiles(config, psphotInput)) {
-            psError(PSPHOT_ERR_CONFIG, false,
+            psError(psErrorCodeLast(), false,
                     "Trouble defining the additional input/output files for psphot");
             return false;
@@ -341 +496 @@
         pmFPAfile *outPSF = pmFPAfileDefineOutputFromFile(config, output, "PSPHOT.PSF.SAVE");
         if (!outPSF) {
-            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.PSF.SAVE"));
+            psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.PSF.SAVE"));
             return false;
         }
         if (outPSF->type != PM_FPA_FILE_PSF) {
-            psError(PS_ERR_IO, true, "PSPHOT.PSF.SAVE is not of type PSF");
+            psError(PPSTACK_ERR_CONFIG, true, "PSPHOT.PSF.SAVE is not of type PSF");
             return false;
         }

branches/tap_branches/ppStack/src/ppStackCleanup.c

@@ -4 +4 @@
 
 #include <stdio.h>
-#include <unistd.h>
 #include <pslib.h>
 #include <psmodules.h>
+#include <ppStats.h>
 
 #include "ppStack.h"
 #include "ppStackLoop.h"
-
-#define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
 
@@ -22 +20 @@
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
     psAssert(recipe, "We've thrown an error on this before.");
-
-    bool mdok;                          // Status of MD lookup
-    bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files?
 
 #if 0
@@ -37 +32 @@
 #endif
 
-    // Close up
-    bool wcsDone = false;           // Have we done the WCS?
-    for (int i = 0; i < options->num; i++) {
-        if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            continue;
-        }
-
-        ppStackThread *thread = stack->threads->data[0]; // Representative stack
-        pmReadout *inRO = thread->readouts->data[i]; // Template readout
-        if (inRO && !wcsDone) {
-            // Copy astrometry over
-            wcsDone = true;
-            pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
-            pmHDU *outHDU = pmHDUFromCell(options->outRO->parent); // Output HDU
-            pmChip *outChip = options->outRO->parent->parent; // Output chip
-            pmFPA *outFPA = outChip->parent; // Output FPA
-            if (!outHDU || !inHDU) {
-                psWarning("Unable to find HDU at FPA level to copy astrometry.");
-            } else {
-                if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
-                    psErrorClear();
-                    psWarning("Unable to read WCS astrometry from input FPA.");
-                    wcsDone = false;
-                } else {
-                    if (!outHDU->header) {
-                        outHDU->header = psMetadataAlloc();
-                    }
-                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
-                        psErrorClear();
-                        psWarning("Unable to write WCS astrometry to output FPA.");
-                        wcsDone = false;
-                    }
-                }
-            }
-        }
-
-        if (tempDelete) {
-            psString imageResolved = pmConfigConvertFilename(options->imageNames->data[i],
-                                                             config, false, false);
-            psString maskResolved = pmConfigConvertFilename(options->maskNames->data[i],
-                                                            config, false, false);
-            psString varianceResolved = pmConfigConvertFilename(options->varianceNames->data[i],
-                                                                config, false, false);
-            if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
-                unlink(varianceResolved) == -1) {
-                psWarning("Unable to delete temporary files for image %d", i);
-            }
-            psFree(imageResolved);
-            psFree(maskResolved);
-            psFree(varianceResolved);
-        }
-    }
-
-    ppStackMemDump("cleanup");
-
     // Generate binned JPEGs
     {
@@ -103 +43 @@
         pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2");
         psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
+        psFree(view);
 
         pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
         if (!pmReadoutRebin(ro1, options->outRO, maskValue, bin1, bin1) ||
             !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to bin output.");
+            psError(PPSTACK_ERR_DATA, false, "Unable to bin output.");
             psFree(ro1);
             psFree(ro2);
@@ -116 +57 @@
     }
 
-    // Put version information into the header
-    pmHDU *hdu = pmHDUFromCell(options->outRO->parent);
-    if (!hdu) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
+    // Statistics on output
+    if (options->stats) {
+        psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
+        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
+        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+        pmFPAview *view = pmFPAviewAlloc(0); // View to readout
+        view->chip = view->cell = view->readout = 0;
+
+        ppStatsFPA(options->stats, options->outRO->parent->parent->parent, view, maskBad, config);
+
+        psFree(view);
+    }
+
+    if (!ppStackFilesIterateUp(config)) {
+        psError(psErrorCodeLast(), false, "Unable to close files.");
         return false;
     }
-    if (!hdu->header) {
-        hdu->header = psMetadataAlloc();
+    ppStackFileActivation(config, PPSTACK_FILES_STACK, false);
+    ppStackFileActivation(config, PPSTACK_FILES_PHOT, false);
+
+    // Ensure files are freed
+    {
+        options->outRO->data_exists = false;
+        options->outRO->parent->data_exists = false;
+        options->outRO->parent->parent->data_exists = false;
+        psFree(options->outRO);
+        options->outRO = NULL;
+
+        options->expRO->data_exists = false;
+        options->expRO->parent->data_exists = false;
+        options->expRO->parent->parent->data_exists = false;
+        psFree(options->expRO);
+        options->expRO = NULL;
+
+        pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
+        view->chip = view->cell = 0;        // pmFPAviewFreeData doesn't want to deal with readouts
+        pmFPAfile *jpeg1 = pmFPAfileSelectSingle(config->files, "PPSTACK.OUTPUT.JPEG1", 0); // JPEG file
+        pmFPAviewFreeData(view, jpeg1);
+        pmFPAfile *jpeg2 = pmFPAfileSelectSingle(config->files, "PPSTACK.OUTPUT.JPEG2", 0); // JPEG file
+        pmFPAviewFreeData(view, jpeg2);
+        pmFPAfile *phot = NULL;         // Photometry file
+        if (options->photometry) {
+            phot = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // Photometry file
+            pmFPAviewFreeData(view, phot);
+        }
+
+        view->readout = 0;
+        pmReadout *ro1 = pmFPAviewThisReadout(view, jpeg1->fpa); // JPEG readout
+        ro1->data_exists = ro1->parent->data_exists = ro1->parent->parent->data_exists = false;
+        pmReadout *ro2 = pmFPAviewThisReadout(view, jpeg2->fpa); // JPEG readout
+        ro2->data_exists = ro2->parent->data_exists = ro2->parent->parent->data_exists = false;
+        if (options->photometry) {
+            pmReadout *ro = pmFPAviewThisReadout(view, phot->fpa); // Photometry readout
+            ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false;
+        }
+        psFree(view);
     }
-    ppStackVersionHeader(hdu->header);
 
     return true;

branches/tap_branches/ppStack/src/ppStackCombineFinal.c

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

branches/tap_branches/ppStack/src/ppStackCombineInitial.c

@@ -10 +10 @@
 #include "ppStackLoop.h"
 
+//#define TESTING                         // Enable test output
+
 bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
 {
@@ -18 +20 @@
     if (!options->convolve) {
         // No need to do initial combination when we haven't convolved
+        // XXX either allocate inspect and rejected here, or do not require them downstream
         return true;
     }
@@ -34 +37 @@
         if (!status) {
             // Something went wrong
-            psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
+            psError(psErrorCodeLast(), false, "Unable to read chunk %d", numChunk);
             return false;
         }
@@ -54 +57 @@
 
     if (!psThreadPoolWait(false)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination.");
+        psError(psErrorCodeLast(), false, "Unable to do initial combination.");
         return false;
     }
@@ -60 +63 @@
     // Harvest the jobs, gathering the inspection lists
     options->inspect = psArrayAlloc(options->num);
+    options->rejected = psArrayAlloc(options->num);
     for (int i = 0; i < options->num; i++) {
         if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
@@ -65 +69 @@
         }
         options->inspect->data[i] = psArrayAllocEmpty(numChunk);
+        options->rejected->data[i] = psArrayAllocEmpty(numChunk);
     }
     psThreadJob *job;               // Completed job
@@ -71 +76 @@
                  "Job has incorrect type: %s", job->type);
         psArray *results = job->results; // Results of job
+        psAssert(results->n == 2, "Results array has wrong size!");
+        psArray *inspect = results->data[0]; // Pixels to inspect
+        psArray *reject = results->data[1];  // Pixels to reject
         for (int i = 0; i < options->num; i++) {
             if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
                 continue;
             }
-            options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, results->data[i]);
+            options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, inspect->data[i]);
+            options->rejected->data[i] = psArrayAdd(options->rejected->data[i], 1, reject->data[i]);
         }
         psFree(job);
@@ -83 +92 @@
 
 #ifdef TESTING
-    ppStackWriteImage("combined_initial.fits", NULL, outRO->image, config);
-    pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits");
+    ppStackWriteImage("combined_image_initial.fits", NULL, options->outRO->image, config);
+    ppStackWriteImage("combined_mask_initial.fits", NULL, options->outRO->mask, config);
+    ppStackWriteImage("combined_variance_initial.fits", NULL, options->outRO->variance, config);
+
+    pmStackVisualPlotTestImage(options->outRO->image, "combined_image_initial.fits");
 #endif
-
-    psFree(options->matchChi2); options->matchChi2 = NULL;
-
 
     if (options->stats) {

branches/tap_branches/ppStack/src/ppStackCombinePrepare.c

@@ -10 +10 @@
 #include "ppStackLoop.h"
 
-bool ppStackCombinePrepare(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
+bool ppStackCombinePrepare(const char *outName, const char *expName,
+                           ppStackFileList files, ppStackThreadData *stack,
+                           ppStackOptions *options, pmConfig *config)
 {
     psAssert(stack, "Require stack");
@@ -20 +22 @@
     ppStackThread *thread = stack->threads->data[0]; // Representative thread
     if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
-        psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
+        psError(PPSTACK_ERR_ARGUMENTS, false, "problem setting output readout size.");
         return false;
     }
 
     pmFPAfileActivate(config->files, false, NULL);
-    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true);
+    ppStackFileActivation(config, files, true);
     pmFPAview *view = ppStackFilesIterateDown(config); // View to readout
     if (!view) {
@@ -31 +33 @@
     }
 
-    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
-    options->outRO = pmReadoutAlloc(outCell); // Output readout
+    pmCell *cell = pmFPAfileThisCell(config->files, view, outName); // Output cell
+    options->outRO = pmReadoutAlloc(cell); // Output readout
+
+    pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell
+    options->expRO = pmReadoutAlloc(expCell); // Output readout
+
     psFree(view);
 
@@ -39 +45 @@
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
     if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
+        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
+        return false;
+    }
+
+    if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) {
+        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
         return false;
     }

branches/tap_branches/ppStack/src/ppStackConvolve.c

@@ -9 +9 @@
 #include "ppStack.h"
 #include "ppStackLoop.h"
+
+//#define TESTING
+
 
 // Update the value of a concept
@@ -39 +42 @@
     options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
     psVectorInit(options->weightings, 0.0);
-    options->covariances = psArrayAlloc(num); // Covariance matrices
+    options->origCovars = psArrayAlloc(num);
+    options->convCovars = psArrayAlloc(num); // Covariance matrices
+
+    psVector *renorms = psVectorAlloc(num, PS_TYPE_F32); // Renormalisation values for variances
+    psVectorInit(renorms, NAN);
 
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
@@ -52 +59 @@
         pmFPAfileActivate(config->files, false, NULL);
         ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i);
+        if (options->convolve) {
+            // XXX PPSTACK.CONV.KERNEL not defined unless convolve
+            // pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
+            pmFPAfileActivateSingle(config->files, true, "PPSTACK.CONV.KERNEL", i); // Activated file
+        }
+
         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
         pmFPAview *view = ppStackFilesIterateDown(config);
@@ -69 +82 @@
         } else if (options->numCols != readout->image->numCols ||
                    options->numRows != readout->image->numRows) {
-            psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d",
+            psError(PPSTACK_ERR_ARGUMENTS, true, "Sizes of input images don't match: %dx%d vs %dx%d",
                     readout->image->numCols, readout->image->numRows, options->numCols, options->numRows);
             psFree(rng);
@@ -79 +92 @@
         // Background subtraction, scaling and normalisation is performed automatically by the image matching
         psTimerStart("PPSTACK_MATCH");
+        options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance);
         if (!ppStackMatch(readout, options, i, config)) {
-            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
-            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH;
-            psErrorClear();
-            continue;
-        }
-        options->covariances->data[i] = psMemIncrRefCounter(readout->covariance);
+            // XXX many things can cause a failure of ppStackMatch -- should some be handled differently?
+            psErrorCode error = psErrorCodeLast(); // Error code
+            switch (error) {
+                // Fatal errors
+              case PM_ERR_CONFIG:
+              case PPSTACK_ERR_CONFIG:
+              case PPSTACK_ERR_IO:
+                psError(error, false, "Unable to match image %d due to fatal error.", i);
+                return false;
+                // Non-fatal errors
+              case PM_ERR_STAMPS:
+              case PM_ERR_SMALL_AREA:
+              case PPSTACK_ERR_DATA:
+              default:
+                psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
+                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH;
+                psErrorClear();
+                continue;
+            }
+        }
+        options->convCovars->data[i] = psMemIncrRefCounter(readout->covariance);
+
+        float renorm = psMetadataLookupF32(NULL, readout->analysis, PM_READOUT_ANALYSIS_RENORM);
+        if (!isfinite(renorm)) {
+            renorm = 1.0;
+        }
+        renorms->data.F32[i] = renorm;
 
         if (options->stats) {
@@ -114 +149 @@
         pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
         assert(hdu);
-        ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config);
+        if (!ppStackWriteImage(options->convImages->data[i], hdu->header, readout->image, config)) {
+            psError(PPSTACK_ERR_IO, false, "Unable to write convolved image %d", i);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(rng);
+            return false;
+        }
         psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
         pmConfigMaskWriteHeader(config, maskHeader);
-        ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config);
+        if (!ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config)) {
+            psError(PPSTACK_ERR_IO, false, "Unable to write convolved mask %d", i);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(rng);
+            psFree(maskHeader);
+            return false;
+        }
         psFree(maskHeader);
-        psImageCovarianceTransfer(readout->variance, readout->covariance);
-        ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config);
+        if (!ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config)) {
+            psError(PPSTACK_ERR_IO, false, "Unable to write convolved variance %d", i);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(rng);
+            return false;
+        }
 #ifdef TESTING
         {
@@ -129 +182 @@
             psFree(name);
         }
+        {
+            int numCols = readout->image->numCols, numRows = readout->image->numRows;
+            psImage *sn = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            for (int y = 0; y < numRows; y++) {
+                for (int x = 0; x < numCols; x++) {
+                    sn->data.F32[y][x] = readout->image->data.F32[y][x] /
+                        sqrtf(readout->variance->data.F32[y][x]);
+                }
+            }
+            psString name = NULL;
+            psStringAppend(&name, "signoise_%d.fits", i);
+            ppStackWriteImage(name, hdu->header, sn, config);
+            psFree(name);
+            psFree(sn);
+        }
 #endif
 
@@ -136 +204 @@
         psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent);
 
+        // Correct ZP
+        if (options->matchZPs) {
+            // I think I need to take off the exposure time because we're going to set the new exposure time
+            psMetadataItem *zpItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.ZP");
+            zpItem->data.F32 += options->norm->data.F32[i] + 2.5*log10(options->sumExposure);
+
+            psMetadataItem *expItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.EXPOSURE");
+            expItem->data.F32 = options->sumExposure;
+
+            expItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE");
+            expItem->data.F32 = options->sumExposure;
+        }
+
         options->cells->data[i] = psMemIncrRefCounter(inCell);
         if (!ppStackFilesIterateUp(config)) {
@@ -149 +230 @@
     psFree(rng);
 
-    psFree(options->norm); options->norm = NULL;
     psFree(options->sourceLists); options->sourceLists = NULL;
     psFree(options->psf); options->psf = NULL;
@@ -155 +235 @@
 
     if (numGood == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
+        psError(PPSTACK_ERR_REJECTED, false, "No good images to combine.");
         psFree(fpaList);
         psFree(cellList);
@@ -167 +247 @@
         pmFPAview view;                 // View for output
         view.chip = view.cell = view.readout = 0;
+
         pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell
         pmFPA *outFPA = outCell->parent->parent; // Output FPA
+
+        pmCell *unconvCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.UNCONV"); // Unconvolved cell
+        pmFPA *unconvFPA = unconvCell->parent->parent;                                  // Unconvolved FPA
+
         pmConceptsAverageFPAs(outFPA, fpaList);
         pmConceptsAverageCells(outCell, cellList, NULL, NULL, true);
+
+        pmConceptsAverageFPAs(unconvFPA, fpaList);
+        pmConceptsAverageCells(unconvCell, cellList, NULL, NULL, true);
+
         psFree(fpaList);
         psFree(cellList);
@@ -177 +266 @@
         UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure);
         UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(outFPA,  "FPA.ZP",        options->zp);
+
+        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",  options->sumExposure);
+        UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure);
+        UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",        options->zp);
     }
 
@@ -191 +286 @@
         assert(values->n == numGood);
         if (!psVectorSortInPlace(values)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to sort vector.");
+            psError(PPSTACK_ERR_PROG, false, "Unable to sort vector.");
             psFree(values);
             return false;
@@ -211 +306 @@
             numGood = 0;                    // Number of good images
             for (int i = 0; i < num; i++) {
-              if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+                if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
                     continue;
                 }
@@ -229 +324 @@
 
     if (numGood == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
+        psError(PPSTACK_ERR_REJECTED, false, "No good images to combine.");
         return false;
     }
+
+    // Correct chi^2 for renormalisation
+    psBinaryOp(options->matchChi2, options->matchChi2, "/", renorms);
+    for (int i = 0; i < num; i++) {
+        psLogMsg("ppStack", PS_LOG_INFO, "Additional variance for image %d: %f\n",
+                 i, options->matchChi2->data.F32[i]);
+    }
+    psFree(renorms);
 
     return true;

branches/tap_branches/ppStack/src/ppStackFiles.c

@@ -17 +17 @@
 
 /// Files required for the convolution
-static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE",
-                                 "PPSTACK.CONV.KERNEL", NULL };
-
-/// Output files for the combination
-static char *filesCombine[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
-                                "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
+static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
+
+//                                 "PPSTACK.CONV.KERNEL", NULL };
+
+/// Regular (convolved) stack files
+static char *filesStack[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
+                              "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT",
+                              "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2",
+                              NULL };
+/// Unconvolved stack files
+static char *filesUnconv[] = { "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE",
+                               "PPSTACK.UNCONV.EXP", "PPSTACK.UNCONV.EXPNUM", "PPSTACK.UNCONV.EXPWT",
+                               NULL };
 
 /// Files for photometry
@@ -35 +42 @@
       case PPSTACK_FILES_PREPARE:  return filesPrepare;
       case PPSTACK_FILES_CONVOLVE: return filesConvolve;
-      case PPSTACK_FILES_COMBINE:  return filesCombine;
+      case PPSTACK_FILES_STACK:    return filesStack;
+      case PPSTACK_FILES_UNCONV:   return filesUnconv;
       case PPSTACK_FILES_PHOT:     return filesPhot;
       default:
@@ -113 +121 @@
     pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return NULL;
     }
     view->chip = 0;
     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return NULL;
     }
     view->cell = 0;
     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return NULL;
     }
     view->readout = 0;
     if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return NULL;
     }
@@ -139 +151 @@
     view->chip = view->cell = view->readout = 0;
     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return false;
     }
     view->readout = -1;
     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return false;
     }
     view->cell = -1;
     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return false;
     }
     view->chip = -1;
     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        psError(PPSTACK_ERR_IO, false, "File checks failed.");
         return false;
     }
@@ -170 +186 @@
     psFits *fits = psFitsOpen(resolved, "w");
     if (!fits) {
-        psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved);
+        psError(PPSTACK_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved);
         psFree(resolved);
         return false;
     }
     if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
-        psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved);
+        psError(PPSTACK_ERR_IO, false, "Unable to write FITS image %s.", resolved);
         psFitsClose(fits);
         psFree(resolved);
         return false;
     }
-    psFitsClose(fits);
+    if (!psFitsClose(fits)) {
+        psError(PPSTACK_ERR_IO, false, "Unable to close FITS image %s.", resolved);
+        psFree(resolved);
+        return false;
+    }
     psFree(resolved);
     return true;

branches/tap_branches/ppStack/src/ppStackFinish.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <unistd.h>
 #include <pslib.h>
 #include <psmodules.h>
-#include <ppStats.h>
+#include <psphot.h>
 
 #include "ppStack.h"
@@ -21 +22 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
-    // Statistics on output
-    if (options->stats) {
-        psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
-        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
-        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+    bool mdok;                          // Status of MD lookup
+    bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files?
 
-        pmFPAview *view = pmFPAviewAlloc(0); // View to readout
-        view->chip = view->cell = view->readout = 0;
+    // Delete temporary images
+    if (tempDelete && options->convolve) {
+        for (int i = 0; i < options->num; i++) {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                continue;
+            }
 
-        ppStatsFPA(options->stats, options->outRO->parent->parent->parent, view, maskBad, config);
-
-        psFree(view);
-    }
-
-    ppStackMemDump("stats");
-
-    psFree(options->outRO); options->outRO = NULL;
-
-    // Write out the output files
-    if (!ppStackFilesIterateUp(config)) {
-        return false;
+            psString imageResolved = pmConfigConvertFilename(options->convImages->data[i],
+                                                             config, false, false);
+            psString maskResolved = pmConfigConvertFilename(options->convMasks->data[i],
+                                                            config, false, false);
+            psString varianceResolved = pmConfigConvertFilename(options->convVariances->data[i],
+                                                                config, false, false);
+            if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
+                unlink(varianceResolved) == -1) {
+                psWarning("Unable to delete temporary files for image %d", i);
+            }
+            psFree(imageResolved);
+            psFree(maskResolved);
+            psFree(varianceResolved);
+        }
     }
 
 
-    // Write out summary statistics
-    if (options->stats) {
-        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_STACK", 0,
-                         "Total time", psTimerClear("PPSTACK_TOTAL"));
+    return true;
+}
 
-        const char *statsMDC = psMetadataConfigFormat(options->stats);
-        if (!statsMDC || strlen(statsMDC) == 0) {
-            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
-        } else {
-            fprintf(options->statsFile, "%s", statsMDC);
-        }
-        psFree((void *)statsMDC);
-        fclose(options->statsFile); options->statsFile = NULL;
-        pmConfigRunFilenameAddWrite(config, "STATS", psMetadataLookupStr(NULL, config->arguments, "STATS"));
+
+psExit ppStackExitCode(psExit exitValue)
+{
+    if (exitValue != PS_EXIT_SUCCESS) {
+        return exitValue;
     }
 
-
-    // Dump configuration
-    bool mdok;                          // Status of MD lookup
-    psString dump = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); // File for config
-    if (dump) {
-        if (!pmConfigDump(config, dump)) {
-            psError(PS_ERR_IO, false, "Unable to dump configuration.");
-            return false;
+    psErrorCode errorCode = psErrorCodeLast(); // Error code
+    if (errorCode != PS_ERR_NONE) {
+        psErrorStackPrint(stderr, "Unable to perform stack.");
+        switch (errorCode) {
+          case PPSTACK_ERR_UNKNOWN:
+          case PS_ERR_UNKNOWN:
+            psLogMsg("ppStack", PS_LOG_WARN, "Unknown error code: %x", errorCode);
+            exitValue = PS_EXIT_UNKNOWN_ERROR;
+            break;
+          case PS_ERR_IO:
+          case PS_ERR_DB_CLIENT:
+          case PS_ERR_DB_SERVER:
+          case PS_ERR_BAD_FITS:
+          case PS_ERR_OS_CALL_FAILED:
+          case PM_ERR_SYS:
+          case PPSTACK_ERR_IO:
+            psLogMsg("ppStack", PS_LOG_WARN, "I/O error code: %x", errorCode);
+            exitValue = PS_EXIT_SYS_ERROR;
+            break;
+          case PS_ERR_BAD_PARAMETER_VALUE:
+          case PS_ERR_BAD_PARAMETER_TYPE:
+          case PS_ERR_BAD_PARAMETER_NULL:
+          case PS_ERR_BAD_PARAMETER_SIZE:
+          case PPSTACK_ERR_ARGUMENTS:
+          case PPSTACK_ERR_CONFIG:
+            psLogMsg("ppStack", PS_LOG_WARN, "Configuration error code: %x", errorCode);
+            exitValue = PS_EXIT_CONFIG_ERROR;
+            break;
+          case PPSTACK_ERR_PSF:
+          case PSPHOT_ERR_PSF:
+          case PM_ERR_STAMPS:
+          case PM_ERR_SMALL_AREA:
+          case PPSTACK_ERR_REJECTED:
+          case PPSTACK_ERR_DATA:
+            psLogMsg("ppStack", PS_LOG_WARN, "Data error code: %x", errorCode);
+            exitValue = PS_EXIT_DATA_ERROR;
+            break;
+          case PS_ERR_UNEXPECTED_NULL:
+          case PS_ERR_PROGRAMMING:
+          case PPSTACK_ERR_NOT_IMPLEMENTED:
+          case PPSTACK_ERR_PROG:
+            psLogMsg("ppStack", PS_LOG_WARN, "Programming error code: %x", errorCode);
+            exitValue = PS_EXIT_PROG_ERROR;
+            break;
+          default:
+            // It's a programming error if we're not dealing with the error correctly
+            psLogMsg("ppStack", PS_LOG_WARN, "Unrecognised error code: %x", errorCode);
+            exitValue = PS_EXIT_PROG_ERROR;
+            break;
         }
     }
-
-    return true;
+    return exitValue;
 }

branches/tap_branches/ppStack/src/ppStackLoop.c

@@ -10 +10 @@
 #include "ppStackLoop.h"
 
-bool ppStackLoop(pmConfig *config)
+/// Print a summary of the inputs, and return the number of good inputs
+static int stackSummary(const ppStackOptions *options, // Stack options, with input mask
+                        const char *place              // Place in code
+    )
+{
+    int numGood = 0;                // Number of good inputs
+    psString summary = NULL;        // Summary of images
+    for (int i = 0; i < options->num; i++) {
+        psString reason = NULL;         // Reason for rejecting
+        if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) {
+            psStringAppend(&reason, " Good.");
+            numGood++;
+        } else {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CAL) {
+                psStringAppend(&reason, " Calibration failed.");
+            }
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_PSF) {
+                psStringAppend(&reason, " PSF measurement failed.");
+            }
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_MATCH) {
+                psStringAppend(&reason, " PSF matching failed.");
+            }
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CHI2) {
+                psStringAppend(&reason, " PSF matching chi^2 deviant.");
+            }
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_REJECT) {
+                psStringAppend(&reason, " Rejection exceeded threshold.");
+            }
+        }
+        psStringAppend(&summary, "Image %d: %s\n", i, reason);
+        psFree(reason);
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Summary of images for %s:\n%s", place, summary);
+    psFree(summary);
+
+    return numGood;
+}
+
+
+
+bool ppStackLoop(pmConfig *config, ppStackOptions *options)
 {
     assert(config);
@@ -16 +56 @@
     psTimerStart("PPSTACK_TOTAL");
     psTimerStart("PPSTACK_STEPS");
-
-    ppStackOptions *options = ppStackOptionsAlloc(); // Options for stacking
 
     // Setup
     psTrace("ppStack", 1, "Setup....\n");
     if (!ppStackSetup(options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to setup.");
-        psFree(options);
+        psError(psErrorCodeLast(), false, "Unable to setup.");
         return false;
     }
@@ -33 +70 @@
     psTrace("ppStack", 1, "Preparation for stacking: merging sources, determining target PSF....\n");
     if (!ppStackPrepare(options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to prepare for stacking.");
-        psFree(options);
+        psError(psErrorCodeLast(), false, "Unable to prepare for stacking.");
         return false;
     }
@@ -40 +76 @@
              psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("prepare");
-
+    if (options->quality) {
+        // Can't do anything else
+        return true;
+    }
 
     // Convolve inputs
     psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
     if (!ppStackConvolve(options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
-        psFree(options);
+        psError(psErrorCodeLast(), false, "Unable to convolve images.");
         return false;
     }
@@ -53 +91 @@
     ppStackMemDump("convolve");
 
+    // Ensure sufficient inputs
+    {
+        int numGood = stackSummary(options, "initial combination");
+        psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+        bool safe = psMetadataLookupBool(NULL, recipe, "SAFE"); // Be safe when combining
+        if (safe && numGood <= 1) {
+            psError(PPSTACK_ERR_REJECTED, true, "Insufficient inputs for combination with safety on");
+            return false;
+        }
+    }
 
     // Start threading
     ppStackThreadInit();
-    ppStackThreadData *stack = ppStackThreadDataSetup(options, config);
+    ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true);
     if (!stack) {
-        psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
-        psFree(options);
-        return false;
-    }
-    psFree(options->cells); options->cells = NULL;
+        psError(psErrorCodeLast(), false, "Unable to initialise stack threads.");
+        return false;
+    }
 
     // Prepare for combination
-    if (!ppStackCombinePrepare(stack, options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to prepare for combination.");
-        psFree(stack);
-        psFree(options);
+    if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK,
+                               stack, options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
+        psFree(stack);
         return false;
     }
@@ -75 +121 @@
     psTrace("ppStack", 1, "Initial stack of convolved images....\n");
     if (!ppStackCombineInitial(stack, options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to perform initial combination.");
-        psFree(stack);
-        psFree(options);
-        return false;
-    }
-    ppStackMemDump("convolve");
+        psError(psErrorCodeLast(), false, "Unable to perform initial combination.");
+        psFree(stack);
+        return false;
+    }
+    ppStackMemDump("initial");
     psLogMsg("ppStack", PS_LOG_INFO, "Stage 3: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
 
+    // Done with stack inputs for now
+    for (int i = 0; i < options->num; i++) {
+        pmCellFreeData(options->cells->data[i]);
+    }
+    psFree(stack);
 
     // Pixel rejection
     psTrace("ppStack", 1, "Reject pixels....\n");
     if (!ppStackReject(options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to reject pixels.");
-        psFree(stack);
-        psFree(options);
+        psError(psErrorCodeLast(), false, "Unable to reject pixels.");
+        psFree(stack);
         return false;
     }
     psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS"));
-    ppStackMemDump("reject");
-
+
+    // Check inputs
+    {
+        int numGood = stackSummary(options, "final combination");
+        if (numGood <= 0) {
+            psError(PPSTACK_ERR_REJECTED, true, "Insufficient inputs for combination");
+            return false;
+        }
+    }
+
+    stack = ppStackThreadDataSetup(options, config, true);
+    if (!stack) {
+        psError(psErrorCodeLast(), false, "Unable to initialise stack threads.");
+        return false;
+    }
 
     // Final combination
     psTrace("ppStack", 2, "Final stack of convolved images....\n");
-    if (!ppStackCombineFinal(stack, options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination.");
-        psFree(stack);
-        psFree(options);
+    if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true)) {
+        psError(psErrorCodeLast(), false, "Unable to perform final combination.");
+        psFree(stack);
         return false;
     }
@@ -107 +168 @@
     ppStackMemDump("final");
 
+    // Photometry
+    psTrace("ppStack", 1, "Photometering stacked image....\n");
+    if (!ppStackPhotometry(options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to perform photometry.");
+        return false;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("photometry");
 
     // Clean up
     psTrace("ppStack", 2, "Cleaning up after combination....\n");
     if (!ppStackCleanup(stack, options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to clean up.");
-        psFree(stack);
-        psFree(options);
-        return false;
-    }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
+        psError(psErrorCodeLast(), false, "Unable to clean up.");
+        psFree(stack);
+        return false;
+    }
+    for (int i = 0; i < options->num; i++) {
+        pmCellFreeData(options->cells->data[i]);
+    }
+    psFree(stack);
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("cleanup");
 
-    psFree(stack);
-
-
-    // Photometry
-    psTrace("ppStack", 1, "Photometering stacked image....\n");
-    if (!ppStackPhotometry(options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry.");
-        psFree(options);
-        return false;
-    }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
-    ppStackMemDump("photometry");
-
+#if 1
+    // Unconvolved stack --- it's cheap to calculate, compared to everything else!
+    if (options->convolve) {
+        // Start threading
+        ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false);
+        if (!stack) {
+            psError(psErrorCodeLast(), false, "Unable to initialise stack threads.");
+            return false;
+        }
+
+        // Prepare for combination
+        if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV,
+                                   stack, options, config)) {
+            psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
+            psFree(stack);
+            return false;
+        }
+
+        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
+        if (!ppStackCombineFinal(stack, options->origCovars, options, config,
+                                 false, true, false)) {
+            psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination.");
+            psFree(stack);
+            return false;
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
+        ppStackMemDump("unconv");
+
+        if (!ppStackFilesIterateUp(config)) {
+            psError(psErrorCodeLast(), false, "Unable to close files.");
+            psFree(stack);
+            return false;
+        }
+        ppStackFileActivation(config, PPSTACK_FILES_UNCONV, false);
+        options->outRO->data_exists = false;
+        options->outRO->parent->data_exists = false;
+        options->outRO->parent->parent->data_exists = false;
+        psFree(options->outRO);
+        options->outRO = NULL;
+        options->expRO->data_exists = false;
+        options->expRO->parent->data_exists = false;
+        options->expRO->parent->parent->data_exists = false;
+        psFree(options->expRO);
+        options->expRO = NULL;
+
+        for (int i = 0; i < options->num; i++) {
+            pmCellFreeData(options->cells->data[i]);
+        }
+        psFree(stack);
+    }
+    psFree(options->cells); options->cells = NULL;
+#endif
 
     // Finish up
     psTrace("ppStack", 1, "Finishing up....\n");
     if (!ppStackFinish(options, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to finish up.");
-        psFree(options);
-        return false;
-    }
-    psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Final output: %f sec", psTimerClear("PPSTACK_STEPS"));
+        psError(psErrorCodeLast(), false, "Unable to finish up.");
+        return false;
+    }
     ppStackMemDump("finish");
 
-    psFree(options);
     return true;
 }

branches/tap_branches/ppStack/src/ppStackLoop.h

@@ -7 +7 @@
 #include "ppStackOptions.h"
 #include "ppStackThread.h"
+#include "ppStack.h"
 
 
 // Loop over the inputs, doing the combination
 bool ppStackLoop(
-    pmConfig *config                    // Configuration
+    pmConfig *config,                    // Configuration
+    ppStackOptions *options             // Options for stacking
     );
 
@@ -36 +38 @@
 // Prepare for combination
 bool ppStackCombinePrepare(
+    const char *outName,                // Name of output file
+    const char *expName,                // Name of exposure file
+    ppStackFileList files,              // Files of interest
     ppStackThreadData *stack,           // Stack
     ppStackOptions *options,            // Options
@@ -57 +62 @@
 bool ppStackCombineFinal(
     ppStackThreadData *stack,           // Stack
+    psArray *covariances,               // Covariances
     ppStackOptions *options,            // Options
-    pmConfig *config                    // Configuration
+    pmConfig *config,                   // Configuration
+    bool safe,                          // Allow safe combination?
+    bool norm,                          // Normalise images?
+    bool grow                           // Grow rejection masks?
     );
 

branches/tap_branches/ppStack/src/ppStackMatch.c

@@ -18 +18 @@
 #define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
 
-//#define TESTING                         // Enable debugging output
+// #define TESTING                         // Enable debugging output
 
 #ifdef TESTING
@@ -31 +31 @@
     psFree(resolved);
     if (!fits) {
-        psError(PS_ERR_IO, false, "Unable to open previously produced image: %s", name);
+        psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name);
         return false;
     }
     psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
     if (!image) {
-        psError(PS_ERR_IO, false, "Unable to read previously produced image: %s", name);
+        psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name);
         psFitsClose(fits);
         return false;
@@ -115 +115 @@
     psFree(coords);
     psFree(tree);
+    psFree(x);
+    psFree(y);
 
     psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
@@ -144 +146 @@
     psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
 
-    psImage *binned = psphotBackgroundModel(ro, config); // Binned background model
+    psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model
     psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis,
                                                   "PSPHOT.BACKGROUND.BINNING"); // Binning for model
@@ -150 +152 @@
     psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
     if (!psImageUnbin(unbinned, binned, binning)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to unbin background model");
+        psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model");
+        psFree(binned);
         psFree(unbinned);
         return NULL;
     }
-
-    // XXX should these really be here?? (probably not...)
-    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
-    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
-
-     return unbinned;
+    psFree(binned);
+
+    return unbinned;
 }
 
@@ -167 +167 @@
     )
 {
+#if 1
     bool mdok; // Status of metadata lookups
 
@@ -176 +177 @@
     int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
     if (!mdok) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.NUM is not set in the recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe");
         return false;
     }
     float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
     if (!mdok) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MIN is not set in the recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe");
         return false;
     }
     float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
     if (!mdok) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MAX is not set in the recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe");
         return false;
     }
@@ -192 +193 @@
     psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
 
+    psImageCovarianceTransfer(readout->variance, readout->covariance);
     return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
+#else
+    return true;
+#endif
 }
 
@@ -212 +217 @@
     int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
@@ -224 +229 @@
 
     if (!pmReadoutMaskNonfinite(readout, maskVal)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in readout.");
+        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout.");
         return false;
     }
@@ -251 +256 @@
             psFree(resolved);
             if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
+                psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel");
                 psFitsClose(fits);
                 return false;
@@ -257 +262 @@
             psFitsClose(fits);
 
-            if (!readImage(&readout->image, options->imageNames->data[index], config) ||
-                !readImage(&readout->mask, options->maskNames->data[index], config) ||
-                !readImage(&readout->variance, options->varianceNames->data[index], config)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced image.");
+            if (!readImage(&readout->image, options->convImages->data[index], config) ||
+                !readImage(&readout->mask, options->convMasks->data[index], config) ||
+                !readImage(&readout->variance, options->convVariances->data[index], config)) {
+                psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image.");
                 return false;
             }
@@ -269 +274 @@
                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);
 
-            pmSubtractionAnalysis(readout->analysis, NULL, kernels, region,
+            pmSubtractionAnalysis(conv->analysis, NULL, kernels, region,
                                   readout->image->numCols, readout->image->numRows);
 
             psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+            bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
             psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
+            psImageCovarianceSetThreads(oldThreads);
             psFree(readout->covariance);
             readout->covariance = covar;
@@ -292 +299 @@
             int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations
             float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
-            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel
+            float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel
+            float normFrac = psMetadataLookupF32(NULL, ppsub, "NORM.FRAC"); // Fraction of window for normalisn windw
+            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images
+            float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky
+            float covarFrac = psMetadataLookupF32(NULL, ppsub, "COVAR.FRAC"); // Fraction for covariance calculation
+
             const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
             pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
@@ -309 +321 @@
             float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
 
+            bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE");        // Scale kernel parameters?
+            float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling
+            float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling
+            float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling
+            if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
+                psError(PPSTACK_ERR_CONFIG, false,
+                        "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
+                        scaleRef, scaleMin, scaleMax);
+                return false;
+            }
+
+
             // These values are specified specifically for stacking
             const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
@@ -322 +346 @@
             psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
             if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-                psError(PS_ERR_UNKNOWN, false, "Can't measure background for image.");
+                psError(PPSTACK_ERR_DATA, false, "Can't measure background for image.");
                 psFree(fake);
                 psFree(optWidths);
@@ -338 +362 @@
                                                        footprint); // Filtered list of sources
 
+            bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
                                           stampSources, SOURCE_MASK, NULL, NULL, options->psf,
                                           minFlux, footprint + size, false, true)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
+                psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
                 psFree(optWidths);
@@ -347 +372 @@
                 return false;
             }
+            pmReadoutFakeThreads(oldThreads);
 
             fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
 
+#if 1
             // Add the background into the target image
             psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
             psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
+#endif
 
 #ifdef TESTING
@@ -379 +407 @@
 
             if (threads > 0) {
-                pmSubtractionThreadsInit(readout, fake);
+                pmSubtractionThreadsInit();
             }
 
@@ -387 +415 @@
             if (kernel) {
                 if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
-                                               stride, sysError, maskVal, maskBad, maskPoor,
+                                               stride, kernelError, covarFrac, maskVal, maskBad, maskPoor,
                                                poorFrac, badFrac)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
+                    psError(psErrorCodeLast(), false, "Unable to convolve images.");
                     psFree(fake);
                     psFree(optWidths);
@@ -395 +423 @@
                     psFree(conv);
                     if (threads > 0) {
-                        pmSubtractionThreadsFinalize(readout, fake);
+                        pmSubtractionThreadsFinalize();
                     }
                     return false;
                 }
             } else {
-                if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
-                                        threshold, stampSources, stampsName, type, size, order, widths,
-                                        orders, inner, ringsOrder, binning, penalty,
-                                        optimum, optWidths, optOrder, optThresh, iter, rej, sysError,
-                                        maskVal, maskBad, maskPoor, poorFrac, badFrac,
-                                        PM_SUBTRACTION_MODE_2)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+                // Scale the input parameters
+                psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
+                if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy,
+                                                       options->inputSeeing->data.F32[index],
+                                                       options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
+                    psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
                     psFree(fake);
                     psFree(optWidths);
                     psFree(stampSources);
                     psFree(conv);
+                    psFree(widthsCopy);
                     if (threads > 0) {
-                        pmSubtractionThreadsFinalize(readout, fake);
+                        pmSubtractionThreadsFinalize();
                     }
                     return false;
                 }
-            }
+
+                if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
+                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
+                                        orders, inner, ringsOrder, binning, penalty,
+                                        optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
+                                        sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor,
+                                        poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
+                    psError(psErrorCodeLast(), false, "Unable to match images.");
+                    psFree(fake);
+                    psFree(optWidths);
+                    psFree(stampSources);
+                    psFree(conv);
+                    psFree(widthsCopy);
+                    if (threads > 0) {
+                        pmSubtractionThreadsFinalize();
+                    }
+                    return false;
+                }
+                psFree(widthsCopy);
+            }
+
 
 #ifdef TESTING
@@ -447 +495 @@
 
             if (threads > 0) {
-                pmSubtractionThreadsFinalize(readout, fake);
+                pmSubtractionThreadsFinalize();
             }
 
@@ -487 +535 @@
             while ((item = psMetadataGetAndIncrement(iter))) {
                 assert(item->type == PS_DATA_UNKNOWN);
-                // Set the normalisation dimensions, since these will be otherwise unavailable when reading
-                // the images by scans.
                 pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
-                kernel->numCols = readout->image->numCols;
-                kernel->numRows = readout->image->numRows;
-
                 kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
             }
@@ -518 +561 @@
         }
 
+        // Kernel normalisation
+        {
+            double sum = 0.0;           // Sum of chi^2
+            int num = 0;                // Number of measurements of chi^2
+            psString regex = NULL;      // Regular expression
+            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
+            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
+            psFree(regex);
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_TYPE_F32);
+                float norm = item->data.F32; // Normalisation
+                sum += norm;
+                num++;
+            }
+            psFree(iter);
+            float conv = sum/num;       // Mean normalisation from convolution
+            float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
+            float renorm =  stars / conv; // Renormalisation to apply
+            psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n",
+                     index, renorm, conv, stars);
+            psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
+            psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
+        }
+
         // Reject image completely if the maximum deconvolution fraction exceeds the limit
         float deconv = psMetadataLookupF32(NULL, conv->analysis,
                                            PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
         if (deconv > deconvLimit) {
-            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
-                      deconv, deconvLimit);
+            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n",
+                      deconv, deconvLimit, index);
             psFree(conv);
             return NULL;
@@ -542 +610 @@
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-      psWarning("Can't measure background for image.");
-      psErrorClear();
+        psWarning("Can't measure background for image.");
+        psErrorClear();
     } else {
-      if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
-        psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
-                 psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
-        (void)psBinaryOp(readout->image, readout->image, "-",
-                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
-      }
+        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
+            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
+                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
+            (void)psBinaryOp(readout->image, readout->image, "-",
+                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
+        }
     }
 
@@ -560 +628 @@
 
     // Measure the variance level for the weighting
-    if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image.");
-        psFree(rng);
-        psFree(bg);
-        return false;
-    }
-    options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
-                                                  psImageCovarianceFactor(readout->covariance));
-    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
-                     "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]);
+    if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) {
+        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
+            psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
+            psFree(rng);
+            psFree(bg);
+            return false;
+        }
+        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
+                                                      psImageCovarianceFactor(readout->covariance));
+    } else {
+        options->weightings->data.F32[index] = 1.0;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
+             index, options->weightings->data.F32[index]);
 
     psFree(rng);

branches/tap_branches/ppStack/src/ppStackOptions.c

@@ -12 +12 @@
         fclose(options->statsFile);
     }
-    psFree(options->imageNames);
-    psFree(options->maskNames);
-    psFree(options->varianceNames);
+    psFree(options->origImages);
+    psFree(options->origMasks);
+    psFree(options->origVariances);
+    psFree(options->origCovars);
+    psFree(options->convImages);
+    psFree(options->convMasks);
+    psFree(options->convVariances);
     psFree(options->psf);
+    psFree(options->inputSeeing);
+    psFree(options->exposures);
     psFree(options->inputMask);
     psFree(options->sourceLists);
     psFree(options->norm);
+    psFree(options->sources);
     psFree(options->cells);
     psFree(options->kernels);
@@ -24 +31 @@
     psFree(options->matchChi2);
     psFree(options->weightings);
-    psFree(options->covariances);
+    psFree(options->convCovars);
     psFree(options->outRO);
+    psFree(options->expRO);
     psFree(options->inspect);
     psFree(options->rejected);
@@ -37 +45 @@
 
     options->convolve = true;
+    options->matchZPs = true;
+    options->photometry = false;
     options->stats = NULL;
     options->statsFile = NULL;
-    options->imageNames = NULL;
-    options->maskNames = NULL;
-    options->varianceNames = NULL;
+    options->origImages = NULL;
+    options->origMasks = NULL;
+    options->origVariances = NULL;
+    options->origCovars = NULL;
+    options->convImages = NULL;
+    options->convMasks = NULL;
+    options->convVariances = NULL;
     options->num = 0;
+    options->quality = 0;
     options->psf = NULL;
+    options->sumExposure = NAN;
+    options->zp = NAN;
+    options->inputSeeing = NULL;
+    options->exposures = NULL;
+    options->targetSeeing = NAN;
     options->inputMask = NULL;
     options->sourceLists = NULL;
     options->norm = NULL;
+    options->sources = NULL;
     options->cells = NULL;
     options->kernels = NULL;
@@ -54 +75 @@
     options->matchChi2 = NULL;
     options->weightings = NULL;
-    options->covariances = NULL;
+    options->convCovars = NULL;
     options->outRO = NULL;
+    options->expRO = NULL;
     options->inspect = NULL;
     options->rejected = NULL;

branches/tap_branches/ppStack/src/ppStackOptions.h

@@ -9 +9 @@
     // Setup
     bool convolve;                      // Convolve images?
+    bool matchZPs;                      // Adjust relative fluxes based on transparency analysis?
+    bool photometry;                    // Perform photometry?
     psMetadata *stats;                  // Statistics for output
     FILE *statsFile;                    // File to which to write statistics
-    psArray *imageNames, *maskNames, *varianceNames; // Filenames for the temporary convolved images
+    psArray *origImages, *origMasks, *origVariances; // Filenames of the original images
+    psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images
+    psArray *origCovars;                // Original covariances matrices
     int num;                            // Number of inputs
+    int quality;                        // Bad data quality flag
     // Prepare
     pmPSF *psf;                         // Target PSF
+    psVector *inputSeeing;              // Input seeing FWHMs
+    float targetSeeing;                 // Target seeing FWHM
+    psVector *exposures;                // Exposure times
     float sumExposure;                  // Sum of exposure times
+    float zp;                           // Zero point for output
     psVector *inputMask;                // Mask for inputs
     psArray *sourceLists;               // Individual lists of sources for matching
     psVector *norm;                     // Normalisation for each image
+    psArray *sources;                   // Matched sources
     // Convolve
     psArray *cells;                     // Cells for convolved images --- a handle for reading again
@@ -26 +36 @@
     psVector *matchChi2;                // chi^2 for stamps from matching
     psVector *weightings;               // Combination weightings for images (1/noise^2)
-    psArray *covariances;               // Covariance matrices
+    psArray *convCovars;                // Convolved covariance matrices
     // Combine initial
     pmReadout *outRO;                   // Output readout
+    pmReadout *expRO;                   // Exposure readout
     psArray *inspect;                   // Array of arrays of pixels to inspect
     // Rejection

branches/tap_branches/ppStack/src/ppStackPSF.c

@@ -14 +14 @@
                   const psArray *psfs, const psVector *inputMask)
 {
+    bool mdok = false;
+
 #ifndef TESTING
     // Get the recipe values
@@ -24 +26 @@
     int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF
 
+    psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
+    if (!mdok || !maskValStr) {
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find MASK.VAL in recipe");
+        return false;
+    }
+    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
+
     for (int i = 0; i < psfs->n; i++) {
         if (inputMask->data.U8[i]) {
@@ -32 +41 @@
 
     // Solve for the target PSF
-    pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel,
-                               psfOrder, psfOrder);
+    pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal);
     if (!psf) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+        psError(PPSTACK_ERR_PSF, false, "Unable to determine output PSF.");
         return NULL;
     }
@@ -42 +50 @@
     pmPSF *psf = pmPSFBuildSimple("PS_MODEL_PS1_V1", 4.0, 4.0, 0.0, 1.0);
     if (!psf) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to build dummy PSF.");
+        psError(PPSTACK_ERR_PSF, false, "Unable to build dummy PSF.");
         return NULL;
     }

branches/tap_branches/ppStack/src/ppStackPhotometry.c

@@ -19 +19 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
-    bool mdok;                          // Status of MD lookup
-    if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
+    if (!options->photometry) {
         // Nothing to do
         return true;
@@ -27 +26 @@
     psTimerStart("PPSTACK_PHOT");
 
-    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, false);
+    pmFPAfileActivate(config->files, false, NULL);
     ppStackFileActivation(config, PPSTACK_FILES_PHOT, true);
     pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout
-    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true);
+    ppStackFileActivation(config, PPSTACK_FILES_STACK, true);
 
     pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry
@@ -60 +59 @@
                            "Bits to use for mark", markValue);
 
-    pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT");
-    psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
+    psArray *inSources = options->sources;
     if (!inSources) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find input sources");
+        psError(PPSTACK_ERR_PROG, false, "Unable to find input sources");
         psFree(photView);
         return false;
     }
 
+    pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
     if (!psphotReadoutKnownSources(config, photView, inSources)) {
         // This is likely a data quality issue
@@ -73 +72 @@
         psErrorStackPrint(stderr, "Unable to perform photometry on image");
         psWarning("Unable to perform photometry on image --- suspect bad data quality.");
-        if (options->stats && psMetadataLookupS32(NULL, options->stats, "QUALITY") == 0) {
-            psMetadataAddS32(options->stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
-                             "Unable to perform photometry on image", psErrorCodeLast());
+        if (options->quality == 0) {
+            options->quality = psErrorCodeLast();
         }
         psErrorClear();
@@ -84 +82 @@
         !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
         !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
+        psError(PPSTACK_ERR_PROG, false, "Unable to drop PSPHOT internal files.");
         return false;
     }
@@ -92 +90 @@
     if (options->stats) {
         pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources
-        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
-        psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0,
-                         "Number of sources detected", sources ? sources->n : 0);
-        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0,
-                         "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
+        pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // detections
+        if (detections && detections->allSources) {
+            psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", detections->allSources->n);
+        } else {
+            psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", 0);
+        }
+        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
     }
+
     psFree(photView);
 

branches/tap_branches/ppStack/src/ppStackPrepare.c

@@ -26 +26 @@
     float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
     if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.RADIUS in recipe");
         return false;
     }
     float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
     if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.SIGMA in recipe");
         return false;
     }
     float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
     if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
+        psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.FRAC in recipe");
         return false;
     }
@@ -42 +42 @@
     psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
     if (!mdok || !maskValStr) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find MASK.VAL in recipe");
         return false;
     }
@@ -53 +53 @@
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
     if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
+        psError(PPSTACK_ERR_DATA, false, "Unable to measure background for image");
         psFree(stats);
         psFree(rng);
@@ -126 +126 @@
     options->inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
     psVectorInit(options->inputMask, 0);
+    options->exposures = psVectorAlloc(options->num, PS_TYPE_F32);
+    psVectorInit(options->exposures, NAN);
 
     pmFPAfileActivate(config->files, false, NULL);
@@ -134 +136 @@
     }
 
-    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, "^PPSTACK.INPUT$");
-    psMetadataItem *fileItem; // Item from iteration
     psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
     int numCols = 0, numRows = 0;   // Size of image
-    int index = 0;                  // Index for file
-    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
-        assert(fileItem->type == PS_DATA_UNKNOWN);
-        pmFPAfile *inputFile = fileItem->data.V; // An input file
+    options->sumExposure = 0.0;
+    for (int i = 0; i < num; i++) {
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
+        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
+
+        options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE");
+        options->sumExposure += options->exposures->data.F32[i];
 
         // Get list of PSFs, to determine target PSF
         if (options->convolve) {
-            pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
+            pmChip *chip = pmFPAviewThisChip(view, file->fpa); // The chip: holds the PSF
             pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
             if (!psf) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
-                psFree(view);
-                psFree(fileIter);
+                psError(PPSTACK_ERR_PROG, false, "Unable to find PSF.");
+                psFree(view);
                 psFree(psfs);
                 return false;
             }
-            psfs->data[index] = psMemIncrRefCounter(psf);
+            psfs->data[i] = psMemIncrRefCounter(psf);
             psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF");
 
-            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
+            pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
             pmHDU *hdu = pmHDUFromCell(cell);
             assert(hdu && hdu->header);
@@ -163 +165 @@
             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(fileIter);
+                psError(PPSTACK_ERR_PROG, false, "Unable to determine size of image from PSF.");
+                psFree(view);
                 psFree(psfs);
                 return false;
@@ -176 +177 @@
 
 
-        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
-        psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
-        if (!sources) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
-            return NULL;
-        }
-        options->sourceLists->data[index] = psMemIncrRefCounter(sources);
+        bool redoPhot = psMetadataLookupBool(NULL, recipe, "PHOT");
+
+        pmDetections *detections = NULL;
+        if (options->convolve || options->matchZPs || options->photometry || redoPhot) {
+            pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout with sources
+            detections = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.DETECTIONS"); // Sources
+            if (!detections || !detections->allSources) {
+                psWarning("No detections found for image %d --- rejecting.", i);
+                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_CAL;
+                continue;
+            }
+            psAssert (detections->allSources, "missing sources?");
+
+            options->sourceLists->data[i] = psMemIncrRefCounter(detections->allSources);
+        }
 
         // Re-do photometry if we don't trust the source lists
-        if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
-            psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
+        if (redoPhot) {
+            psTrace("ppStack", 2, "Photometering input %d of %d....\n", i, num);
             pmFPAfileActivate(config->files, false, NULL);
-            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
-            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
+            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i);
+            if (options->convolve) {
+                pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
+            }
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File
             pmFPAview *photView = ppStackFilesIterateDown(config);
             if (!photView) {
@@ -198 +210 @@
             pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
 
-            if (!ppStackInputPhotometer(ro, sources, config)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
+            if (!ppStackInputPhotometer(ro, detections->allSources, config)) {
+                psError(psErrorCodeLast(), false, "Unable to do photometry on input sources");
                 psFree(view);
                 psFree(photView);
@@ -213 +225 @@
             ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
         }
-
-        index++;
-    }
-    psFree(fileIter);
+    }
+
+    psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
+    bool havePSFs = false;                                // Do we have any PSFs?
+    options->inputSeeing = psVectorAlloc(num, PS_TYPE_F32);
+    psVectorInit(options->inputSeeing, NAN);
+    for (int i = 0; i < num; i++) {
+        pmPSF *psf = psfs->data[i];     // PSF for image
+        if (!psf) {
+            continue;
+        }
+        havePSFs = true;
+
+        int xNum = PS_MAX(psf->trendNx, 1), yNum = PS_MAX(psf->trendNy, 1); // Number of realisations
+        double sumFWHM = 0.0;           // FWHM for image
+        int numFWHM = 0;                // Number of FWHM measurements
+        for (int y = 0; y < yNum; y++) {
+            float yPos = ((float)y + 0.5) / (float)yNum * numRows;
+            for (int x = 0; x < xNum; x++) {
+                float xPos = ((float)x + 0.5) / (float)xNum * numCols;
+                float fwhm = pmPSFtoFWHM(psf, xPos, yPos); // FWHM for image
+                if (isfinite(fwhm)) {
+                    sumFWHM += fwhm;
+                    numFWHM++;
+                }
+            }
+        }
+        if (numFWHM == 0) {
+            options->inputSeeing->data.F32[i] = NAN;
+            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+            psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", i);
+        } else {
+            options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM;
+        }
+
+        psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
+    }
+    if (havePSFs) {
+        psLogMsg("ppStack", PS_LOG_INFO, "%s", log);
+    }
+    psFree(log);
 
     // Generate target PSF
@@ -223 +272 @@
         psFree(psfs);
         if (!options->psf) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psError(psErrorCodeLast(), false, "Unable to determine output PSF.");
             psFree(view);
             return false;
@@ -229 +278 @@
         psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
                          "Target PSF for stack", options->psf);
-
-        pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.OUTPUT"); // Output chip
+        options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * numCols, 0.5 * numRows); // FWHM for target
+        psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
+
+        pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
         psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
                          "Target PSF", options->psf);
@@ -238 +289 @@
     // Zero point calibration
     if (!ppStackSourcesTransparency(options, view, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
+        psError(PPSTACK_ERR_DATA, false, "Unable to calculate transparency differences");
         psFree(view);
         return false;

branches/tap_branches/ppStack/src/ppStackReadout.c

@@ -23 +23 @@
     psVector *mask = options->inputMask; // Mask for inputs
     psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
+    psVector *exposures = options->exposures;   // Exposure times for each image
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
 
-    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
-                                             weightings, addVariance);
-
-    job->results = inspect;
+    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
+                                         weightings, exposures, addVariance);
     thread->busy = false;
 
-    return inspect ? true : false;
+    return job->results ? true : false;
 }
 
@@ -40 +39 @@
     psArray *args = job->args;          // Arguments
     ppStackThread *thread = args->data[0]; // Thread
-    ppStackOptions *options = args->data[1]; // Options
-    pmConfig *config = args->data[2];   // Configuration
-
-    pmReadout *outRO = options->outRO;  // Output readout
+    psArray *reject = args->data[1];    // Rejected pixels for each image
+    ppStackOptions *options = args->data[2]; // Options
+    pmConfig *config = args->data[3];   // Configuration
+    bool safety = PS_SCALAR_VALUE(args->data[4], U8);    // Safety switch on?
+    bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images?
+
     psVector *mask = options->inputMask; // Mask for inputs
-    psArray *rejected = options->rejected; // Rejected pixels
     psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
+    psVector *exposures = options->exposures;   // Exposure times for each image
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
-
-    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
-                                      weightings, addVariance); // Status of operation
+    psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
+
+    bool status = ppStackReadoutFinal(config, options->outRO, options->expRO, thread->readouts, mask, reject,
+                                      weightings, exposures, addVariance, safety, norm); // Status of operation
 
     thread->busy = false;
 
+    psAssert(status, "Stacking failed.");
+
     return status;
 }
@@ -63 +67 @@
 
     psArray *args = job->args;  // Input arguments
-    psArray *inspect = args->data[0]; // Array of pixel arrays
-    int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest
-
-    psArray *inputs = inspect->data[index]; // Array of interest
-    psPixels *output = NULL;    // Output pixel list
-    for (int i = 0; i < inputs->n; i++) {
-        psPixels *input = inputs->data[i]; // Input pixel list
-        if (!input || input->n == 0) {
-            continue;
-        }
-        output = psPixelsConcatenate(output, input);
-    }
-
-    if (!output) {
-        // If there are no pixels to inspect, then just fake it
-        output = psPixelsAllocEmpty(0);
-    }
-
-    psFree(inputs);
-    inspect->data[index] = output;
+    psArray *inspects = args->data[0]; // Array of pixel arrays
+    psArray *rejects = args->data[1];  // Array of pixel arrays
+    int index = PS_SCALAR_VALUE(args->data[2], S32); // Index of interest
+
+    psArray *inInspects = inspects->data[index]; // Array of interest
+    psArray *inRejects = rejects->data[index]; // Array of interest
+    psAssert(inInspects->n == inRejects->n, "Size should be the same");
+    psPixels *outInspect = NULL, *outReject = NULL; // Output pixel lists
+    for (int i = 0; i < inInspects->n; i++) {
+        psPixels *inInspect = inInspects->data[i]; // Input pixel list
+        if (inInspect && inInspect->n > 0) {
+            outInspect = psPixelsConcatenate(outInspect, inInspect);
+        }
+        psPixels *inReject = inRejects->data[i]; // Input pixel list
+        if (inReject && inReject->n > 0) {
+            outReject = psPixelsConcatenate(outReject, inReject);
+        }
+    }
+
+    // If there are no pixels to inspect, then just fake it
+    if (!outInspect) {
+        outInspect = psPixelsAllocEmpty(0);
+    }
+    if (!outReject) {
+        outReject = psPixelsAllocEmpty(0);
+    }
+
+    psFree(inspects->data[index]);
+    inspects->data[index] = outInspect;
+    psFree(rejects->data[index]);
+    rejects->data[index] = outReject;
 
     return true;
@@ -89 +104 @@
 
 psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
-                               const psVector *mask, const psVector *weightings, const psVector *addVariance)
+                               const psVector *mask, const psVector *weightings, const psVector *exposures,
+                               const psVector *addVariance)
 {
     assert(config);
@@ -104 +120 @@
 
     bool mdok;                          // Status of MD lookup
-    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
+    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
     float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
     float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
@@ -116 +132 @@
     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
@@ -129 +147 @@
         }
 
-#if 0
-        // This doesn't seem to work, so getting the weightings directly from a vector
-        float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
-        if (!mdok || !isfinite(weighting)) {
-            psWarning("No weighting supplied for image %d --- set to unity.", i);
-            weighting = 1.0;
-        } else {
-            psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f", i, weighting);
-        }
-#endif
-
         // Ensure there is a mask, or pmStackCombine will complain
         if (!ro->mask) {
@@ -146 +153 @@
         }
 
-        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
-    }
-
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
-                        combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) {
+        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
+                                          addVariance->data.F32[i]);
+    }
+
+    if (!pmStackCombine(outRO, NULL, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
+                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
         psFree(stack);
@@ -156 +164 @@
     }
 
-    // Save list of pixels to inspect
+    // Save lists of pixels
     psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
+    psArray *reject = psArrayAlloc(num);  // List of pixels rejected
     for (int i = 0; i < num; i++) {
         pmStackData *data = stack->data[i]; // Data for this image
@@ -168 +177 @@
         }
         inspect->data[i] = psMemIncrRefCounter(data->inspect);
+        reject->data[i] = psMemIncrRefCounter(data->reject);
     }
     psFree(stack);
@@ -173 +183 @@
     sectionNum++;
 
-    return inspect;
-}
-
-
-
-bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+    psArray *results = psArrayAlloc(2); // Array of results
+    results->data[0] = inspect;
+    results->data[1] = reject;
+
+    return results;
+}
+
+
+
+bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts,
                          const psVector *mask, const psArray *rejected, const psVector *weightings,
-                         const psVector *addVariance)
+                         const psVector *exposures, const psVector *addVariance, bool safety,
+                         const psVector *norm)
 {
     assert(config);
     assert(outRO);
+    assert(expRO);
     assert(readouts);
     assert(!rejected || readouts->n == rejected->n);
@@ -196 +212 @@
 
     bool mdok;                          // Status of MD lookup
-    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
-    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
-    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
-    float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction
     bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
     bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
@@ -205 +217 @@
     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
@@ -211 +225 @@
     psArray *stack = psArrayAlloc(num); // Array for stacking
 
-    bool entire = true;                 // Combine entire image?
-    if (rejected) {
-        // We have rejection from a previous combination: combine without flagging pixels to inspect
-        entire = false;
-        safe = false;
-        iter = 0;
-        combineRej = NAN;
-        combineSys = NAN;
-    }
+    // We have rejection from a previous combination: combine without flagging pixels to inspect
+    safe &= safety;
+    int iter = 0;
+    float combineRej = NAN;
+    float combineSys = NAN;
+    float combineDiscard = NAN;
 
     for (int i = 0; i < num; i++) {
         pmReadout *ro = readouts->data[i];
-        if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) {
-            // Image completely rejected since previous combination
-            entire = true;
-            continue;
-        } else if (mask->data.U8[i]) {
-            // Image completely rejected before original combination
-            continue;
-        }
-
-#if 0
-        // This doesn't seem to work, so getting the weightings directly from a vector
-        bool mdok;                      // Status of MD lookup
-        float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
-        if (!mdok || !isfinite(weighting)) {
-            psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
-            weighting = 1.0;
-        }
-#endif
+        if (mask->data.U8[i]) {
+            // Image completely rejected
+            continue;
+        }
 
         // Ensure there is a mask, or pmStackCombine will complain
@@ -248 +245 @@
         }
 
-        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i],
+        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
                                              addVariance ? addVariance->data.F32[i] : NAN);
         data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
         stack->data[i] = data;
-    }
-
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0,
-                        iter, combineRej, combineSys, combineDiscard,
-                        entire, useVariance, safe, !rejected)) {
+
+        if (norm) {
+            float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
+            psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(normalise, PS_TYPE_F32));
+            psBinaryOp(ro->variance, ro->variance, "*", psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32));
+        }
+    }
+
+    if (!pmStackCombine(outRO, expRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
+                        combineSys, combineDiscard, useVariance, safe, rejected)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
         psFree(stack);
@@ -264 +266 @@
     pmCell *outCell = outRO->parent;    // Output cell
     pmChip *outChip = outCell->parent;  // Output chip
-
     outRO->data_exists = true;
     outCell->data_exists = true;
     outChip->data_exists = true;
 
+    pmCell *expCell = expRO->parent;    // Exposure cell
+    pmChip *expChip = expCell->parent;  // Exposure chip
+    expRO->data_exists = true;
+    expCell->data_exists = true;
+    expChip->data_exists = true;
+
     psFree(stack);
 

branches/tap_branches/ppStack/src/ppStackReject.c

@@ -10 +10 @@
 #include "ppStackLoop.h"
 
-//#define TESTING
+// #define TESTING
 
 bool ppStackReject(ppStackOptions *options, pmConfig *config)
@@ -23 +23 @@
 
     int num = options->num;             // Number of inputs
-    options->rejected = psArrayAlloc(num);
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
@@ -32 +31 @@
 
     float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution
-    float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor"
     float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject
                                                                      // before rejecting entire image
@@ -53 +51 @@
         psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start
         psArrayAdd(job->args, 1, options->inspect);
+        psArrayAdd(job->args, 1, options->rejected);
         PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
         if (!psThreadJobAddPending(job)) {
@@ -61 +60 @@
     }
     if (!psThreadPoolWait(true)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists.");
+        psError(psErrorCodeLast(), false, "Unable to concatenate inspection lists.");
         return false;
     }
@@ -92 +91 @@
 #endif
 
-        psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
-                                         threshold, poorFrac, stride, options->regions->data[i],
-                                         options->kernels->data[i]); // Rejected pixels
-
 #ifdef TESTING
         {
-            psImage *mask = psPixelsToMask(NULL, reject,
+            psImage *mask = psPixelsToMask(NULL, options->rejected->data[i],
                                            psRegionSet(0, options->numCols - 1, 0, options->numRows - 1),
                                            0xff); // Mask image
             psString name = NULL;           // Name of image
-            psStringAppend(&name, "reject_%03d.fits", i);
+            psStringAppend(&name, "pre_reject_%03d.fits", i);
             pmStackVisualPlotTestImage(mask, name);
             psFits *fits = psFitsOpen(name, "w");
@@ -111 +106 @@
         }
 #endif
+
+        psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
+                                         threshold, stride, options->regions->data[i],
+                                         options->kernels->data[i]); // Rejected pixels
 
         psFree(options->inspect->data[i]);
@@ -127 +126 @@
                           "exceeds limit (%.3f)", i, frac, imageRej);
                 psFree(reject);
-                // reject == NULL means reject image completely
                 reject = NULL;
                 options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_BAD;
                 numRejected++;
+            } else {
+                // Add to list of pixels already rejected
+                reject = psPixelsConcatenate(reject, options->rejected->data[i]);
+                options->rejected->data[i] = psPixelsDuplicates(options->rejected->data[i], reject);
             }
         }
 
-        // Images without a list of rejected pixels (the list may be empty) are rejected completely
-        options->rejected->data[i] = reject;
+#ifdef TESTING
+        {
+            psImage *mask = psPixelsToMask(NULL, options->rejected->data[i],
+                                           psRegionSet(0, options->numCols - 1, 0, options->numRows - 1),
+                                           0xff); // Mask image
+            psString name = NULL;           // Name of image
+            psStringAppend(&name, "reject_%03d.fits", i);
+            pmStackVisualPlotTestImage(mask, name);
+            psFits *fits = psFitsOpen(name, "w");
+            psFree(name);
+            psFitsWriteImage(fits, NULL, mask, 0, NULL);
+            psFree(mask);
+            psFitsClose(fits);
+        }
+#endif
 
         if (options->stats) {
@@ -143 +158 @@
                              "Number of pixels rejected", reject ? reject->n : 0);
         }
+
+        psFree(reject);
         psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
                  psTimerClear("PPSTACK_REJECT"));
@@ -148 +165 @@
 
     psFree(options->inspect); options->inspect = NULL;
-    psFree(options->kernels); options->kernels = NULL;
-    psFree(options->regions); options->regions = NULL;
 
-    if (numRejected >= num - 1) {
-        psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed.");
+    if (numRejected >= num) {
+        psError(PPSTACK_ERR_REJECTED, true, "All inputs completely rejected; unable to proceed.");
         return false;
     }

branches/tap_branches/ppStack/src/ppStackSetup.c

@@ -22 +22 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
+    // XXX : switch to this name? options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis?
+    options->matchZPs = psMetadataLookupBool(NULL, recipe, "ZP"); // Adjust zero points?
+
+    options->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); // Perform photometry?
+
     options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images?
     if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
@@ -37 +42 @@
         options->statsFile = fopen(resolved, "w");
         if (!options->statsFile) {
-            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
+            psError(PPSTACK_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
             psFree(resolved);
             return false;
@@ -43 +48 @@
         psFree(resolved);
         options->stats = psMetadataAlloc();
-        psMetadataAddS32(options->stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
     }
 
     // Generate temporary names for convolved images
-    const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
+    const char *tempDir = psMetadataLookupStr(NULL, config->arguments, "-temp-dir"); // Directory for temps
     if (!tempDir) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
+        tempDir = psMetadataLookupStr(NULL, config->site, "TEMP.DIR");
+    }
+    if (!tempDir) {
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find TEMP.DIR in site configuration");
         return false;
     }
@@ -57 +64 @@
     const char *tempName = basename(outputName);
     if (!tempName) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
+        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to construct basename for temporary files.");
         psFree(outputName);
         return false;
@@ -66 +73 @@
     const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for var maps
     if (!tempImage || !tempMask || !tempVariance) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+        psError(PPSTACK_ERR_CONFIG, false,
                 "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
         psFree(outputName);
@@ -72 +79 @@
     }
 
-    options->imageNames = psArrayAlloc(num);
-    options->maskNames = psArrayAlloc(num);
-    options->varianceNames = psArrayAlloc(num);
+    options->convImages = psArrayAlloc(num);
+    options->convMasks = psArrayAlloc(num);
+    options->convVariances = psArrayAlloc(num);
     for (int i = 0; i < num; i++) {
         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
@@ -81 +88 @@
         psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
         psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
-        options->imageNames->data[i] = imageName;
-        options->maskNames->data[i] = maskName;
-        options->varianceNames->data[i] = varianceName;
+        options->convImages->data[i] = imageName;
+        options->convMasks->data[i] = maskName;
+        options->convVariances->data[i] = varianceName;
     }
     psFree(outputName);
 
+    // Original images
+    options->origImages = psArrayAlloc(num);
+    options->origMasks = psArrayAlloc(num);
+    options->origVariances = psArrayAlloc(num);
+    pmFPAview *view = pmFPAviewAlloc(0);
+    for (int i = 0; i < num; i++) {
+        {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i);
+            options->origImages->data[i] = pmFPAfileName(file, view, config);
+        }
+        {
+            // We want the convolved mask, since that defines the area that has been tested for outliers
+            options->origMasks->data[i] = psMemIncrRefCounter(options->convMasks->data[i]);
+        }
+        {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i);
+            options->origVariances->data[i] = pmFPAfileName(file, view, config);
+        }
+    }
+    psFree(view);
+
     if (!pmConfigMaskSetBits(NULL, NULL, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to determine mask value.");
         return false;
     }

branches/tap_branches/ppStack/src/ppStackSources.c

@@ -8 +8 @@
 
 //#define TESTING                         // Enable debugging output
+
+//#define ASTROMETRY                    // Correct astrometry?
 
 #ifdef TESTING
@@ -61 +63 @@
     PS_ASSERT_PTR_NON_NULL(config, false);
 
+    if (!options->matchZPs && !options->photometry) {
+        options->norm = psVectorAlloc(options->num, PS_TYPE_F32);
+        psVectorInit(options->norm, 0.0);
+        return true;
+    }
+
     psArray *sourceLists = options->sourceLists; // Source lists for each input
     psVector *inputMask = options->inputMask; // Mask for inputs
@@ -79 +87 @@
             }
             source->psfMag += 1.0;
+#ifdef ASTROMETRY
+            if (source->modelPSF) {
+                source->modelPSF->params->data.F32[PM_PAR_XPOS] += 1.0;
+                source->modelPSF->params->data.F32[PM_PAR_YPOS] += 1.0;
+            }
+            if (source->peak) {
+                source->peak->xf += 1.0;
+                source->peak->yf += 1.0;
+            }
+#endif
         }
     }
@@ -105 +123 @@
     psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
     if (!airmassZP) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe.");
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.AIRMASS in recipe.");
         return false;
     }
-
-    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    psMetadata *zpTargetMenu = psMetadataLookupMetadata(NULL, recipe, "ZP.TARGET"); // Target zero point terms
+    if (!zpTargetMenu) {
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.TARGET in recipe.");
+        return false;
+    }
+
+    int num = options->num;             // Number of inputs
     psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
 
-    psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Zero points for each image
+    psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Relative zero points for each image
+    psVector *zpExp = psVectorAlloc(num, PS_TYPE_F32); // Measured zero points for each image (maybe)
+    int zpExpNum = 0;                                  // Number of measured zero points
     const char *filter = NULL;          // Filter name
     float airmassTerm = NAN;            // Airmass term
-    float sumExpTime = 0.0;             // Sum of the exposure time
+    float zpTarget = NAN;               // Target zero point
+    int numGoodImages = 0;              // Number of good images
     for (int i = 0; i < num; i++) {
+        psArray *sources = sourceLists->data[i]; // Source list
+        if (!sources || sources->n == 0) {
+            psLogMsg("ppStack", PS_LOG_WARN, "Image %d has no sources for transparency measurement.", i);
+            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CAL;
+            zp->data.F32[i] = NAN;
+            continue;
+        }
+        numGoodImages++;
+
         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
-        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
 
 #if defined(TESTING) && 0
@@ -135 +169 @@
 #endif
 
-
-        float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
+        float exptime = options->exposures->data.F32[i]; // Exposure time
         float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass
         const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
+        zpExp->data.F32[i] = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.ZP"); // Exposure zero point
+        psLogMsg("ppStack", PS_LOG_INFO,
+                 "Image %d: %.2f sec exposure in %s at airmass %.2f with zero point %.2f",
+                 i, exptime, expFilter, airmass, zpExp->data.F32[i]);
         if (!isfinite(exptime) || exptime == 0 || !isfinite(airmass) || airmass == 0 ||
             !expFilter || strlen(expFilter) == 0) {
-            psError(PS_ERR_UNEXPECTED_NULL, false,
+            psError(PPSTACK_ERR_CONFIG, false,
                     "Unable to find exposure time (%f), airmass (%f) or filter (%s)",
                     exptime, airmass, expFilter);
@@ -147 +184 @@
             return false;
         }
+        if (isfinite(zpExp->data.F32[i])) {
+            zpExp->data.F32[i] += 2.5 * log10(exptime);
+            zpExpNum++;
+        }
 
         if (!filter) {
             filter = expFilter;
-            airmassTerm = psMetadataLookupF32(NULL, airmassZP, filter);
-            if (!isfinite(airmassTerm)) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+            bool mdok;
+            airmassTerm = psMetadataLookupF32(&mdok, airmassZP, filter);
+            if (!mdok || !isfinite(airmassTerm)) {
+                psError(PPSTACK_ERR_CONFIG, false,
                         "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter);
                 psFree(zp);
                 return false;
             }
+            zpTarget = psMetadataLookupF32(&mdok, zpTargetMenu, filter);
+            if (!mdok || !isfinite(zpTarget)) {
+                psError(PPSTACK_ERR_CONFIG, false,
+                        "Unable to find target zero point (ZP.TARGET) for filter %s", filter);
+                psFree(zp);
+                return false;
+            }
         } else if (strcmp(filter, expFilter) != 0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter);
+            psError(PPSTACK_ERR_CONFIG, false, "Filters don't match: %s vs %s", filter, expFilter);
             psFree(zp);
             return false;
         }
 
-        zp->data.F32[i] = airmassTerm * airmass - 2.5 * log10(exptime);
-        sumExpTime += exptime;
-    }
-
-    options->sumExposure = sumExpTime;
+        zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
+    }
+
+    if (numGoodImages == 0) {
+        psLogMsg("ppStack", PS_LOG_WARN, "No images with sources to measure transparency.");
+        options->quality = PPSTACK_ERR_REJECTED;
+        psFree(zp);
+        psFree(zpExp);
+        return true;
+    }
+    if (numGoodImages == 1) {
+        psArray *sources = NULL;        // Sources
+        for (int i = 0; i < num && !sources; i++) {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                continue;
+            }
+            sources = sourceLists->data[i];
+        }
+        options->quality = PPSTACK_ERR_REJECTED;
+        options->sources = psMemIncrRefCounter(sources);
+        options->norm = psVectorAlloc(num, PS_TYPE_F32);
+        psVectorInit(options->norm, 1.0);
+        options->zp = NAN;
+        psLogMsg("ppStack", PS_LOG_WARN, "Single image with sources --- cannot match transparency.");
+        psFree(zp);
+        psFree(zpExp);
+        return true;
+    }
+
+    if (zpExpNum == numGoodImages) {
+        for (int i = 0; i < num; i++) {
+            zp->data.F32[i] = zpExp->data.F32[i];
+        }
+    }
 
     psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
     if (!matches) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
+        psError(PPSTACK_ERR_DATA, false, "Unable to match sources");
         psFree(zp);
         return false;
@@ -180 +258 @@
 #endif
 
-    psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
-                                           iter2, starRej2, starSys2, starLimit,
-                                           transIter, transRej, transThresh); // Transparencies for each image
-    if (!trans) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
-        return false;
-    }
+    if (options->matchZPs) {
+        psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+                                               iter2, starRej2, starSys2, starLimit,
+                                               transIter, transRej, transThresh); // Transparencies per image
+        if (!trans) {
+            psError(PPSTACK_ERR_DATA, false, "Unable to measure transparencies");
+            return false;
+        }
+        for (int i = 0; i < trans->n; i++) {
+            if (!isfinite(trans->data.F32[i])) {
+                inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
+            }
+        }
+        // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
+        // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
+        // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
+        // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
+        // We don't need to know the magnitude zero point for the filter, since it cancels out
+        if (options->matchZPs) {
+            options->norm = psVectorAlloc(num, PS_TYPE_F32);
+            for (int i = 0; i < num; i++) {
+                if (inputMask->data.U8[i] || !isfinite(trans->data.F32[i])) {
+                    continue;
+                }
+                psArray *sources = sourceLists->data[i]; // Sources of interest
+                float magCorr = zp->data.F32[i] + trans->data.F32[i] - 2.5*log10(options->sumExposure);
+                if (zpExpNum == numGoodImages) {
+                    // Using measured zero points, so attempt to set target zero point
+                    magCorr -= zpTarget;
+                }
+                options->norm->data.F32[i] = -magCorr;
+                psLogMsg("ppStack", PS_LOG_INFO,
+                         "Applying scale correction to image %d: %f mag (%f)\n",
+                         i, magCorr, trans->data.F32[i]);
+
+                for (int j = 0; j < sources->n; j++) {
+                    pmSource *source = sources->data[j]; // Source of interest
+                    if (!source) {
+                        continue;
+                    }
+                    source->psfMag -= magCorr;
+                }
+            }
+        }
+
+        if (zpExpNum == numGoodImages) {
+            // Producing image with target zero point
+            options->zp = zpTarget;
+        } else {
+            options->zp = NAN;
+        }
+
 
 #ifdef TESTING
-    dumpMatches("source_mags.dat", num, matches, zp, trans);
-#endif
-
-    for (int i = 0; i < trans->n; i++) {
-        if (!isfinite(trans->data.F32[i])) {
-            inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
-        }
-    }
-
-    // Save best matches SOMEWHERE for future photometry
-    // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack
-    pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
-    psArray *sourcesBest = psArrayAllocEmpty(matches->n);
-
-    // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
-    int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
-    for (int i = 0; i < matches->n; i++) {
-        pmSourceMatch *match = matches->data[i]; // Match of interest
-        if (match->num < minMatches) {
-            continue;
-        }
-
-        // We need to grab a single instance of this source: just take the first available
-        int image = match->image->data.S32[0]; // Index of image
-        int index = match->index->data.S32[0]; // Index of source within image
-        psArray *sources = sourceLists->data[image]; // Sources for image
-        pmSource *source = sources->data[index]; // Source of interest
-
-        psArrayAdd(sourcesBest, sourcesBest->n, source);
-    }
-    psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE,
-                  "psphot sources", sourcesBest);
-    psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
-    psFree(sourcesBest);
+        dumpMatches("source_mags.dat", num, matches, zp, trans);
+#endif
+        psFree(trans);
+
+#ifdef TESTING
+        // Double check: all transparencies should be zero
+        {
+            psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
+            if (!matches) {
+                psError(PPSTACK_ERR_DATA, false, "Unable to match sources");
+                psFree(zp);
+                return false;
+            }
+            psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+                                                   iter2, starRej2, starSys2, starLimit,
+                                                   transIter, transRej, transThresh); // Transparencies
+            for (int i = 0; i < num; i++) {
+                fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
+            }
+            psFree(trans);
+            psFree(matches);
+        }
+#endif
+    }
+
+    psFree(zp);
+    psFree(zpExp);
+
+#ifdef ASTROMETRY
+    // Position offsets
+    {
+        psArray *offsets = pmSourceMatchRelastro(matches, num, tol, iter1, starRej1,
+                                                  iter2, starRej2, starLimit); // Shifts for each image
+        if (!offsets) {
+            psError(PPSTACK_ERR_DATA, false, "Unable to measure offsets");
+            return false;
+        }
+        for (int i = 0; i < num; i++) {
+            if (options->inputMask->data.U8[i]) {
+                continue;
+            }
+            psArray *sources = sourceLists->data[i]; // Sources of interest
+            psVector *offset = offsets->data[i];                      // Offsets for image
+            float dx = offset->data.F32[0], dy = offset->data.F32[1]; // Offsets to apply
+            if (!isfinite(dx) || !isfinite(dy)) {
+                continue;
+            }
+            psLogMsg("ppStack", PS_LOG_INFO, "Applying astrometric correction to image %d: %f,%f\n",
+                     i, dx, dy);
+            for (int j = 0; j < sources->n; j++) {
+                pmSource *source = sources->data[j]; // Source of interest
+                if (!source) {
+                    continue;
+                }
+                if (source->modelPSF) {
+                    source->modelPSF->params->data.F32[PM_PAR_XPOS] -= dx;
+                    source->modelPSF->params->data.F32[PM_PAR_YPOS] -= dy;
+                }
+                if (source->peak) {
+                    source->peak->xf -= dx;
+                    source->peak->yf -= dy;
+                }
+            }
+        }
+        psFree(offsets);
+    }
+#endif
+
+#if (defined TESTING && defined ASTROMETRY)
+        // Double check: all offsets should be zero
+        {
+            psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
+            if (!matches) {
+                psError(PPSTACK_ERR_DATA, false, "Unable to match sources");
+                psFree(zp);
+                return false;
+            }
+            psArray *offsets = pmSourceMatchRelastro(matches, num, tol, iter1, starRej1,
+                                                     iter2, starRej2, starLimit); // Shifts for each image
+            for (int i = 0; i < num; i++) {
+                psVector *offset = offsets->data[i]; // Offsets for image
+                fprintf(stderr, "Offset of image %d: %f,%f\n", i, offset->data.F32[0], offset->data.F32[1]);
+            }
+            psFree(offsets);
+            psFree(matches);
+        }
+#endif
+
+
+    if (options->photometry) {
+        // Save best matches for future photometry
+        psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n);
+
+        // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
+        int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
+        for (int i = 0; i < matches->n; i++) {
+            pmSourceMatch *match = matches->data[i]; // Match of interest
+            if (match->num < minMatches) {
+                continue;
+            }
+
+            // We need to grab a single instance of this source: just take the first available
+            int image = match->image->data.S32[0]; // Index of image
+            int index = match->index->data.S32[0]; // Index of source within image
+            psArray *sources = sourceLists->data[image]; // Sources for image
+            pmSource *source = sources->data[index]; // Source of interest
+
+            psArrayAdd(sourcesBest, sourcesBest->n, source);
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
+    }
 
     psFree(matches);
-
-    // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
-    // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
-    // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
-    // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
-    // We don't need to know the magnitude zero point for the filter, since it cancels out
-    options->norm = psVectorAlloc(num, PS_TYPE_F32);
-    for (int i = 0; i < num; i++) {
-        if (!isfinite(trans->data.F32[i])) {
-            continue;
-        }
-        psArray *sources = sourceLists->data[i]; // Sources of interest
-        float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
-        options->norm->data.F32[i] = magCorr;
-        psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
-
-        for (int j = 0; j < sources->n; j++) {
-            pmSource *source = sources->data[j]; // Source of interest
-            if (!source) {
-                continue;
-            }
-            source->psfMag += magCorr;
-        }
-    }
-    psFree(trans);
-
-#ifdef TESTING
-    // Double check: all transparencies should be zero
-    {
-        psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
-        if (!matches) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
-            psFree(zp);
-            return false;
-        }
-        psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
-                                               transThresh, starRej, starSys);
-        for (int i = 0; i < num; i++) {
-            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
-        }
-        psFree(trans);
-    }
-#endif
 
     return true;

branches/tap_branches/ppStack/src/ppStackThread.c

@@ -56 +56 @@
 }
 
-ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config)
+ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config, bool conv)
 {
     psAssert(options, "Require options");
@@ -62 +62 @@
 
     const psArray *cells = options->cells; // Array of input cells
-    const psArray *imageNames = options->imageNames; // Names of images to read
-    const psArray *maskNames = options->maskNames; // Names of masks to read
-    const psArray *varianceNames = options->varianceNames; // Names of variance maps to read
-    const psArray *covariances = options->covariances; // Covariance matrices (already read)
+    const psArray *imageNames = conv ? options->convImages : options->origImages; // Names of images to read
+    const psArray *maskNames = conv ? options->convMasks : options->origMasks; // Names of masks to read
+    const psArray *varianceNames = conv ? options->convVariances : options->origVariances; // Variance names
+    const psArray *covariances = conv ? options->convCovars : options->origCovars; // Covariance matrices
 
     PS_ASSERT_ARRAY_NON_NULL(cells, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
+    if (imageNames) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
+    }
+    if (maskNames) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
+    }
+    if (varianceNames) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
+    }
+    if (covariances) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
+    }
 
     ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
@@ -87 +95 @@
             continue;
         }
-        // Resolved names
-        psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
-        psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
-        psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
-        stack->imageFits->data[i] = psFitsOpen(imageResolved, "r");
-        stack->maskFits->data[i] = psFitsOpen(maskResolved, "r");
-        stack->varianceFits->data[i] = psFitsOpen(varianceResolved, "r");
-        psFree(imageResolved);
-        psFree(maskResolved);
-        psFree(varianceResolved);
-        if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->varianceFits->data[i]) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s",
-                    (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)varianceNames->data[i]);
-            return NULL;
-        }
+
+// Open an image
+#define IMAGE_OPEN(NAMES, FITS, INDEX)          \
+        if (NAMES) { \
+            psString resolved = pmConfigConvertFilename((NAMES)->data[INDEX], config, false, false); \
+            (FITS)->data[INDEX] = psFitsOpen(resolved, "r");                            \
+            if (!(FITS)->data[INDEX]) { \
+                psError(PPSTACK_ERR_IO, false, "Unable to open file %s", (char*)(NAMES)->data[INDEX]); \
+                psFree(resolved); \
+                return NULL; \
+            } \
+            psFree(resolved); \
+        }
+
+        IMAGE_OPEN(imageNames, stack->imageFits, i);
+        IMAGE_OPEN(maskNames, stack->maskFits, i);
+        IMAGE_OPEN(varianceNames, stack->varianceFits, i);
     }
 
@@ -116 +126 @@
             }
             pmReadout *ro = pmReadoutAlloc(cell); // Readout for thread
-            ro->covariance = psMemIncrRefCounter(covariances->data[j]);
+            if (covariances) {
+                ro->covariance = psMemIncrRefCounter(covariances->data[j]);
+            }
             readouts->data[j] = ro;
         }
@@ -140 +152 @@
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe
     if (!recipe) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE);
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find recipe %s", PPSTACK_RECIPE);
         return NULL;
     }
     int rows = psMetadataLookupS32(NULL, recipe, "ROWS"); // Number of rows to read per chunk
     if (rows <= 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "ROWS is not set in the recipe.");
+        psError(PPSTACK_ERR_CONFIG, false, "ROWS is not set in the recipe.");
         return NULL;
     }
@@ -186 +198 @@
                 psFits *varianceFits = stack->varianceFits->data[i]; // FITS file for variance
 
-
-                int zMax = 0;
+                int zMax = 0;
                 bool keepReading = false;
-                if (pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) {
+
+                if (imageFits && pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunk(ro, imageFits, 0, NULL, rows, overlap, config)) {
-                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT %d",
+                        psError(PPSTACK_ERR_IO, false,
+                                "Unable to read chunk %d for file PPSTACK.INPUT %d",
                                 numChunk, i);
                         *status = false;
@@ -199 +212 @@
                 }
 
-                if (pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) {
+                if (maskFits && pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunkMask(ro, maskFits, 0, NULL, rows, overlap, config)) {
-                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.MASK %d",
+                        psError(PPSTACK_ERR_IO, false,
+                                "Unable to read chunk %d for file PPSTACK.INPUT.MASK %d",
                                 numChunk, i);
                         *status = false;
@@ -209 +223 @@
                 }
 
-                if (pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) {
+                if (varianceFits && pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) {
                     keepReading = true;
                     if (!pmReadoutReadChunkVariance(ro, varianceFits, 0, NULL, rows, overlap, config)) {
-                        psError(PS_ERR_IO, false,
+                        psError(PPSTACK_ERR_IO, false,
                                 "Unable to read chunk %d for file PPSTACK.INPUT.VARIANCE %d",
                                 numChunk, i);
@@ -263 +277 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 2);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 3);
         task->function = &ppStackInspect;
         psThreadTaskAdd(task);
@@ -270 +284 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 3);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6);
         task->function = &ppStackReadoutFinalThread;
         psThreadTaskAdd(task);

branches/tap_branches/ppStack/src/ppStackThread.h

@@ -35 +35 @@
 ppStackThreadData *ppStackThreadDataSetup(
     const ppStackOptions *options,      // Options
-    const pmConfig *config              // Configuration
+    const pmConfig *config,             // Configuration
+    bool conv                           // Use convolved products?
     );
 

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h