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

Timestamp:

Jun 8, 2008, 2:42:04 PM (18 years ago)

Author:

eugene

Message:

add some trace entries; sort input exptime and counts to pmShutterCorrectionGuess (needed for valid result); add test of improved scatter for pmShutterCorectionFit

File:

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

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -80 +80 @@
     corr->num = 0;
     corr->stdev = NAN;
+    corr->valid = true;
 
     return corr;
@@ -118 +119 @@
 
     long index;
-    for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++)
-        ; // Iterate only
+
+    // 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,
@@ -336 +339 @@
     corr->offset = params->data.F32[1];
     corr->offref = params->data.F32[2];
+
+    // apply the correction and measure the residual scatter
+    psVector *resid = psVectorAlloc (exptime->n, PS_TYPE_F32);
+    for (int i = 0; i < exptime->n; i++) {
+        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;
+    }
+    
+    psStats *rawStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    psStats *resStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    psVectorStats (rawStats, counts, NULL, NULL, 0);
+    psVectorStats (resStats, resid, NULL, NULL, 0);
+
+    // XXX temporary hard-wired minimum stdev improvement factor 
+    psTrace("psModules.detrend", 3, "raw scatter %f vs res scatter %f\n", rawStats->sampleStdev, resStats->sampleStdev);
+    if (rawStats->sampleStdev / resStats->sampleStdev < 1.5) corr->valid = false;
+
+    psFree (rawStats);
+    psFree (resStats);
+    psFree (resid);
 
     psFree(minInfo);
@@ -409 +432 @@
             numRows = image->numRows;
             numCols = image->numCols;
-            // Set up the sample regions
+            // define the reference region : a box of size 'size' at the center
             refRegion = psRegionForSquare(0.5 * numCols, 0.5 * numRows, size);
+            // Set up the sample regions : boxes of size 'size' at the 4 image corners
             for (int j = 0; j < MEASURE_SAMPLES; j++) {
                 int x = (j % 2) ? size : image->numCols - size;
@@ -823 +847 @@
     psFree(rng);
     float refValue = psStatsGetValue(stats, meanStat); // Reference value
-    psTrace("psModules.detrend", 3, "Reference value for shutter image = %f\n", refValue);
+    psTrace("psModules.detrend", 3, "Reference value & exptime for shutter image : %f cnts %f sec\n", refValue, exptime);
     if (refValue <= 0.0) {
         psError(PS_ERR_UNKNOWN, true, "Measured non-positive reference value.\n");
@@ -860 +884 @@
         mean->data.F32[mean->n] = psStatsGetValue(stats, meanStat) * refValue;
         stdev->data.F32[stdev->n] = psStatsGetValue(stats, stdevStat) * refValue;
+
+        psTrace("psModules.detrend", 5, "input shutter image sample value %d: %f +/- %f  ->  %f +/- %f\n", j,
+                psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat),
+                mean->data.F32[mean->n], stdev->data.F32[stdev->n]);
+
         data->mean->data[j] = psVectorExtend(mean, IMAGES_BUFFER, 1);
         data->stdev->data[j] = psVectorExtend(stdev, IMAGES_BUFFER, 1);
@@ -869 +898 @@
 }
 
-float pmShutterCorrectionReference(const pmShutterCorrectionData *data)
+float pmShutterCorrectionReference(pmShutterCorrectionData *data)
 {
     PS_ASSERT_PTR_NON_NULL(data, NAN);
     PS_ASSERT_INT_POSITIVE(data->num, NAN);
+
+    // supply counts sorted by exptime
+
+    // generate the index for the exptimes vector
+    psVector *index = psVectorSortIndex (NULL, data->exptimes);
+    psVector *newtimes = psVectorAlloc (data->exptimes->n, PS_TYPE_F32);
+
+    // reshuffle exptimes to new sequence (this is only a local value)
+    for (int j = 0; j < newtimes->n; j++) {
+        newtimes->data.F32[j] = data->exptimes->data.F32[index->data.S32[j]];
+    }
 
     float meanRef = 0.0;                // Mean reference offset
     int numGood = 0;                    // Number of good measurements
     for (int i = 0; i < MEASURE_SAMPLES; i++) {
-        psVector *counts = data->mean->data[i]; // Mean for each image
-        psVector *errors = data->stdev->data[i]; // Stdev for each image
-        pmShutterCorrection *guess = pmShutterCorrectionGuess(data->exptimes, counts); // Guess at correction
-        pmShutterCorrection *corr = pmShutterCorrectionFullFit(data->exptimes, counts,
-                                                               errors, guess); // The actual correction
-        psFree(guess);
-        if (corr && isfinite(corr->offref)) {
+        psVector *newcounts = psVectorAlloc (data->exptimes->n, PS_TYPE_F32);
+        psVector *newerrors = psVectorAlloc (data->exptimes->n, PS_TYPE_F32);
+        psVector *counts = data->mean->data[i];
+        psVector *errors = data->stdev->data[i];
+
+        for (int j = 0; j < newcounts->n; j++) {
+            newcounts->data.F32[j] = counts->data.F32[index->data.S32[j]];
+            newerrors->data.F32[j] = errors->data.F32[index->data.S32[j]];
+        }
+
+        // use the sorted exptime, counts, and errors for the measurements
+        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);
+
+        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);
+
+        if (corr && isfinite(corr->offref) && corr->valid) {
             psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref);
             meanRef += corr->offref;
             numGood++;
         }
+
         psFree(corr);
-    }
+        psFree(guess);
+        psFree (newcounts);
+        psFree (newerrors);
+    }
+    psFree (newtimes);
+    psFree (index);
 
     if (numGood == 0) {