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

Timestamp:

Aug 8, 2008, 8:09:07 AM (18 years ago)

Author:

Paul Price

Message:

Fixing to match slight differences in psThread

File:

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

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -122 +122 @@
 
     // Iterate only
-    for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++); 
+    for (index = 0; !isfinite(exptime->data.F32[index]) && index < N - 1; index++);
 
     if (index == N - 1) {
@@ -344 +344 @@
     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;
-    }
-    
+        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);
@@ -353 +353 @@
     psVectorStats (resStats, resid, NULL, NULL, 0);
 
-    // XXX temporary hard-wired minimum stdev improvement factor 
+    // 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;
@@ -433 +433 @@
             numRows = image->numRows;
             numCols = image->numCols;
-            // define the reference region : a box of size 'size' at the center
+            // 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
@@ -656 +656 @@
 
 
-bool pmShutterCorrectionApplyScan_Threaded (psThreadJob *job) {
+bool pmShutterCorrectionApplyScan_Threaded(const psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
     psImage *image        = job->args->data[0];
     psImage *shutterImage = job->args->data[1];
@@ -672 +675 @@
 {
     for (int y = rowStart; y < rowStop; y++) {
-        for (int x = 0; x < image->numCols; x++) {
-            if (mask && !isfinite(shutterImage->data.F32[y][x])) {
-                mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
-                image->data.F32[y][x] = NAN;
-                continue;
-            }
-            image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]);
-        }
+        for (int x = 0; x < image->numCols; x++) {
+            if (mask && !isfinite(shutterImage->data.F32[y][x])) {
+                mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
+                image->data.F32[y][x] = NAN;
+                continue;
+            }
+            image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]);
+        }
     }
     return true;
@@ -732 +735 @@
     bool threaded = true;
     int scanRows = pmDetrendGetScanRows();
-    if (scanRows == 0) { 
-        threaded = false;
-        scanRows = image->numRows;
+    if (scanRows == 0) {
+        threaded = false;
+        scanRows = image->numRows;
     }
 
@@ -756 +759 @@
         psFree (line);
     } else {
-        for (int rowStart = 0; rowStart < image->numRows; rowStart += scanRows) {
-            int rowStop = PS_MIN (rowStart + scanRows, image->numRows);
-
-# define PS_ARRAY_ADD_SCALAR(ARRAY, VALUE, TYPE) { \
-      psScalar *scalar = psScalarAlloc(VALUE, TYPE); \
-      psArrayAdd(ARRAY, 1, scalar); \
-      psFree (scalar); }
-      
-            if (threaded) {
-                // allocate a job, construct the arguments for this job
-                psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_SHUTTER");
-                psArrayAdd (job->args, 1, image);
-                psArrayAdd (job->args, 1, shutterImage);
-                psArrayAdd (job->args, 1, mask);
-                PS_ARRAY_ADD_SCALAR (job->args, exptime, PS_TYPE_F32);
-                PS_ARRAY_ADD_SCALAR (job->args, blank, PS_TYPE_MASK);
-                PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
-                PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
-
-                // ppImageDetrendReadout(config, options, view)
-                if (!psThreadJobAddPending (job)) {
-                    return false;
-                }
-            } else {
-                pmShutterCorrectionApplyScan (image, shutterImage, mask, exptime, blank, rowStart, rowStop);
-            }
-        }
-        if (threaded) {
-            // wait here for the threaded jobs to finish
-            if (!psThreadPoolWait ()) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
-                return false;
-            }
-            fprintf (stderr, "success for threaded jobs\n");
-
-            // free the done jobs
-            psThreadJob *job = NULL;
-            while ((job = psThreadJobGetDone()) != NULL) {
-                psFree (job);
-            }
-        }
+        for (int rowStart = 0; rowStart < image->numRows; rowStart += scanRows) {
+            int rowStop = PS_MIN (rowStart + scanRows, image->numRows);
+
+            if (threaded) {
+                // allocate a job, construct the arguments for this job
+                psThreadJob *job = psThreadJobAlloc ("PSMODULES_DETREND_SHUTTER");
+                psArrayAdd (job->args, 1, image);
+                psArrayAdd (job->args, 1, shutterImage);
+                psArrayAdd (job->args, 1, mask);
+                PS_ARRAY_ADD_SCALAR (job->args, exptime, PS_TYPE_F32);
+                PS_ARRAY_ADD_SCALAR (job->args, blank, PS_TYPE_MASK);
+                PS_ARRAY_ADD_SCALAR (job->args, rowStart, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR (job->args, rowStop, PS_TYPE_S32);
+
+                // ppImageDetrendReadout(config, options, view)
+                if (!psThreadJobAddPending (job)) {
+                    psFree(job);
+                    return false;
+                }
+                psFree(job);
+            } else {
+                pmShutterCorrectionApplyScan (image, shutterImage, mask, exptime, blank, rowStart, rowStop);
+            }
+        }
+        if (threaded) {
+            // wait here for the threaded jobs to finish
+            if (!psThreadPoolWait(false)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
+                return false;
+            }
+            fprintf (stderr, "success for threaded jobs\n");
+
+            // free the done jobs
+            psThreadJob *job = NULL;
+            while ((job = psThreadJobGetDone()) != NULL) {
+                psFree (job);
+            }
+        }
     }
     psFree(shutterImage);
@@ -952 +952 @@
         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]);
+        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);
@@ -978 +978 @@
     // 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]];
+        newtimes->data.F32[j] = data->exptimes->data.F32[index->data.S32[j]];
     }
 
@@ -984 +984 @@
     int numGood = 0;                    // Number of good measurements
     for (int i = 0; i < MEASURE_SAMPLES; i++) {
-        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
+        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);
+        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);
+        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) {
@@ -1009 +1009 @@
         psFree(corr);
         psFree(guess);
-        psFree (newcounts);
-        psFree (newerrors);
+        psFree (newcounts);
+        psFree (newerrors);
     }
     psFree (newtimes);
@@ -1040 +1040 @@
     pmReadoutStackDefineOutput(shutter, col0, row0, numCols, numRows, false, false, maskVal);
     if (pattern) {
-        pmReadoutStackDefineOutput(pattern, col0, row0, numCols, numRows, false, false, maskVal);
+        pmReadoutStackDefineOutput(pattern, col0, row0, numCols, numRows, false, false, maskVal);
     }
 
@@ -1072 +1072 @@
 
     pattern->data_exists = true;
-    if (pattern->parent) { 
-        pattern->parent->data_exists = true;
-        if (pattern->parent->parent) {
-            pattern->parent->parent->data_exists = true;
-        }
+    if (pattern->parent) {
+        pattern->parent->data_exists = true;
+        if (pattern->parent->parent) {
+            pattern->parent->parent->data_exists = true;
+        }
     }
 
@@ -1137 +1137 @@
                     errors->data.F32[r] = sqrtf(readout->weight->data.F32[yIn][xIn]) * ref;
                 } else {
-                    // XXX guess that the input data is Poisson distributed; if we go negative, force high
+                    // XXX guess that the input data is Poisson distributed; if we go negative, force high
                     errors->data.F32[r] = sqrtf(fabs(image->data.F32[yIn][xIn])) * ref;
                 }