Index: /trunk/psModules/src/detrend/pmShutterCorrection.c
===================================================================
--- /trunk/psModules/src/detrend/pmShutterCorrection.c	(revision 9314)
+++ /trunk/psModules/src/detrend/pmShutterCorrection.c	(revision 9315)
@@ -121,5 +121,5 @@
         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                 "Not enough good values to guess shutter correction.\n");
-        goto ERROR;
+        goto GUESS_ERROR;
     }
     tmpX->data.F32[0] = exptime->data.F32[index];
@@ -133,5 +133,5 @@
         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                 "Exposure times are all identical --- cannot guess shutter correction.\n");
-        goto ERROR;
+        goto GUESS_ERROR;
     }
     tmpY->data.F32[1] = counts->data.F32[index];
@@ -141,5 +141,5 @@
     if (!psVectorFitPolynomial1D(line, NULL, 0, tmpY, NULL, tmpX)) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the time offset.\n");
-        goto ERROR;
+        goto GUESS_ERROR;
     }
     float ratio = psPolynomial1DEval(line, 0.0) / corr->scale;
@@ -170,5 +170,5 @@
     if (!psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX)) {
         psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the reference offset.\n");
-        goto ERROR;
+        goto GUESS_ERROR;
     }
     corr->offref = psPolynomial1DEval(line, value);
@@ -181,5 +181,5 @@
     return corr;
 
-ERROR:
+GUESS_ERROR:
     psFree(tmpX);
     psFree(tmpY);
