Changeset 8866 for trunk/psModules/src/detrend/pmShutterCorrection.c

Timestamp:

Sep 20, 2006, 4:18:19 PM (20 years ago)

Author:

magnier

Message:

working on shutter corrections

File:

• trunk/psModules/src/detrend/pmShutterCorrection.c (modified) (4 diffs)

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -51 +51 @@
 */
 
+static void pmShutterCorrParsFree (pmShutterCorrPars *pars)
+{
+    if (pars == NULL)
+        return;
+    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)
@@ -62 +80 @@
 
     N = exptime->n;
+    // ASSERT (N >> 5)
 
-    // 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.
+    // 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.
 
     // choose the highest exptime point as the guess for the scale:
@@ -72 +91 @@
 
     // fit a line to the lowest three points and extrapolate to 0.0
-    tmpX = psVectorRecycle (tmpX, 3, PS_TYPE_F64);
-    tmpY = psVectorRecycle (tmpY, 3, PS_TYPE_F64);
+    tmpX = psVectorAlloc (2, PS_TYPE_F64);
+    tmpY = psVectorAlloc (2, PS_TYPE_F64);
 
-    tmpX->data.F64[0] = exptime->data.F64[N-3];
-    tmpX->data.F64[1] = exptime->data.F64[N-2];
-    tmpX->data.F64[2] = exptime->data.F64[N-1];
-    tmpY->data.F64[0] = counts->data.F64[N-3];
-    tmpY->data.F64[1] = counts->data.F64[N-2];
-    tmpY->data.F64[2] = counts->data.F64[N-1];
+    tmpX->data.F64[0] = exptime->data.F64[0];
+    tmpX->data.F64[1] = exptime->data.F64[1];
+    tmpY->data.F64[0] = counts->data.F64[0];
+    tmpY->data.F64[1] = counts->data.F64[1];
 
-    // fit a line to use for interpolation
+    // fit a line and extrapolate the fit to 0.0
     psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     line = psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX);
+    ratio = psPolynomial1DEval (line, 0.0) / pars->scale;
 
     // XXX we need a sanity check:
@@ -92 +110 @@
     // then the slope should be positive.
 
-    // extrapolate the fit to 0.0
-    ratio = psPolynomial1DEval (line, 0.0);
+    // find two points bracketing the value counts = A (1 + dTk/dTo) / 2 = pars->scale (1 + ratio) / 2
+    value = pars->scale * (1 + ratio) / 2.0;
 
+    Nm = psVectorBracket (exptime, value, (ratio < 1.0));
+    Np = (Nm == N - 1) ? Nm - 1 : Nm + 1;
+
+    tmpX->data.F64[0] = counts->data.F64[Nm];
+    tmpX->data.F64[1] = counts->data.F64[Np];
+    tmpY->data.F64[0] = exptime->data.F64[Nm];
+    tmpY->data.F64[1] = exptime->data.F64[Np];
+
+    // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
+    psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+    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);
 }