Index: /trunk/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 9297) +++ /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 9298) @@ -63,74 +63,102 @@ #include "psVectorBracket.h" -static void pmShutterCorrParsFree (pmShutterCorrPars *pars) -{ - if (pars == NULL) - return; +static void pmShutterCorrectionFree(pmShutterCorrection *pars) +{ + // Nothing to free 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) -pmShutterCorrPars *pmShutterCorrectionGuess (psVector *exptime, psVector *counts) -{ - - psVector *tmpX = NULL; - psVector *tmpY = NULL; - - pmShutterCorrPars *pars = pmShutterCorrParsAlloc (); - psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); - - int N = exptime->n; - // ASSERT (N >> 5) - +pmShutterCorrection *pmShutterCorrectionAlloc() +{ + pmShutterCorrection *corr = (pmShutterCorrection*)psAlloc(sizeof(pmShutterCorrection)); + psMemSetDeallocator(corr, (psFreeFunc)pmShutterCorrectionFree); + + corr->scale = 0.0; + corr->offset = 0.0; + corr->offref = 0.0; + + return corr; +} + +pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, const psVector *counts) +{ // 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. + + PS_ASSERT_VECTOR_NON_NULL(exptime, NULL); + PS_ASSERT_VECTOR_NON_NULL(counts, NULL); + PS_ASSERT_VECTOR_TYPE(exptime, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(counts, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(exptime, counts, NULL); + if (exptime->n <= 2) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Require more than 2 exposures to guess shutter correction.\n"); + return NULL; + } + + long N = exptime->n; // Number of exposures + + // use interpolation to guess shutter correction parameters given a set of exposures times and normalized + // counts (divided by the reference counts for each image) + + pmShutterCorrection *corr = pmShutterCorrectionAlloc(); // Shutter correction, to be returned + psPolynomial1D *line = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Straight line, for extrapolation // 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.F32[N-1]; + corr->scale = counts->data.F32[N-1]; // fit a line to the lowest three points and extrapolate to 0.0 - tmpX = psVectorAlloc (2, PS_TYPE_F32); - tmpY = psVectorAlloc (2, PS_TYPE_F32); + psVector *tmpX = psVectorAlloc(2, PS_TYPE_F32); + psVector *tmpY = psVectorAlloc(2, PS_TYPE_F32); tmpX->n = tmpY->n = 2; - tmpX->data.F32[0] = exptime->data.F32[0]; - tmpX->data.F32[1] = exptime->data.F32[1]; - tmpY->data.F32[0] = counts->data.F32[0]; - tmpY->data.F32[1] = counts->data.F32[1]; + long index; + for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++) + ; // Iterate only + if (index == N - 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Not enough good values to guess shutter correction.\n"); + goto ERROR; + } + tmpX->data.F32[0] = exptime->data.F32[index]; + tmpY->data.F32[0] = counts->data.F32[index]; + + for (index++; + (!isfinite(exptime->data.F32[index]) || exptime->data.F32[index] == exptime->data.F32[0]) && + index < N; index++) + ; // Iterate only + if (index == N) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Exposure times are all identical --- cannot guess shutter correction.\n"); + goto ERROR; + } + tmpY->data.F32[1] = counts->data.F32[index]; + tmpX->data.F32[1] = exptime->data.F32[index]; // fit a line and extrapolate the fit to 0.0 - psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX); - float ratio = psPolynomial1DEval (line, 0.0) / pars->scale; + if (!psVectorFitPolynomial1D(line, NULL, 0, tmpY, NULL, tmpX)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the time offset.\n"); + goto ERROR; + } + float ratio = psPolynomial1DEval(line, 0.0) / corr->scale; // XXX we need a sanity check: - // if the mean value of the three points is higher than pars->scale, + // if the mean value of the three points is higher than corr->scale, // then the slope should be negative. - // if the mean value of the three points is lower than pars->scale, + // if the mean value of the three points is lower than corr->scale, // then the slope should be positive. - // find two points bracketing the value counts = A (1 + dTk/dTo) / 2 = pars->scale (1 + ratio) / 2 - float value = pars->scale * (1 + ratio) / 2.0; - - int Np; + // find two points bracketing the value counts = A (1 + dTk/dTo) / 2 = corr->scale (1 + ratio) / 2 + float value = corr->scale * (1 + ratio) / 2.0; + + int Np; // Index of the value above (positive side) if (ratio < 1.0) { - Np = psVectorBracket (counts, value, true); + Np = psVectorBracket(counts, value, true); } else { - Np = psVectorBracketDescend (counts, value, true); - } - int Nm = (Np == 0) ? 1 : Np - 1; + Np = psVectorBracketDescend(counts, value, true); + } + int Nm = (Np == 0) ? 1 : Np - 1; // Index of the value below (negative side) tmpX->data.F32[0] = counts->data.F32[Nm]; @@ -140,25 +168,55 @@ // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo - 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); + if (!psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the reference offset.\n"); + goto ERROR; + } + corr->offref = psPolynomial1DEval(line, value); + corr->offset = ratio * corr->offref; + + psFree(line); + psFree(tmpX); + psFree(tmpY); + + return corr; + +ERROR: + psFree(tmpX); + psFree(tmpY); + psFree(line); + psFree(corr); + return NULL; } // linear fit to the counts and exptime, given a value for offref -pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, psVector *cntError, float offref) -{ +pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, + const psVector *counts, + const psVector *cntError, + float offref + ) +{ + PS_ASSERT_VECTOR_NON_NULL(exptime, NULL); + PS_ASSERT_VECTOR_NON_NULL(counts, NULL); + PS_ASSERT_VECTOR_TYPE(exptime, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(counts, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(exptime, counts, NULL); + if (exptime->n <= 2) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Require more than 2 exposures to guess shutter correction.\n"); + return NULL; + } + if (cntError) { + PS_ASSERT_VECTOR_TYPE(cntError, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(counts, cntError, NULL); + } + PS_ASSERT_FLOAT_LARGER_THAN(offref, 0.0, NULL); // this step is identical for all pixels: do it once and save? - psVector *x = psVectorAlloc (exptime->n, PS_TYPE_F32); - psVector *y = psVectorAlloc (exptime->n, PS_TYPE_F32); + psVector *x = psVectorAlloc(exptime->n, PS_TYPE_F32); + psVector *y = psVectorAlloc(exptime->n, PS_TYPE_F32); x->n = y->n = exptime->n; - for (int i = 0; i < exptime->n; i++) { + for (long i = 0; i < exptime->n; i++) { + // Should be safe (if expensive) to stick NaNs in --- the fitter deals with them float value = 1.0 / (exptime->data.F32[i] + offref); x->data.F32[i] = exptime->data.F32[i] * value; @@ -174,17 +232,22 @@ line->coeff[1][1] = 0; - psVectorFitPolynomial2D (line, NULL, 0, counts, cntError, x, y); - - pmShutterCorrPars *pars = pmShutterCorrParsAlloc (); - - pars->offref = offref; - pars->scale = line->coeff[1][0]; - pars->offset = line->coeff[0][1] / line->coeff[1][0]; - - psFree (x); - psFree (y); - psFree (line); - - return (pars); + if (!psVectorFitPolynomial2D(line, NULL, 0, counts, cntError, x, y)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit shutter correction.\n"); + psFree(x); + psFree(y); + psFree(line); + return NULL; + } + + pmShutterCorrection *corr = pmShutterCorrectionAlloc(); + corr->offref = offref; + corr->scale = line->coeff[1][0]; + corr->offset = line->coeff[0][1] / line->coeff[1][0]; + + psFree(x); + psFree(y); + psFree(line); + + return corr; } @@ -193,24 +256,45 @@ const psVector *x) { - psF32 A = params->data.F32[0]; - psF32 p = x->data.F32[0] + params->data.F32[1]; - psF32 q = 1.0 / (x->data.F32[0] + params->data.F32[2]); - psF32 f = A*p*q; - - if (deriv != NULL) { - deriv->data.F32[0] = p*q; - deriv->data.F32[1] = A*q; - deriv->data.F32[2] = -A*p*q*q; - } - return(f); + // This is in a tight loop, so we won't assert here. + + psF32 A = params->data.F32[0]; + psF32 p = x->data.F32[0] + params->data.F32[1]; + psF32 q = 1.0 / (x->data.F32[0] + params->data.F32[2]); + psF32 f = A * p * q; + + if (deriv) { + deriv->data.F32[0] = p * q; + deriv->data.F32[1] = A * q; + deriv->data.F32[2] = - f * q; + } + return f; } // non-linear fit to the counts and exptime, given a guess for the three parameters -pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, psVector *cntError, pmShutterCorrPars *guess) -{ - - psMinimization *myMin = psMinimizationAlloc(15, 0.1); - - psVector *params = psVectorAlloc (3, PS_TYPE_F32); +pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, + const psVector *counts, + const psVector *cntError, + const pmShutterCorrection *guess + ) +{ + PS_ASSERT_VECTOR_NON_NULL(exptime, NULL); + PS_ASSERT_VECTOR_NON_NULL(counts, NULL); + PS_ASSERT_VECTOR_TYPE(exptime, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(counts, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(exptime, counts, NULL); + if (exptime->n <= 2) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Require more than 2 exposures to guess shutter correction.\n"); + return NULL; + } + if (cntError) { + PS_ASSERT_VECTOR_TYPE(cntError, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(counts, cntError, NULL); + } + PS_ASSERT_PTR_NON_NULL(guess, NULL); + + psMinimization *minInfo = psMinimizationAlloc(15, 0.1); // Minimization information + + psVector *params = psVectorAlloc (3, PS_TYPE_F32); // Fitting parameters params->data.F32[0] = guess->scale; params->data.F32[1] = guess->offset; @@ -224,15 +308,15 @@ // constrain->paramMax = psVectorAlloc (3, PS_TYPE_F32); // constrain->paramMask = NULL; - psMinConstrain *constrain = NULL; + psMinConstrain *constrain = NULL; // Constraints on the minimization // XXX ignore covariance matrix for the moment // psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psImage *covar = NULL; + psImage *covar = NULL; // Covariance matrix // construct the coordinate and value entries (y is counts) - psArray *x = psArrayAlloc(exptime->n); + psArray *x = psArrayAlloc(exptime->n); // Coordinates x->n = exptime->n; - for (int i = 0; i < exptime->n; i++) { + for (long i = 0; i < exptime->n; i++) { psVector *coord = psVectorAlloc(1, PS_TYPE_F32); coord->n = 1; @@ -241,16 +325,21 @@ } - psMinimizeLMChi2(myMin, covar, params, constrain, x, counts, cntError, pmShutterCorrectionModel); - - pmShutterCorrPars *pars = pmShutterCorrParsAlloc (); - - pars->scale = params->data.F32[0]; - pars->offset = params->data.F32[1]; - pars->offref = params->data.F32[2]; - - psFree (myMin); - psFree (params); - psFree (x); - - return (pars); -} + if (!psMinimizeLMChi2(minInfo, covar, params, constrain, x, counts, cntError, pmShutterCorrectionModel)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for shutter correction.\n"); + psFree(x); + psFree(minInfo); + psFree(params); + return NULL; + } + + pmShutterCorrection *corr = pmShutterCorrectionAlloc(); // Shutter correction + corr->scale = params->data.F32[0]; + corr->offset = params->data.F32[1]; + corr->offref = params->data.F32[2]; + + psFree(minInfo); + psFree(params); + psFree(x); + + return corr; +} Index: /trunk/psModules/src/detrend/pmShutterCorrection.h =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 9297) +++ /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 9298) @@ -31,8 +31,8 @@ * 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)*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. + * values by interpolation and/or extrapolation. * To take the strategy one step further, we could use the above recipe to @@ -48,6 +48,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-09-25 01:06:28 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-05 20:48:56 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -56,14 +56,34 @@ #include "pslib.h" +// Shutter correction parameters, applicable for a single pixel typedef struct { - double scale; // A(k) - double offset; // dTk - double offref; // dTo + double scale; // The normalisation for an exposure, A(k) + double offset; // The time offset, dTk + double offref; // The reference time offset, dTo } -pmShutterCorrPars; +pmShutterCorrection; -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); +// Allocator +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 + ); + +// 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. + float offref // Reference time offset + ); + +// 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 + );