Changeset 27839 for branches/simmosaic_branches/psModules/src/detrend/pmShutterCorrection.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/detrend/pmShutterCorrection.c (modified) (9 diffs)

branches/simmosaic_branches/psModules/src/detrend/pmShutterCorrection.c

@@ -64 +64 @@
                                         // the code (see pmShutterCorrectionDataAlloc).
 
-//int corrRefCount = 1;
-//int corrFitCount = 1;
-//int corrGuessCount = 1;
 
 static void pmShutterCorrectionFree(pmShutterCorrection *pars)
@@ -91 +88 @@
 pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, const psVector *counts)
 {
-//printf("pmShutterCorrectGuess count: %d\n", corrGuessCount++);
     // 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) {
-        printf("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.
-    corr->scale = counts->data.F32[N-1];
-
-    // fit a line to the lowest three points and extrapolate to 0.0
-    psVector *tmpX = psVectorAlloc(2, PS_TYPE_F32);
-    psVector *tmpY = psVectorAlloc(2, PS_TYPE_F32);
-
-    long index;
-
-    // Iterate only
-    for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++);
-
-    if (index == N - 1) {
-        printf("Not enough good values to guess shutter correction.\n");
-        goto GUESS_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) {
-        printf("Exposure times are all identical --- cannot guess shutter correction.\n");
-        goto GUESS_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
-    if (!psVectorFitPolynomial1D(line, NULL, 0, tmpY, NULL, tmpX)) {
-        printf("Unable to fit for the time offset.\n");
-        goto GUESS_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 corr->scale,
-    // then the slope should be negative.
-    // 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 = 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);
-    } else {
-        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];
-    tmpX->data.F32[1] = counts->data.F32[Np];
-    tmpY->data.F32[0] = exptime->data.F32[Nm];
-    tmpY->data.F32[1] = exptime->data.F32[Np];
-
-    // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
-    if (!psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX)) {
-        printf("Unable to fit for the reference offset.\n");
-        printf("tmpX0: %f, tmpX1: %f, tmpY0: %f, tmpY1: %f\n", tmpX->data.F32[0], tmpX->data.F32[1], tmpY->data.F32[0], tmpY->data.F32[1]);
-        goto GUESS_ERROR;
-    }
-    corr->offref = psPolynomial1DEval(line, value);
-    corr->offset = ratio * corr->offref;
-printf("guess scale, offset, offref: %f, %f, %f\n", corr->scale, corr->offset, corr->offref);
-    psFree(line);
-    psFree(tmpX);
-    psFree(tmpY);
-
-    return corr;
-
-GUESS_ERROR:
-    psFree(tmpX);
-    psFree(tmpY);
-    psFree(line);
-//    psFree(corr);
-//    return NULL;
-    corr->scale = 0.5;
-    corr->offset = 0.5;
-    corr->offref = 0.5;
-    return corr;
-}
-
-// linear fit to the counts and exptime, given a value for offref
-pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, const psVector *counts,
-                                               const psVector *cntError, const psVector *mask, float offref,
-                                               int nIter, float rej)
-{
     PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
     PS_ASSERT_VECTOR_NON_NULL(counts, NULL);
@@ -216 +101 @@
         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);
-
-    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;
-        y->data.F32[i] = value;
-    }
-
-    psPolynomial2D *line = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
-
-    // mask out the terms we will not fit
-    line->coeffMask[0][0] = PS_POLY_MASK_SET;
-    line->coeffMask[1][1] = PS_POLY_MASK_SET;
-    line->coeff[0][0] = 0;
-    line->coeff[1][1] = 0;
-
-    // the stats structure determines how the clipping statistic is measured
-    // too few points to use the robust analysis method
-    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    stats->clipSigma = rej;
-    stats->clipIter = nIter;
-
-    if (!psVectorClipFitPolynomial2D(line, stats, mask, 0xff, counts, cntError, x, y)) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit shutter correction.\n");
-        psFree(stats);
-        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];
-    corr->num = stats->clippedNvalues;
-    corr->stdev = stats->clippedStdev;
-
-    psFree(stats);
-    psFree(x);
-    psFree(y);
+
+    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.
+    corr->scale = counts->data.F32[N-1];
+
+    // fit a line to the lowest three points and extrapolate to 0.0
+    psVector *tmpX = psVectorAlloc(2, PS_TYPE_F32);
+    psVector *tmpY = psVectorAlloc(2, PS_TYPE_F32);
+
+    long index;
+
+    // Iterate only
+    for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++);
+
+    if (index == N - 1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Not enough good values to guess shutter correction.\n");
+        goto GUESS_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 GUESS_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
+    if (!psVectorFitPolynomial1D(line, NULL, 0, tmpY, NULL, tmpX)) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the time offset.\n");
+        goto GUESS_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 corr->scale,
+    // then the slope should be negative.
+    // 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 = 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);
+    } else {
+        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];
+    tmpX->data.F32[1] = counts->data.F32[Np];
+    tmpY->data.F32[0] = exptime->data.F32[Nm];
+    tmpY->data.F32[1] = exptime->data.F32[Np];
+
+    // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
+    if (!psVectorFitPolynomial1D (line, NULL, 0, tmpY, NULL, tmpX)) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit for the reference offset.\n");
+        goto GUESS_ERROR;
+    }
+    corr->offref = psPolynomial1DEval(line, value);
+    corr->offset = ratio * corr->offref;
+
     psFree(line);
