Index: /trunk/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 9617) +++ /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 9618) @@ -1,54 +1,2 @@ -/* - * 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 @@ -63,4 +11,46 @@ #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,15 +350,7 @@ #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,10 +358,8 @@ 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,7 +556,5 @@ -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); Index: /trunk/psModules/src/detrend/pmShutterCorrection.h =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 9617) +++ /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 9618) @@ -1,115 +1,134 @@ -/** @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 -// 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 );