Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 27310)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 27400)
@@ -35,6 +35,6 @@
 
 //#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,4 +46,5 @@
     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,4 +58,5 @@
     psFree(buffer->variances);
     psFree(buffer->weights);
+    psFree(buffer->exps);
     psFree(buffer->sources);
     psFree(buffer->limits);
@@ -73,4 +75,5 @@
     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,6 +98,9 @@
 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,4 +127,6 @@
     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,4 +135,8 @@
             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,4 +148,10 @@
     if (var) {
         *var = sumVarianceWeight / PS_SQR(sumWeight);
+    }
+    if (exp) {
+        *exp = sumExp;
+    }
+    if (expWeight) {
+        *expWeight = sumExpWeight / sumWeight;
     }
     return true;
@@ -276,4 +294,5 @@
                            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,4 +311,5 @@
     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,4 +351,5 @@
         }
         pixelWeights->data.F32[numGood] = data->weight;
+        pixelExps->data.F32[numGood] = data->exp;
         pixelSources->data.U16[numGood] = i;
         numGood++;
@@ -347,7 +368,8 @@
     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,4 +384,7 @@
                           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,8 +397,11 @@
     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,4 +421,7 @@
                   varianceValue = pixelVariances->data.F32[0];
               }
+              if (exp) {
+                  expValue = pixelExps->data.F32[0];
+              }
               maskValue = 0;
 #ifdef TESTING
@@ -411,13 +442,11 @@
           // 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,14 +464,12 @@
       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,9 +483,13 @@
         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,5 +834,6 @@
                               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,4 +901,14 @@
     }
 
+    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,5 +979,5 @@
 
 /// 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,4 +988,5 @@
     data->inspect = NULL;
     data->weight = weight;
+    data->exp = exp;
     data->addVariance = addVariance;
 
@@ -952,6 +995,8 @@
 
 /// 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,5 +1004,5 @@
     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,4 +1022,5 @@
     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,7 +1028,9 @@
         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,4 +1093,17 @@
         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,6 +1144,7 @@
             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) {