Changeset 9618

Timestamp:

Oct 17, 2006, 2:47:15 PM (20 years ago)

Author:

Paul Price

Message:

Documenting pmShutterCorrection.[ch]

Location:

trunk/psModules/src/detrend

Files:

• pmShutterCorrection.c (modified) (5 diffs)
• pmShutterCorrection.h (modified) (1 diff)

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -1 @@
-/*
- * 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
- 
-*/
-
 #if HAVE_CONFIG_H
 #include <config.h>
@@ -63 +11 @@
 #include "pmShutterCorrection.h"
 #include "psVectorBracket.h"
+
+/// 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.
+///
+/// The process for generating a shutter correction is as follows:
+/// - 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
+
 
 static void pmShutterCorrectionFree(pmShutterCorrection *pars)
@@ -360 +350 @@
 #define MEASURE_SAMPLES 4
 
-psImage *pmShutterCorrectionMeasure(const psVector *exptimes, // Exposure times
-                                    const psArray *images, // Input images
-                                    const psArray *weights, // Weight 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
-                                    int nIter, // Number of iterations per pixel
-                                    float rej, // Rejection limit
-                                    psMaskType maskVal // Mask value
-                                   )
+psImage *pmShutterCorrectionMeasure(const psVector *exptimes, const psArray *images, const psArray *weights,
+                                    const psArray *masks, unsigned int size, psStatsOptions meanStat,
+                                    psStatsOptions stdevStat, int nIter, float rej, psMaskType maskVal)
 {
     PS_ASSERT_VECTOR_NON_NULL(exptimes, NULL);
@@ -376 +358 @@
     PS_ASSERT_ARRAY_NON_NULL(images, NULL);
     if (masks) {
-        // XXX ASSERT not defined
-        // PS_ASSERT_ARRAYS_SIZE_EQUAL(images, masks, NULL);
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(images, masks, NULL);
     }
     if (weights) {
-        // XXX ASSERT not defined
-        // PS_ASSERT_ARRAYS_SIZE_EQUAL(images, weights, NULL);
+        PS_ASSERT_ARRAYS_SIZE_EQUAL(images, weights, NULL);
     }
     long num = images->n;               // Number of images
@@ -576 +556 @@
 
 
-bool pmShutterCorrectionApply(pmReadout *readout, // Readout to which to apply shutter correction
-                              const pmReadout *shutter // Shutter correction readout
-                             )
+bool pmShutterCorrectionApply(pmReadout *readout, const pmReadout *shutter)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);

trunk/psModules/src/detrend/pmShutterCorrection.h

@@ -1 @@
-/** @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.8 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-06 03:33:49 $
- *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- */
+/// @file pmShutterCorrection.h
+///
+/// @brief Functions to build and apply a shutter exposure-time correction.
+///
+/// @ingroup Detrend
+///
+/// @author Eugene Magnier, IfA
+/// @author Paul Price, IfA
+///
+/// @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2006-10-18 00:47:15 $
+///
+/// Copyright 2006 Institute for Astronomy, University of Hawaii
+///
+
+
+/// 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)
 
 #ifndef PM_SHUTTER_CORRECTION_H
 #define PM_SHUTTER_CORRECTION_H
 
-#include "pslib.h"
+#include <pslib.h>
 
-// Shutter correction parameters, applicable for a single pixel
+/// Shutter correction parameters, applicable for a single pixel
 typedef struct
 {
-    double scale;                       // The normalisation for an exposure, A(k)
-    double offset;                      // The time offset, dTk
-    double offref;                      // The reference time offset, dTo
+    double scale;                       ///< The normalisation for an exposure, A(k)
+    double offset;                      ///< The time offset, dTk
+    double offref;                      ///< The reference time offset, dTo
 }
 pmShutterCorrection;
 
-// Allocator
+/// Allocator for shutter correction parameters
 pmShutterCorrection *pmShutterCorrectionAlloc();
 
