Changeset 27400

Timestamp:

Mar 22, 2010, 8:34:28 PM (16 years ago)

Author:

Paul Price

Message:

Adding exposure map (both exposure time and number of inputs) to the outputs of ppStack. These are calculated both for the regular (convolved) and unconvolved stacks. These should compress pretty well, but I haven't enabled the compression yet. I also plan to add a weighted exposure time output.

Location:

trunk

Files:

• ippconfig/recipes/filerules-mef.mdc (modified) (1 diff)
• ippconfig/recipes/filerules-simple.mdc (modified) (1 diff)
• ippconfig/recipes/filerules-split.mdc (modified) (1 diff)
• ppStack/src/ppStack.h (modified) (2 diffs)
• ppStack/src/ppStackCamera.c (modified) (3 diffs)
• ppStack/src/ppStackCleanup.c (modified) (1 diff)
• ppStack/src/ppStackCombineFinal.c (modified) (8 diffs)
• ppStack/src/ppStackCombineInitial.c (modified) (1 diff)
• ppStack/src/ppStackCombinePrepare.c (modified) (3 diffs)
• ppStack/src/ppStackFiles.c (modified) (1 diff)
• ppStack/src/ppStackLoop.c (modified) (5 diffs)
• ppStack/src/ppStackLoop.h (modified) (2 diffs)
• ppStack/src/ppStackOptions.c (modified) (4 diffs)
• ppStack/src/ppStackOptions.h (modified) (2 diffs)
• ppStack/src/ppStackPrepare.c (modified) (6 diffs)
• ppStack/src/ppStackReadout.c (modified) (8 diffs)
• ppStack/src/ppStackReject.c (modified) (1 diff)
• ppStack/src/ppStackSources.c (modified) (7 diffs)
• ppStack/src/ppStackThread.c (modified) (1 diff)
• psModules/src/imcombine/pmStack.c (modified) (28 diffs)
• psModules/src/imcombine/pmStack.h (modified) (3 diffs)
• psModules/src/imcombine/pmSubtractionMask.c (modified) (1 diff)
• psModules/src/imcombine/pmSubtractionStamps.c (modified) (1 diff)

trunk/ippconfig/recipes/filerules-mef.mdc

@@ -265 +265 @@
 PPSTACK.OUTPUT.MASK     OUTPUT {OUTPUT}.mk.fits                  MASK      COMP_MASK  FPA        TRUE      NONE
 PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.wt.fits                  VARIANCE  COMP_WT    FPA        TRUE      NONE
+PPSTACK.OUTPUT.EXP      OUTPUT {OUTPUT}.exp.fits                 IMAGE     NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.EXPNUM   OUTPUT {OUTPUT}.num.fits                 MASK      NONE       FPA        TRUE      NONE
 PPSTACK.UNCONV          OUTPUT {OUTPUT}.unconv.fits              IMAGE     COMP_IMG   FPA        TRUE      NONE
 PPSTACK.UNCONV.MASK     OUTPUT {OUTPUT}.unconv.mask.fits         MASK      COMP_MASK  FPA        TRUE      NONE
 PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits           VARIANCE  COMP_WT    FPA        TRUE      NONE
+PPSTACK.UNCONV.EXP      OUTPUT {OUTPUT}.unconv.exp.fits          IMAGE     NONE       FPA        TRUE      NONE
+PPSTACK.UNCONV.EXPNUM   OUTPUT {OUTPUT}.unconv.num.fits          MASK      NONE       FPA        TRUE      NONE
 PPSTACK.TARGET.PSF      OUTPUT {OUTPUT}.target.psf               PSF       NONE       CHIP       TRUE      NONE
 PPSTACK.CONV.KERNEL     OUTPUT {OUTPUT}.{FILE.INDEX}.kernel      SUBKERNEL NONE       FPA        TRUE      NONE

trunk/ippconfig/recipes/filerules-simple.mdc