@@ -344,2 +344,146 @@
     return corr;
 }
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+#define MEASURE_SAMPLES 4
+
+psImage *pmShutterCorrectionMeasure(const psVector *exptimes, // Exposure times
+                                    const psArray *images, // Input images
+                                    const psArray *masks, // Mask images
+                                    unsigned int size, // Size of samples
+                                    psStatsOptions meanStat, // Statistic to use for mean
+                                    psStatsOptions stdevStat, // Statistic to use for stdev
+                                    psMaskType maskVal // Mask value
+                                   )
+{
+    PS_ASSERT_VECTOR_NON_NULL(exptimes, NULL);
+    PS_ASSERT_VECTOR_TYPE(exptimes, PS_TYPE_F32, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(images, NULL);
+    if (masks) {
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(images, masks, NULL);
+    }
+    long num = images->n;               // Number of images
+    PS_ASSERT_VECTOR_SIZE(exptimes, num, NULL);
+
+    // Check input sizes, generate first-pass statistics
+    psRegion refRegion;                 // Reference region
+    psVector *refs = psVectorAlloc(num, PS_TYPE_F32); // Reference measurements
+    refs->n = num;
+    psVectorInit(refs, 0);
+    psArray *regions = psArrayAlloc(MEASURE_SAMPLES); // Array of sample regions, made on each image
+    regions->n = MEASURE_SAMPLES;
+    psImage *samplesMean = psImageAlloc(num, MEASURE_SAMPLES, PS_TYPE_F32); // Measurements for each file
+    psImage *samplesStdev = psImageAlloc(num, MEASURE_SAMPLES, PS_TYPE_F32); // Errors for each file
+    psStats *stats = psStatsAlloc(meanStat | stdevStat);
+    int numRows = 0, numCols = 0; // Size of images
+    for (long i = 0; i < images->n; i++) {
+        psImage *image = images->data[i]; // Image of interest
+        if (!image) {
+            continue;
+        }
+        if (numRows == 0 || numCols == 0) {
+            numRows = image->numRows;
+            numCols = image->numCols;
+            // Set up the sample regions
+            refRegion = psRegionForSquare(0.5 * numCols, 0.5 * numRows, size);
+            for (int j = 0; j < MEASURE_SAMPLES; j++) {
+                int x = (j % 2) ? size : image->numCols - size;
+                int y = (j > 1) ? size : image->numRows - size;
+                psRegion region = psRegionForSquare(x, y, size);
+                region = psRegionForImage(image, region);
+                regions->data[j] = psRegionAlloc(region.x0, region.x1, region.y0, region.y1);
+            }
+        } else if (numRows != image->numRows || numCols != image->numCols) {
+            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                    "Image sizes don't match: %dx%d vs %dx%d\n", image->numCols, image->numRows,
+                    numCols, numRows);
+            goto MEASURE_ERROR;
+        }
+
+        // Measure statistics
+        if (!psImageStats(stats, image, masks->data[i], maskVal)) {
+            psWarning("Unable to measure reference statistics.\n");
+        }
+        refs->data.F32[i] = psStatsGetValue(stats, meanStat);
+        if (refs->data.F32[i] <= 0.0) {
+            psError(PS_ERR_UNKNOWN, true, "Measured non-positive reference value.\n");
+            goto MEASURE_ERROR;
+        }
+        refs->data.F32[i] = 1.0 / refs->data.F32[i];
+        for (int j = 0; j < MEASURE_SAMPLES; j++) {
+            psRegion *region = regions->data[j]; // Region of interest
+            psImage *subImage = psImageSubset(image, *region); // Sub-image
+            if (!psImageStats(stats, subImage, masks->data[i], maskVal)) {
+                psString regionString = psRegionToString(*region);
+                psWarning("Unable to measure sample statistics at %s in image %ld.\n",
+                          regionString, i);
+                psFree(regionString);
+            }
+            psFree(subImage);
+            samplesMean->data.F32[j][i] = psStatsGetValue(stats, meanStat) * refs->data.F32[i];
+            samplesStdev->data.F32[j][i] = psStatsGetValue(stats, stdevStat) * refs->data.F32[i];
+        }
+    }
+    psFree(regions);
+    psFree(stats);
+
+    float meanRef = 0.0;                // Mean reference offset
+    int numGood = 0;                    // Number of good measurements
+    psVector *counts = psVectorAlloc(num, PS_TYPE_F32); // Mean for each image
+    psVector *errors = psVectorAlloc(num, PS_TYPE_F32); // Stdev for each image
+    for (int i = 0; i < MEASURE_SAMPLES; i++) {
+        counts = psImageRow(counts, samplesMean, i);
+        errors = psImageRow(errors, samplesStdev, i);
+        pmShutterCorrection *guess = pmShutterCorrectionGuess(exptimes, counts); // Guess at correction
+        pmShutterCorrection *params = pmShutterCorrectionFullFit(exptimes, counts, errors, guess); // Correc'n
+        if (isfinite(params->offref)) {
+            meanRef += params->offref;
+            numGood++;
+        }
+        psFree(params);
+        psFree(guess);
+    }
+    psFree(samplesMean);
+    psFree(samplesStdev);
+
+    if (numGood == 0) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to measure mean reference offset.\n");
+        psFree(counts);
+        psFree(errors);
+        goto MEASURE_ERROR;
+    }
+    meanRef /= (float)numGood;
+
+    psImage *shutter = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Shutter correction image
+    //psImage *pattern = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Illumination pattern
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            for (int i = 0; i < num; i++) {
+                psImage *image = images->data[i]; // Image of interest
+                counts->data.F32[i] = image->data.F32[y][x] * refs->data.F32[i];
+                errors->data.F32[i] = sqrtf(image->data.F32[y][x]) * refs->data.F32[i];
+            }
+
+            pmShutterCorrection *corr = pmShutterCorrectionLinFit(exptimes, counts, errors, meanRef);
+            shutter->data.F32[y][x] = corr->offset;
+            //pattern->data.F32[y][x] = corr->scale;
+            psFree(corr);
+        }
+    }
+    psFree(counts);
+    psFree(errors);
+    psFree(refs);
+
+    return shutter;
+
+
+MEASURE_ERROR:
+    // Clean up after error
+    psFree(refs);
+    psFree(regions);
+    psFree(stats);
+    psFree(samplesMean);
+    psFree(samplesStdev);
+    return NULL;
+}
Index: /trunk/psModules/src/detrend/pmShutterCorrection.h
===================================================================
--- /trunk/psModules/src/detrend/pmShutterCorrection.h	(revision 9314)
+++ /trunk/psModules/src/detrend/pmShutterCorrection.h	(revision 9315)
@@ -48,9 +48,12 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-05 20:48:56 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-10-05 23:36:26 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
+
+#ifndef PM_SHUTTER_CORRECTION_H
+#define PM_SHUTTER_CORRECTION_H
 
 #include "pslib.h"
@@ -88,2 +91,14 @@
         const pmShutterCorrection *guess // Initial guess
                                                );
+
+// Measure a shutter correction image from an array of images
+psImage *pmShutterCorrectionMeasure(const psVector *exptimes, // Exposure times
+                                    const psArray *images, // Input images
+                                    const psArray *masks, // Mask images
+                                    unsigned int size, // Size of samples
+                                    psStatsOptions meanStat, // Statistic to use for mean
+                                    psStatsOptions stdevStat, // Statistic to use for stdev
+                                    psMaskType maskVal // Mask value
+                                   );
+
+#endif