-// Guess a shutter correction, based on plot of counts vs exposure time.
-// Used before doing the full non-linear fit, to get parameters close to the true.
-// Assumes exptime vector is sorted (ascending order; longest is last) prior to input.
-pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, // Exposure times for each exposure
-        const psVector *counts // Counts for each exposure
+/// Guess a shutter correction, based on plot of counts vs exposure time
+///
+/// This function is used before doing the full non-linear fit, to get parameters close to the true.  Assumes
+/// exptime vector is sorted (ascending order; longest is last) prior to input.
+pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, ///< Exposure times for each exposure
+        const psVector *counts ///< Counts for each exposure
                                              );
 
-// Generate shutter correction based on a linear fit
-pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, // Exposure times for each exposure
-        const psVector *counts, // Counts for each exposure
-        const psVector *cntError, // Error in the counts, for each exp.
-        const psVector *mask, // Mask for each exposure
-        float offref, // Reference time offset
-        int nIter, // Number of iterations
-        float rej, // Rejection threshold (sigma)
-        psMaskType maskVal // Mask value
+/// Generate shutter correction based on a linear fit
+///
+/// Performs a linear fit to counts as a function of exposure time, with the reference time offset fixed (so
+/// that the system is linear).  Performs iterative clipping, if nIter > 1.
+pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, ///< Exposure times for each exposure
+        const psVector *counts, ///< Counts for each exposure
+        const psVector *cntError, ///< Error in the counts
+        const psVector *mask, ///< Mask for each exposure
+        float offref, ///< Reference time offset
+        int nIter, ///< Number of iterations
+        float rej, ///< Rejection threshold (sigma)
+        psMaskType maskVal ///< Mask value
                                               );
 
-// Generate shutter correction based on a full non-linear fit
-pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, // Exposure times for each exposure
-        const psVector *counts, // Counts for each exposure
-        const psVector *cntError, // Error in the counts, for each exp
-        const pmShutterCorrection *guess // Initial guess
+/// Generate shutter correction based on a full non-linear fit
+///
+/// Performs a full non-linear fit to counts as a function of exposure time.  The main purpose is to solve for
+/// the reference time offset, so that future fits may be performed using linear fitting with the reference
+/// time offset fixed.
+pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, ///< Exposure times for each exposure
+        const psVector *counts, ///< Counts for each exposure
+        const psVector *cntError, ///< Error in the counts
+        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 *weights, // Weight 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
-                                    int nIter, // Number of iterations
-                                    float rej, // Rejection threshold (sigma)
-                                    psMaskType maskVal // Mask value
+/// Measure a shutter correction image from an array of images
+///
+/// Given an array of images with known exposure times, this function measures the shutter correction (our
+/// principal concern is for the time offset, rather than the normalisation) by measuring the reference time
+/// offset using the full non-linear fit for a small number of representative regions (middle and corners),
+/// and then using that to perform a linear fit to each pixel.
+psImage *pmShutterCorrectionMeasure(const psVector *exptimes, ///< Exposure times
+                                    const psArray *images, ///< Input images
+                                    const psArray *weights, ///< Weight images
+                                    const psArray *masks, ///< Mask images
+                                    unsigned int size, ///< Size of samples for statistics for non-linear fit
+                                    psStatsOptions meanStat, ///< Statistic to use for mean
+                                    psStatsOptions stdevStat, ///< Statistic to use for stdev
+                                    int nIter, ///< Number of iterations
+                                    float rej, ///< Rejection threshold (sigma)
+                                    psMaskType maskVal ///< Mask value
                                    );
 
-// Apply a shutter correction
-bool pmShutterCorrectionApply(pmReadout *readout, // Readout to which to apply shutter correction
-                              const pmReadout *shutter // Shutter correction readout
+/// Apply a shutter correction
+///
+/// Given a shutter correction (with dT for each pixel), applies this correction to an input image.
+bool pmShutterCorrectionApply(pmReadout *readout, ///< Readout to which to apply shutter correction
+                              const pmReadout *shutter ///< Shutter correction readout, with dT for each pixel
                              );