@@ -214 +214 @@
 PPSTACK.OUTPUT.MASK     OUTPUT {OUTPUT}.mask.fits    MASK      NONE       FPA        TRUE      NONE
 PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.weight.fits  VARIANCE  NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.EXP      OUTPUT {OUTPUT}.exp.fits     IMAGE     NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.EXPNUM   OUTPUT {OUTPUT}.num.fits     MASK      NONE       FPA        TRUE      NONE
 PPSTACK.UNCONV          OUTPUT {OUTPUT}.unconv.fits  IMAGE     COMP_IMG   FPA        TRUE      NONE
 PPSTACK.UNCONV.MASK     OUTPUT {OUTPUT}.unconv.mask.fits MASK  COMP_MASK  FPA        TRUE      NONE
 PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits VARIANCE COMP_WT   FPA        TRUE      NONE
+PPSTACK.UNCONV.EXP      OUTPUT {OUTPUT}.unconv.exp.fits IMAGE  NONE       FPA        TRUE      NONE
+PPSTACK.UNCONV.EXPNUM   OUTPUT {OUTPUT}.unconv.num.fits MASK   NONE       FPA        TRUE      NONE
 PPSTACK.TARGET.PSF      OUTPUT {OUTPUT}.target.psf   PSF       NONE       CHIP       TRUE      NONE
 PPSTACK.CONV.KERNEL     OUTPUT {OUTPUT}.{FILE.INDEX}.kernel SUBKERNEL NONE FPA       TRUE      NONE

trunk/ippconfig/recipes/filerules-split.mdc

@@ -233 +233 @@
 PPSTACK.OUTPUT.MASK     OUTPUT {OUTPUT}.mask.fits                MASK      COMP_MASK  FPA        TRUE      NONE
 PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.wt.fits                  VARIANCE  COMP_WT    FPA        TRUE      NONE
+PPSTACK.OUTPUT.EXP      OUTPUT {OUTPUT}.exp.fits                 IMAGE     NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.EXPNUM   OUTPUT {OUTPUT}.num.fits                 MASK      NONE       FPA        TRUE      NONE
 PPSTACK.UNCONV          OUTPUT {OUTPUT}.unconv.fits              IMAGE     COMP_IMG   FPA        TRUE      NONE
 PPSTACK.UNCONV.MASK     OUTPUT {OUTPUT}.unconv.mask.fits         MASK      COMP_MASK  FPA        TRUE      NONE
 PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits           VARIANCE  COMP_WT    FPA        TRUE      NONE
+PPSTACK.UNCONV.EXP      OUTPUT {OUTPUT}.unconv.exp.fits          IMAGE     NONE       FPA        TRUE      NONE
+PPSTACK.UNCONV.EXPNUM   OUTPUT {OUTPUT}.unconv.num.fits          MASK      NONE       FPA        TRUE      NONE
 PPSTACK.TARGET.PSF      OUTPUT {OUTPUT}.target.psf               PSF       NONE       CHIP       TRUE      NONE
 PPSTACK.CONV.KERNEL     OUTPUT {OUTPUT}.{FILE.INDEX}.kernel      SUBKERNEL NONE       FPA        TRUE      NONE

trunk/ppStack/src/ppStack.h

@@ -69 +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
     );
@@ -83 +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 *exposures,  // Exposure times for each image
                          const psVector *addVariance, // Additional variance for rejection
                          bool safety,                 // Enable safety switch?

trunk/ppStack/src/ppStackCamera.c

@@ -283 +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.MASK"));
+        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;
+
+
+
     if (havePSFs) {
         pmFPA *psfFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain PSF
@@ -302 +340 @@
     }
 
-#if 1
     // Unconvolved stack
     pmFPA *unconvFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain unconvolved output
@@ -351 +388 @@
         unconvVariance->save = true;
     }
-#endif
+
+
+    // 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;
+    }
+
+    // Output mask
+    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;
+
 
     // Output JPEGs

trunk/ppStack/src/ppStackCleanup.c

@@ -86 +86 @@
         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

trunk/ppStack/src/ppStackCombineFinal.c