+    psFree(tmpX);
+    psFree(tmpY);
 
     return corr;
-}
-
-static psF32 pmShutterCorrectionModel(psVector *deriv, const psVector *params, const psVector *x)
-{
-    // 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
-pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, const psVector *counts,
-                                                const psVector *cntError, const pmShutterCorrection *guess)
-{
-//printf("pmShutterCorrectionFullFit count: %d\n", corrFitCount++);
+
+GUESS_ERROR:
+    psFree(tmpX);
+    psFree(tmpY);
+    psFree(line);
+    psFree(corr);
+    return NULL;
+}
+
+// linear fit to the counts and exptime, given a value for offref
+pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, const psVector *counts,
+                                               const psVector *cntError, const psVector *mask, float offref,
+                                               int nIter, float rej)
+{
     PS_ASSERT_VECTOR_NON_NULL(exptime, NULL);
     PS_ASSERT_VECTOR_NON_NULL(counts, NULL);
@@ -307 +214 @@
         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);
+
+    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;
+        y->data.F32[i] = value;
+    }
+
+    psPolynomial2D *line = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
+
+    // mask out the terms we will not fit
+    line->coeffMask[0][0] = PS_POLY_MASK_SET;
+    line->coeffMask[1][1] = PS_POLY_MASK_SET;
+    line->coeff[0][0] = 0;
+    line->coeff[1][1] = 0;
+
+    // the stats structure determines how the clipping statistic is measured
+    // too few points to use the robust analysis method
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipSigma = rej;
+    stats->clipIter = nIter;
+
+    if (!psVectorClipFitPolynomial2D(line, stats, mask, 0xff, counts, cntError, x, y)) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit shutter correction.\n");
+        psFree(stats);
+        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];
+    corr->num = stats->clippedNvalues;
+    corr->stdev = stats->clippedStdev;
+
+    psFree(stats);
+    psFree(x);
+    psFree(y);
+    psFree(line);
+
+    return corr;
+}
+
+static psF32 pmShutterCorrectionModel(psVector *deriv, const psVector *params, const psVector *x)
+{
+    // 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
+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);
 
@@ -351 +344 @@
     psVector *resid = psVectorAlloc (exptime->n, PS_TYPE_F32);
     for (int i = 0; i < exptime->n; i++) {
-float divTest = exptime->data.F32[i] + corr->offref;
-if(divTest == 0) {printf("division by 0\n");}
-
         float fitCounts = corr->scale * (exptime->data.F32[i] + corr->offset) / (exptime->data.F32[i] + corr->offref);
         resid->data.F32[i] = counts->data.F32[i] - fitCounts;
@@ -373 +363 @@
     if (rawStats->sampleStdev / resStats->sampleStdev < 1.5) corr->valid = false;
     if (isnan(rawStats->sampleStdev) || isnan(resStats->sampleStdev)) corr->valid = false;
-printf("corr scale, offset, offref: %f, %f, %f\n", corr->scale, corr->offset, corr->offref);
 
     psFree (rawStats);
@@ -994 +983 @@
 float pmShutterCorrectionReference(pmShutterCorrectionData *data)
 {
-//printf("pmShutterCorrectionReference: count %d\n", corrRefCount++);
     PS_ASSERT_PTR_NON_NULL(data, NAN);
     PS_ASSERT_INT_POSITIVE(data->num, NAN);
@@ -1025 +1013 @@
         pmShutterCorrection *guess = pmShutterCorrectionGuess(newtimes, newcounts); // Guess at correction
         psTrace("psModules.detrend", 5, "Shutter correction guess: scale: %f, offset: %f, offref: %f\n", guess->scale, guess->offset, guess->offref);
-//printf("Reference guess scale, offset, offref: %f, %f, %f\n", guess->scale, guess->offset, guess->offref);
 
         pmShutterCorrection *corr = pmShutterCorrectionFullFit(newtimes, newcounts, newerrors, guess); // The actual correction
-        psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref);
-
-//printf("Reference corr scale, offset, offref: %f, %f, %f\n", corr->scale, corr->offset, corr->offref);
-        if (corr && isfinite(corr->offref) && corr->valid) {
-            psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref);
-            meanRef += corr->offref;
-            numGood++;
-        }
+
+        if (corr) {
+            psTrace("psModules.detrend", 5, "Shutter correction fit: scale: %f, offset: %f, offref: %f\n", corr->scale, corr->offset, corr->offref);
+            if (isfinite(corr->offref) && corr->valid) {
+                psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref);
+                meanRef += corr->offref;
+                numGood++;
+            }
+        } else {
+            psTrace("psModules.detrend", 5, "failed Shutter correction fit\n");
+        }
 
         psFree(corr);
@@ -1146 +1136 @@
     psImage *patternImage = pattern->image; // Illumination pattern
 
-    int num = data->num;                // NmaxInputRowsmaxInputRowsumber of images
+    int num = data->num;                // Number of images
     psVector *counts = psVectorAlloc(num, PS_TYPE_F32); // Counts in each image
     psVector *errors = psVectorAlloc(num, PS_TYPE_F32); // Counts in each image