Changeset 9298

Timestamp:

Oct 5, 2006, 10:48:56 AM (20 years ago)

Author:

Paul Price

Message:

Cleaning up code, adding error checks, adding const where appropriate.

Location:

trunk/psModules/src/detrend

Files:

• pmShutterCorrection.c (modified) (6 diffs)
• pmShutterCorrection.h (modified) (3 diffs)

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -63 +63 @@
 #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 +168 @@
 
     // 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 +232 @@
     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 +256 @@
                                       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 +308 @@
     // 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 +325 @@
     }
 
-    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;
+}

trunk/psModules/src/detrend/pmShutterCorrection.h

@@ -31 +31 @@
  * 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 +48 @@
  *  @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 +56 @@
 #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
+                                               );