@@ -14 +14 @@
 //#define TESTING                         // Enable test output
 
-bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances,
-                         ppStackOptions *options, pmConfig *config, bool safe, bool normalise, bool grow)
+bool ppStackCombineFinal(ppStackThreadData *stack, psArray *covariances, ppStackOptions *options,
+                         pmConfig *config, bool safe, bool normalise, bool grow)
 {
     psAssert(stack, "Require stack");
@@ -21 +21 @@
     psAssert(config, "Require configuration");
 
-    int numCols = target->image->numCols, numRows = target->image->numRows; // Size of image
+    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
@@ -43 +45 @@
     }
 
-    if (!target->mask) {
-        target->mask = psImageAlloc(target->image->numCols, target->image->numRows, PS_TYPE_IMAGE_MASK);
+    if (!outRO->mask) {
+        outRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+    }
+    if (!expRO->mask) {
+        expRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     }
 
@@ -62 +67 @@
         }
 
-        // call: ppStackReadoutFinal(config, target, readouts, rejected)
+        // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
         psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
-        psArrayAdd(job->args, 1, target);
         psArrayAdd(job->args, 1, thread);
         psArrayAdd(job->args, 1, reject);
@@ -105 +109 @@
         }
         if (sumWeights > 0.0) {
-            target->covariance = psImageCovarianceSum(covariances);
-            psBinaryOp(target->covariance->image, target->covariance->image, "/",
+            outRO->covariance = psImageCovarianceSum(covariances);
+            psBinaryOp(outRO->covariance->image, outRO->covariance->image, "/",
                        psScalarAlloc(sumWeights, PS_TYPE_F32));
-            psImageCovarianceTransfer(target->variance, target->covariance);
+            psImageCovarianceTransfer(outRO->variance, outRO->covariance);
         }
     } else {
-        target->covariance = psImageCovarianceNone();
+        outRO->covariance = psImageCovarianceNone();
     }
 
@@ -127 +131 @@
             wcsDone = true;
             pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
-            pmHDU *outHDU = pmHDUFromCell(target->parent); // Output HDU
-            pmChip *outChip = target->parent->parent; // Output chip
+            pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
+            pmChip *outChip = outRO->parent->parent; // Output chip
             pmFPA *outFPA = outChip->parent; // Output FPA
             if (!outHDU || !inHDU) {
@@ -151 +155 @@
     }
 
