Index: /trunk/psModules/src/detrend/pmShutterCorrection.c
===================================================================
--- /trunk/psModules/src/detrend/pmShutterCorrection.c	(revision 8865)
+++ /trunk/psModules/src/detrend/pmShutterCorrection.c	(revision 8866)
@@ -51,4 +51,22 @@
 */
 
+static void pmShutterCorrParsFree (pmShutterCorrPars *pars)
+{
+    if (pars == NULL)
+        return;
+    return;
+}
+
+pmShutterCorrPars *pmShutterCorrParsAlloc ()
+{
+
+    pmShutterCorrPars *pars = (pmShutterCorrPars *) psAlloc(sizeof(pmShutterCorrPars));
+    psMemSetDeallocator(pars, (psFreeFunc) pmShutterCorrParsFree);
+    pars->scale  = 0.0;
+    pars->offset = 0.0;
+    pars->offref = 0.0;
+    return (pars);
+}
+
 // use interpolation to guess shutter correction parameters given a set of exposures times and normalized
 // counts (divided by the reference counts for each image)
@@ -62,7 +80,8 @@
 
     N = exptime->n;
+    // ASSERT (N >> 5)
 
-    // NOTE: vectors must be sorted on input it is expensive to sort or check this here, but it is easy to
-    // arrange by sorting the images before generating these vectors.
+    // NOTE: vectors must be sorted on input.  It is expensive to sort or check this here, but
+    // it is easy to arrange by sorting the images before generating these vectors.
 
     // choose the highest exptime point as the guess for the scale:
@@ -72,17 +91,16 @@
 
     // fit a line to the lowest three points and extrapolate to 0.0
-    tmpX = psVectorRecycle (tmpX, 3, PS_TYPE_F64);
-    tmpY = psVectorRecycle (tmpY, 3, PS_TYPE_F64);
+    tmpX = psVectorAlloc (2, PS_TYPE_F64);
+    tmpY = psVectorAlloc (2, PS_TYPE_F64);
 
-    tmpX->data.F64[0] = exptime->data.F64[N-3];
-    tmpX->data.F64[1] = exptime->data.F64[N-2];
-    tmpX->data.F64[2] = exptime->data.F64[N-1];
-    tmpY->data.F64[0] = counts->data.F64[N-3];
-    tmpY->data.F64[1] = counts->data.F64[N-2];
-    tmpY->data.F64[2] = counts->data.F64[N-1];
+    tmpX->data.F64[0] = exptime->data.F64[0];
+    tmpX->data.F64[1] = exptime->data.F64[1];
+    tmpY->data.F64[0] = counts->data.F64[0];
+    tmpY->data.F64[1] = counts->data.F64[1];
 
-    // fit a line to use for interpolation
+    // fit a line and extrapolate the fit to 0.0
     psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     line = psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX);
+    ratio = psPolynomial1DEval (line, 0.0) / pars->scale;
 
     // XXX we need a sanity check:
@@ -92,6 +110,25 @@
     // then the slope should be positive.
 
-    // extrapolate the fit to 0.0
-    ratio = psPolynomial1DEval (line, 0.0);
+    // find two points bracketing the value counts = A (1 + dTk/dTo) / 2 = pars->scale (1 + ratio) / 2
+    value = pars->scale * (1 + ratio) / 2.0;
 
+    Nm = psVectorBracket (exptime, value, (ratio < 1.0));
+    Np = (Nm == N - 1) ? Nm - 1 : Nm + 1;
+
+    tmpX->data.F64[0] = counts->data.F64[Nm];
+    tmpX->data.F64[1] = counts->data.F64[Np];
+    tmpY->data.F64[0] = exptime->data.F64[Nm];
+    tmpY->data.F64[1] = exptime->data.F64[Np];
+
+    // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
+    psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+    line = psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX);
+    pars->offref = psPolynomial1DEval (line, value);
+    pars->offset = ratio * pars->offref;
+
+    psFree (line);
+    psFree (tmpX);
+    psFree (tmpY);
+
+    return (pars);
 }
Index: /trunk/psModules/src/detrend/pmShutterCorrection.h
===================================================================
--- /trunk/psModules/src/detrend/pmShutterCorrection.h	(revision 8865)
+++ /trunk/psModules/src/detrend/pmShutterCorrection.h	(revision 8866)
@@ -28,5 +28,5 @@
  * we avoid directly fitting these values as the process would be a non-linear
  * least-squares problem for every pixel in the image, and thus very time
- * consumning.  There are instead linear options which may be used instead.
+ * consuming.  There are linear options which may be used instead.
  * First, as T_o goes to a large value, s() approaches the value of f'(x,y).
  * Next, as T_o goes to a very small value, s() approaches the value of
@@ -48,6 +48,6 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-09-20 01:03:00 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-09-21 02:18:19 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -64,14 +64,9 @@
  */
 
-bool pmFlatField(
-    pmReadout *in,          ///< Readout with input image
-    const pmReadout *flat   ///< Readout with flat image
-);
-
 typedef struct
 {
-    double scale;
-    double offset;
-    double offref;
+    double scale;     // A(k)
+    double offset;   // dTk
+    double offref;   // dTo
 }
 pmShutterCorrPars;