Index: trunk/psModules/src/detrend/pmShutterCorrection.c
===================================================================
--- trunk/psModules/src/detrend/pmShutterCorrection.c	(revision 8842)
+++ trunk/psModules/src/detrend/pmShutterCorrection.c	(revision 8842)
@@ -0,0 +1,97 @@
+/*
+ * measure shutter correction:
+ 
+ * input  : collection of shutter correction exposures (pre-processed)
+ * output : a shutter correction image
+ 
+ * the measurement could be performed on any focal-plane unit at a time. for
+ * GPC, the obvious scale is to measure the effect on the entire focal plane at
+ * once, with a single reference point in the field.  this is a little more
+ * complex than just measuring the effect for a single 2D image array.  the
+ * reference point and the detailed analysis points need to be defined for the
+ * entire hierarchy rather than just as coordinate pairs or regions.  a
+ * pmFPAview would be a natural element with which to define these points, but
+ * at the moment, the pmFPAview structure defines a band in the CCD, not a
+ * coordinate.  An option is to instead specify the reference locations as a
+ * pmFPAview coupled with a psRegion, though we need to be careful not to
+ * over-specify the pixels (ie, conflict between pmFPAview and psRegion).
+ 
+ * at each point in an image with exposure time T, we measure f(k;T) = F(k;T) /
+ * F(0;T) where k is the coordinate of the point of interest, 0 is the reference
+ * coordinate, and F(k;T) is the measured number of counts at the point of
+ * interest in this image.  given a collection of f(k;T) values, we need to
+ * determine the model f(k;T) = A(k) (T + dTk) / (T + dTo) where dTk is the shutter
+ * error at the given position, dTo is the shutter error at the reference
+ * position, and A(k) is the scaling factor for the given position.
+ 
+ - load the configuration information
+ - load the input file information, build pmFPAfile's for all inputs
+ 
+ - determine the reference point and detailed analysis regions from the config information
+ - for each image
+ -- measure the reference point counts
+ - for each analysis region:
+ -- measure shutter parameters (dTo, dTk, A):
+ --- for each image:
+ ---- measure counts at the region
+ ---- divide by the reference counts
+ --- linear extrapolation to find f(inf) = A(k)
+ --- linear extrapolation to find f(0) = A(k) dTk / dTo
+ --- linear interpolation to find coordinate where f(dTo) = A (1 + dTk/dTo) / 2
+ --- non-linear fit of T, f(T) to f(k;T) = A(k) (T + dTk) / (T + dTo)
+ - use the collection of dTo values to choose a best value for dTo (median)
+ - for each image pixel
+ -- divide by the reference counts
+ -- generate the vectors T, f(T)
+ -- linear fit of T, f(T) to f(k;T) = A(k) (T + dTk) / (T + dTo) using dTo above
+ -- save dTk, A(k) in output image pixels
+ -- apply dTk, A(k) to measure residual images
+ -- generate residual FITS/JPEG images
+ 
+*/
+
+// use interpolation to guess shutter correction parameters given a set of exposures times and normalized
+// counts (divided by the reference counts for each image)
+pmShutterCorrPars *pmShutterCorrectionGuess (psVector *exptime, psVector *counts)
+{
+
+    psVector *tmpX = NULL;
+    psVector *tmpY = NULL;
+
+    pmShutterCorrPars *pars = pmShutterCorrParsAlloc ();
+
+    N = exptime->n;
+
+    // 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:
+    // XXX we could examine the top 2 or 3 values and decide if we
+    // extended exptime enough or median clip.
+    pars->scale = counts->data.F64[N-1];
+
+    // 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->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];
+
+    // fit a line to use for interpolation
+    psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+    line = psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX);
+
+    // XXX we need a sanity check:
+    // if the mean value of the three points is higher than pars->scale,
+    // then the slope should be negative.
+    // if the mean value of the three points is lower than pars->scale,
+    // then the slope should be positive.
+
+    // extrapolate the fit to 0.0
+    ratio = psPolynomial1DEval (line, 0.0);
+
+}
Index: trunk/psModules/src/detrend/pmShutterCorrection.h
===================================================================
--- trunk/psModules/src/detrend/pmShutterCorrection.h	(revision 8842)
+++ trunk/psModules/src/detrend/pmShutterCorrection.h	(revision 8842)
@@ -0,0 +1,77 @@
+/** @file  pmShutterCorrection.h
+ *
+ *  @brief Functions to build and apply a shutter exposure-time correction.
+ *
+ *  A mechanical shutter may not yield uniform exposure times as a function of
+ *  position on the detector.  The typical error consists of a constant
+ *  exposure-time offset relative to the requested value, ie exposure time is
+ *  T_o + dT(x,y).  The exposure error, dT, may be measured with the following
+ *  scheme.  Obtain a set of exposures with different exposures times taken of
+ *  the same flat-field source; the source must be spatially stable between the
+ *  exposures, but need not have a stable amplitude.  For an illuminating flux
+ *  of intensity F(x,y) = F_o f(x,y), the signal recorded by any pixel in the
+ *  detector is given by: S(t,x,y) = F_o(t) f(x,y) (T_o + dT(x,y)) where F_o is
+ *  the F_o(t) is the (variable) overall intensity of the illuminating source
+ *  and f(x,y) is the spatial illumination patter times the flat-field response.
+ *  Choose a reference location in the image (eg, the detector center) and
+ *  divide by the value of that region (ie, mean or median):
+ 
+ *  s(t,x,y) = S(t,x,y) / S(t,0,0)
+ *  s(t,x,y) = F_o(t) f(x,y) (T_o + dT(x,y)) / F_o(t) f(0,0) (T_o + dT(0,0))
+ *  s(t,x,y) = f(x,y) (T_o + dT(x,y)) / f(0,0) (T_o + dT(0,0))
+ 
+ * we can absorb the term f(0,0) into f(x,y) as we have no motivation for the
+ * scale of f(x,y).  For any single pixel, over the set of exposures, we thus
+ * need to solve for dT(x,y), dT(0,0), and f'(x,y) in the equation:
+ * s(t,x,y) = f'(x,y) (T_o + dT(x,y)) / (T_o + dT(0,0))
+ 
+ * 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.
+ * 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
+ * f'(x,y)*dT(x,y)/dT(0,0).  Finally, when s() has the value of 
+ * f'(x,y)*(1 + dT(x,y)/dT(0,0))/2, T_o has the value of dT(0,0).  with data
+ * points covering a reasonable dynamic range, we can solve for these three
+ * values by interpolation and/or extrapolation.  
+ 
+ * To take the strategy one step further, we could use the above recipe to
+ * obtain a guess for the three parameters and then apply non-linear fitting to
+ * solve more accurately for the parameters.  If we limit this operation to a
+ * handful of positions in the image (user defined, but the obvious choice would
+ * be positions near the center, edges, and corners), then we may determine a
+ * good value for dT(0,0).  Since there is only one dT(0,0) for the image, we
+ * can apply the resulting measurement to the rest of the pixels in the image.
+ * If dT(0,0) is not a free parameter, then the fitting process is linear in
+ * terms of dT(x,y) and f'(x,y)
+ 
+ *  @author Eugene Magnier, IfA
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-09-20 01:03:00 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#include "pslib.h"
+
+/** Execute flat field module.
+ *
+ *  Given an input image and a flat-field image, pmFlatField shall divide the input image by the flat field
+ *  image.
+ *
+ *  @return  bool: True or false for success or failure
+ */
+
+bool pmFlatField(
+    pmReadout *in,          ///< Readout with input image
+    const pmReadout *flat   ///< Readout with flat image
+);
+
+typedef struct
+{
+    double scale;
+    double offset;
+    double offref;
+}
+pmShutterCorrPars;