+    // Set exposure time correctly
+    {
+        float exptime = 0.0;            // Summed exposure time
+        for (int i = 0; i < options->num; i++) {
+            if (options->inputMask) {
+                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(target->parent);
+    pmHDU *hdu = pmHDUFromCell(outRO->parent);
     if (!hdu) {
         psError(PPSTACK_ERR_PROG, false, "Unable to find HDU for output.");
@@ -168 +192 @@
     psStringAppend(&name, "combined_image_final_%d.fits", pass);
     pass++;
-    ppStackWriteImage(name, NULL, target->image, config);
+    ppStackWriteImage(name, NULL, outRO->image, config);
     psStringSubstitute(&name, "mask", "image");
-    ppStackWriteImage(name, NULL, target->mask, config);
+    ppStackWriteImage(name, NULL, outRO->mask, config);
     psStringSubstitute(&name, "variance", "mask");
-    ppStackWriteImage(name, NULL, target->variance, config);
+    ppStackWriteImage(name, NULL, outRO->variance, config);
     psFree(name);
 
-    pmStackVisualPlotTestImage(target->image, "combined_image_final.fits");
+    pmStackVisualPlotTestImage(outRO->image, "combined_image_final.fits");
 #endif
 

trunk/ppStack/src/ppStackCombineInitial.c

@@ -89 +89 @@
     }
 
+    ppStackMemDump("initial");
+
 #ifdef TESTING
     ppStackWriteImage("combined_image_initial.fits", NULL, options->outRO->image, config);

trunk/ppStack/src/ppStackCombinePrepare.c

@@ -10 +10 @@
 #include "ppStackLoop.h"
 
-bool ppStackCombinePrepare(pmReadout **ro, const char *name, ppStackFileList files, ppStackThreadData *stack,
+bool ppStackCombinePrepare(const char *outName, const char *expName,
+                           ppStackFileList files, ppStackThreadData *stack,
                            ppStackOptions *options, pmConfig *config)
 {
@@ -32 +33 @@
     }
 
-    pmCell *cell = pmFPAfileThisCell(config->files, view, name); // Output cell
-    *ro = pmReadoutAlloc(cell); // 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);
@@ -42 +46 @@
     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
 
-    if (!pmReadoutStackDefineOutput(*ro, col0, row0, numCols, numRows, true, true, maskBad)) {
+    if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
+        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;

trunk/ppStack/src/ppStackFiles.c

@@ -23 +23 @@
 /// Regular (convolved) stack files
 static char *filesStack[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
-                              "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
+                              "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM",
+                              "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2",
+                              NULL };
 /// Unconvolved stack files
-static char *filesUnconv[] = { "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE", NULL };
+static char *filesUnconv[] = { "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE",
+                               "PPSTACK.UNCONV.EXP", "PPSTACK.UNCONV.EXPNUM", NULL };
 
 /// Files for photometry

trunk/ppStack/src/ppStackLoop.c

@@ -111 +111 @@
 
     // Prepare for combination
-    if (!ppStackCombinePrepare(&options->outRO, "PPSTACK.OUTPUT", PPSTACK_FILES_STACK,
+    if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK,
                                stack, options, config)) {
         psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
@@ -160 +160 @@
     // Final combination
     psTrace("ppStack", 2, "Final stack of convolved images....\n");
-    if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config,
-                             false, false, true)) {
+    if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true)) {
         psError(psErrorCodeLast(), false, "Unable to perform final combination.");
         psFree(stack);
@@ -203 +202 @@
 
         // Prepare for combination
-        if (!ppStackCombinePrepare(&options->unconvRO, "PPSTACK.UNCONV", PPSTACK_FILES_UNCONV,
+        if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV,
                                    stack, options, config)) {
             psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
@@ -211 +210 @@
 
         psTrace("ppStack", 2, "Stack of unconvolved images....\n");
-        if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config,
+        if (!ppStackCombineFinal(stack, options->origCovars, options, config,
                                  false, true, false)) {
             psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination.");
@@ -226 +225 @@
         }
         ppStackFileActivation(config, PPSTACK_FILES_UNCONV, false);
-        options->unconvRO->data_exists = false;
-        options->unconvRO->parent->data_exists = false;
-        options->unconvRO->parent->parent->data_exists = false;
-        psFree(options->unconvRO);
-        options->unconvRO = NULL;
+        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++) {

trunk/ppStack/src/ppStackLoop.h

@@ -38 +38 @@
 // Prepare for combination
 bool ppStackCombinePrepare(
-    pmReadout **readout,                // Readout to set
-    const char *name,                   // Name of file
+    const char *outName,                // Name of output file
+    const char *expName,                // Name of exposure file
     ppStackFileList files,              // Files of interest
     ppStackThreadData *stack,           // Stack
@@ -61 +61 @@
 // Final combination
 bool ppStackCombineFinal(
-    pmReadout *target,                  // Target readout
     ppStackThreadData *stack,           // Stack
     psArray *covariances,               // Covariances

trunk/ppStack/src/ppStackOptions.c

@@ -21 +21 @@
     psFree(options->psf);
     psFree(options->inputSeeing);
+    psFree(options->exposures);
     psFree(options->inputMask);
     psFree(options->sourceLists);
@@ -32 +33 @@
     psFree(options->convCovars);
     psFree(options->outRO);
-    psFree(options->unconvRO);
+    psFree(options->expRO);
     psFree(options->inspect);
     psFree(options->rejected);
@@ -61 +62 @@
     options->zp = NAN;
     options->inputSeeing = NULL;
+    options->exposures = NULL;
     options->targetSeeing = NAN;
     options->inputMask = NULL;
@@ -75 +77 @@
     options->convCovars = NULL;
     options->outRO = NULL;
-    options->unconvRO = NULL;
+    options->expRO = NULL;
     options->inspect = NULL;
     options->rejected = NULL;

trunk/ppStack/src/ppStackOptions.h

@@ -22 +22 @@
     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
@@ -38 +39 @@
     // Combine initial
     pmReadout *outRO;                   // Output readout
-    pmReadout *unconvRO;                // Unconvolved readout
+    pmReadout *expRO;                   // Exposure readout
     psArray *inspect;                   // Array of arrays of pixels to inspect
     // Rejection

trunk/ppStack/src/ppStackPrepare.c

@@ -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(PPSTACK_ERR_PROG, false, "Unable to find PSF.");
                 psFree(view);
-                psFree(fileIter);
                 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);
@@ -165 +167 @@
                 psError(PPSTACK_ERR_PROG, false, "Unable to determine size of image from PSF.");
                 psFree(view);
-                psFree(fileIter);
                 psFree(psfs);
                 return false;
@@ -180 +181 @@
         pmDetections *detections = NULL;
         if (options->convolve || options->matchZPs || options->photometry || redoPhot) {
-            pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
+            pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout with sources
             detections = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.DETECTIONS"); // Sources
             if (!detections) {
@@ -188 +189 @@
             psAssert (detections->allSources, "missing sources?");
 
-            options->sourceLists->data[index] = psMemIncrRefCounter(detections->allSources);
+            options->sourceLists->data[i] = psMemIncrRefCounter(detections->allSources);
         }
 
         // Re-do photometry if we don't trust the source lists
         if (redoPhot) {
-            psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
+            psTrace("ppStack", 2, "Photometering input %d of %d....\n", i, num);
             pmFPAfileActivate(config->files, false, NULL);
-            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index);
+            ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i);
             if (options->convolve) {
                 pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL");
             }
-            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File
             pmFPAview *photView = ppStackFilesIterateDown(config);
             if (!photView) {
@@ -223 +224 @@
             ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true);
         }
-
-        index++;
-    }
-    psFree(fileIter);
+    }
 
     psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message

trunk/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
 
     job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
-                                         weightings, addVariance);
+                                         weightings, exposures, addVariance);
     thread->busy = false;
 
@@ -37 +38 @@
 
     psArray *args = job->args;          // Arguments
-    pmReadout *target = args->data[0];  // Output readout
-    ppStackThread *thread = args->data[1]; // Thread
-    psArray *reject = args->data[2];    // Rejected pixels for each image
-    ppStackOptions *options = args->data[3]; // Options
-    pmConfig *config = args->data[4];   // Configuration
-    bool safety = PS_SCALAR_VALUE(args->data[5], U8);    // Safety switch on?
-    bool normalise = PS_SCALAR_VALUE(args->data[6], U8); // Normalise images?
+    ppStackThread *thread = args->data[0]; // Thread
+    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
     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
     psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
 
-    bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, reject,
-                                      weightings, addVariance, safety, norm); // Status of operation
+    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;
@@ -101 +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);
@@ -149 +153 @@
         }
 
