/** @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
 * 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
 * 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.4 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-09-25 01:06:28 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

#include "pslib.h"

typedef struct
{
    double scale;     // A(k)
    double offset;   // dTk
    double offref;   // dTo
}
pmShutterCorrPars;

pmShutterCorrPars *pmShutterCorrParsAlloc ();
pmShutterCorrPars *pmShutterCorrectionGuess (psVector *exptime, psVector *counts);
pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, psVector *cntError, float offref);
pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, psVector *cntError, pmShutterCorrPars *guess);