-        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
-    }
-
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
+        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.");
@@ -187 +192 @@
 
 
-bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts,
                          const psVector *mask, const psArray *rejected, const psVector *weightings,
-                         const psVector *addVariance, bool safety, const psVector *norm)
+                         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);
@@ -238 +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;
@@ -250 +257 @@
     }
 
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
+    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.");
@@ -259 +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);
 

trunk/ppStack/src/ppStackReject.c

@@ -162 +162 @@
         psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
                  psTimerClear("PPSTACK_REJECT"));
-
-        ppStackMemDump("reject");
     }
 

trunk/ppStack/src/ppStackSources.c

@@ -64 +64 @@
 
     if (!options->matchZPs && !options->photometry) {
-        int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-        options->norm = psVectorAlloc(num, PS_TYPE_F32);
-        psVectorInit (options->norm, 0.0);
-
-        // XXX do I need to set this?
-        // options->sumExposure = sumExpTime;
-
+        options->norm = psVectorAlloc(options->num, PS_TYPE_F32);
+        psVectorInit(options->norm, 0.0);
         return true;
     }
@@ -137 +132 @@
     }
 
-    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    int num = options->num;             // Number of inputs
     psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
 
@@ -146 +141 @@
     float airmassTerm = NAN;            // Airmass term
     float zpTarget = NAN;               // Target zero point
-    float sumExpTime = 0.0;             // Sum of the exposure time
     int numGoodImages = 0;              // Number of good images
     for (int i = 0; i < num; i++) {
@@ -160 +154 @@
 
         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
-        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
 
 #if defined(TESTING) && 0
@@ -177 +170 @@
 #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
@@ -221 +214 @@
 
         zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
-        sumExpTime += exptime;
-
-    }
-
-    options->sumExposure = sumExpTime;
+    }
 
     if (numGoodImages == 0) {
@@ -291 +280 @@
                 }
                 psArray *sources = sourceLists->data[i]; // Sources of interest
-                float magCorr = zp->data.F32[i] + trans->data.F32[i] - 2.5*log10(sumExpTime);
+                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

trunk/ppStack/src/ppStackThread.c

@@ -284 +284 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7);
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6);
         task->function = &ppStackReadoutFinalThread;
         psThreadTaskAdd(task);

trunk/psModules/src/imcombine/pmStack.c

@@ -35 +35 @@
 
 //#define TESTING                         // Enable test output
-//#define TEST_X 2148-1                     // x coordinate to examine
-//#define TEST_Y 248-1                     // y coordinate to examine
+//#define TEST_X 843-1                     // x coordinate to examine
+//#define TEST_Y 813-1                     // y coordinate to examine
 //#define TEST_RADIUS 0                    // Radius to examine
 
@@ -46 +46 @@
     psVector *variances;                // Pixel variances
     psVector *weights;                  // Pixel weightings
+    psVector *exps;                     // Pixel exposures
     psVector *sources;                  // Pixel sources (which image did they come from?)
     psVector *limits;                   // Rejection limits
@@ -57 +58 @@
     psFree(buffer->variances);
     psFree(buffer->weights);
+    psFree(buffer->exps);
     psFree(buffer->sources);
     psFree(buffer->limits);
@@ -73 +75 @@
     buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
+    buffer->exps = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
     buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
@@ -95 +98 @@
 static bool combinationMeanVariance(float *mean, // Mean value, to return
                                     float *var, // Variance value, to return
+                                    float *exp, // Exposure time, to return
+                                    float *expWeight,          // Weighted exposure time, to return
                                     const psVector *values, // Values to combine
                                     const psVector *variances, // Pixel variances to combine
+                                    const psVector *exps,      // Exposure times to combine
                                     const psVector *weights // Weights to apply
                                     )
@@ -121 +127 @@
     float sumVarianceWeight = 0.0;     // Sum of the pixel variances multiplied by the global weights
     float sumWeight = 0.0;              // Sum of the image weights
+    float sumExp = 0.0;                 // Sum of the exposure time
+    float sumExpWeight = 0.0;           // Sum of the exposure time multiplied by the global weights
     for (int i = 0; i < values->n; i++) {
         sumValueWeight += values->data.F32[i] * weights->data.F32[i];
@@ -127 +135 @@
             sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
         }
+        if (exps) {
+            sumExp += exps->data.F32[i];
+            sumExpWeight += exps->data.F32[i] * weights->data.F32[i];
+        }
     }
 
@@ -136 +148 @@
     if (var) {
         *var = sumVarianceWeight / PS_SQR(sumWeight);
+    }
+    if (exp) {
+        *exp = sumExp;
+    }
+    if (expWeight) {
+        *expWeight = sumExpWeight / sumWeight;
     }
     return true;
@@ -276 +294 @@
                            const psArray *inputs, // Stack data
                            const psVector *weights, // Global (single value) weights for data, or NULL
+                           const psVector *exps,    // Exposures for data, or NULL
                            const psVector *addVariance, // Additional variance for data
                            const psVector *reject, // Indices of pixels to reject, or NULL
@@ -292 +311 @@
     psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
     psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
+    psVector *pixelExps = buffer->exps;       // Exposure times
     psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
     psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
@@ -331 +351 @@
         }
         pixelWeights->data.F32[numGood] = data->weight;
+        pixelExps->data.F32[numGood] = data->exp;
         pixelSources->data.U16[numGood] = i;
         numGood++;
@@ -347 +368 @@
     if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
         for (int i = 0; i < numGood; i++) {
-            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %d\n",
+            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d\n",
                     i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
-                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]);
+                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
+                    pixelSuspects->data.U8[i]);
         }
     }
@@ -362 +384 @@
                           psImage *mask, // Combined mask, for output
                           psImage *variance, // Combined variance map, for output
+                          psImage *exp,   // Exposure map (time), for output
+                          psImage *expnum,       // Exposure map (number) for output
+                          psImage *expweight,    // Exposure map (weighted time) for output
                           int num,      // Number of good pixels
                           combineBuffer *buffer, // Buffer with vectors
@@ -372 +397 @@
     psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
     psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
+    psVector *pixelExps = buffer->exps;       // Exposure times
 
     // Default option is that the pixel is bad
     float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
     psImageMaskType maskValue = bad;    // Value for combined mask
+    float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted)
+
     switch (num) {
       case 0: {
@@ -393 +421 @@
                   varianceValue = pixelVariances->data.F32[0];
               }
+              if (exp) {
+                  expValue = pixelExps->data.F32[0];
+              }
               maskValue = 0;
 #ifdef TESTING
@@ -411 +442 @@
           // Automatically accept the mean of the pixels only if we're not playing safe
           if (!safe) {
-              float mean, var;   // Mean and variance from combination
-              if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
-                  imageValue = mean;
-                  varianceValue = var;
+              if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
+                                          pixelData, pixelVariances, pixelExps, pixelWeights)) {
                   maskValue = 0;
 #ifdef TESTING
                   if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
                       fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n",
-                              x, y, mean, var);
+                              x, y, imageValue, varianceValue);
                   }
 #endif
@@ -435 +464 @@
       default: {
           // Can combine without too much worrying
-          float mean, var;           // Mean and variance of the combination
-          if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
+          if (!combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
+                                       pixelData, pixelVariances, pixelExps, pixelWeights)) {
               break;
           }
-          imageValue = mean;
-          varianceValue = var;
           maskValue = 0;
 #ifdef TESTING
           if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-              fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, mean, var);
+              fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
           }
 #endif
@@ -456 +483 @@
         variance->data.F32[y][x] = varianceValue;
     }
-
-#ifdef TESTING
-    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
-#endif
-
+    if (exp) {
+        exp->data.F32[y][x] = expValue;
+    }
+    if (expnum) {
+        expnum->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
+    }
+    if (expweight) {
+        expweight->data.F32[y][x] = expWeightValue;
+    }
 
     return;
@@ -803 +834 @@
                               int *numCols, int *numRows, // Size of (sky) images
                               const psArray *input, // Input array of pmStackData to validate
-                              const pmReadout *output // Output readout
+                              const pmReadout *output, // Output readout
+                              const pmReadout *exp    // Exposure map
     )
 {
@@ -869 +901 @@
     }
 
+    if (exp) {
+        PM_ASSERT_READOUT_NON_NULL(exp, false);
+        if (exp->image) {
+            PS_ASSERT_IMAGES_SIZE_EQUAL(exp->image, output->image, false);
+        }
+        if (exp->mask) {
+            PS_ASSERT_IMAGES_SIZE_EQUAL(exp->mask, output->image, false);
+        }
+    }
+
     return true;
 }
@@ -937 +979 @@
 
 /// Constructor
-pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance)
+pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float exp, float addVariance)
 {
     pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return
@@ -946 +988 @@
     data->inspect = NULL;
     data->weight = weight;
+    data->exp = exp;
     data->addVariance = addVariance;
 
@@ -952 +995 @@
 
 /// Stack input images
-bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
-                    psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic,
+bool pmStackCombine(pmReadout *combined, pmReadout *expmaps, psArray *input,
+                    psImageMaskType maskVal, psImageMaskType maskSuspect,
+                    psImageMaskType bad, int kernelSize,
+                    float iter, float rej, float sys, float olympic,
                     bool useVariance, bool safe, bool rejection)
 {
@@ -959 +1004 @@
     int num;                            // Number of inputs
     int numCols, numRows;               // Size of (sky) images
-    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
+    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined, expmaps)) {
         return false;
     }
@@ -977 +1022 @@
     psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image
     psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image
+    psVector *exps = psVectorAlloc(num, PS_TYPE_F32);    // Exposure times for each image
     psArray *stack = psArrayAlloc(num); // Stack of readouts
     for (int i = 0; i < num; i++) {
@@ -982 +1028 @@
         if (!data) {
             weights->data.F32[i] = 0.0;
+            exps->data.F32[i] = NAN;
             continue;
         }
         weights->data.F32[i] = data->weight;
+        exps->data.F32[i] = data->exp;
         pmReadout *ro = data->readout;  // Readout of interest
         stack->data[i] = psMemIncrRefCounter(ro);
@@ -1045 +1093 @@
         combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
         combinedVariance = combined->variance;
+    }
+
+    psImage *exp = NULL, *expnum = NULL; // Exposure map and exposure number
+    if (expmaps) {
+        if (!expmaps->image) {
+            expmaps->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        exp = expmaps->image;
+
+        if (!expmaps->mask) {
+            expmaps->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+        }
+        expnum = expmaps->mask;
     }
 
@@ -1083 +1144 @@
             bool suspect;               // Suspect pixels in stack?
             combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
-                           input, weights, addVariance, reject, x, y, maskVal, maskSuspect);
-            combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe);
+                           input, weights, exps, addVariance, reject, x, y, maskVal, maskSuspect);
+            combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, NULL,
+                          num, buffer, x, y, bad, safe);
 
             if (iter > 0) {

trunk/psModules/src/imcombine/pmStack.h

@@ -31 +31 @@
     psPixels *inspect;                  ///< Pixels to inspect
     float weight;                       ///< Relative weighting for image
+    float exp;                          ///< Exposure for image
     float addVariance;                  ///< Additional variance when rejecting
 } pmStackData;
@@ -37 +38 @@
 pmStackData *pmStackDataAlloc(pmReadout *readout, ///< Warped readout (sky cell)
                               float weight, ///< Weight to apply
+                              float exp,    ///< Exposure time for input
                               float addVariance ///< Additional variance when rejecting
     );
@@ -42 +44 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined,///< Combined readout (output)
+                    pmReadout *expmap,  ///< Exposure map (output)
                     psArray *input,     ///< Input array of pmStackData
                     psImageMaskType maskVal, ///< Mask value of bad pixels

trunk/psModules/src/imcombine/pmSubtractionMask.c

@@ -92 +92 @@
                     continue;
                 }
+                if (imageData1[y][x] > 50000) {
+                    maskData1[y][x] = maskVal;
+                    numBad++;
+                    continue;
+                }
+                if (imageData2[y][x] > 50000) {
+                    maskData2[y][x] = maskVal;
+                    numBad++;
+                    continue;
+                }
                 xMin = PS_MIN(xMin, x);
                 xMax = PS_MAX(xMax, x);

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -107 +107 @@
             if ((image1 && image1->data.F32[y][x] < thresh1) ||
                 (image2 && image2->data.F32[y][x] < thresh2)) {
+                continue;
+            }
+            if ((image1 && image1->data.F32[y][x] > 30000) ||
+                (image2 && image2->data.F32[y][x] > 30000)) {
                 continue;
             }