Changeset 20622

Timestamp:

Nov 9, 2008, 8:32:10 PM (18 years ago)

Author:

Paul Price

Message:

Another different algorithm: sum over the column in stages and compare with the background stdev.

Location:

trunk/psModules/src/detrend

Files:

• pmRemnance.c (modified) (2 diffs)
• pmRemnance.h (modified) (1 diff)

trunk/psModules/src/detrend/pmRemnance.c

@@ -10 +10 @@
 
 #define SIZE 30                         // Size of accumulation patch
-
-#define THRESHOLD 4.0                   // Threshold above background
-
+#define THRESHOLD 20.0                   // Threshold above background
 
 bool pmRemnance(pmReadout *ro,           ///< Readout with input image
                 psMaskType maskVal,      ///< Value of mask
-                psMaskType maskRem       ///< Value to give remance
+                psMaskType maskRem,       ///< Value to give remance
+                int size,               ///< Size of accumulation patches
+                float threshold         ///< Threshold for masking
     )
 {
@@ -36 +36 @@
     }
     psFree(rng);
-    float bg = stats->robustMedian;     // Background level
+    float bgMean = stats->robustMedian; // Background level
+    float bgStdev = stats->robustStdev; // Background stdev
 
-    // Starting at the bottom of the detector, mask pixels that are significant
-    psVector *colSums = psVectorAlloc(numCols, PS_TYPE_F32); // Sums for each column
-    psVectorInit(colSums, 0.0);
-    psVector *colNums = psVectorAlloc(numCols, PS_TYPE_S32); // Number for each column
-    psVectorInit(colNums, 0);
-    psVector *colAvgs = psVectorAlloc(numCols, PS_TYPE_F32); // Average for each column
-    psVector *colMask = psVectorAlloc(numCols, PS_TYPE_MASK); // Mask for each column
+    stats->options = PS_STAT_SAMPLE_MEDIAN;
+
     int numMasked = 0;                  // Number of pixels masked
-    int numAccumulations = ceil(numRows / (float)SIZE); // Number of accumulations up the columns
-    for (int i = 0; i < numAccumulations; i++) {
-        int min = i * SIZE;             // Minimum coordinate
-        int max = PS_MIN(min + SIZE, numRows); // Maximum coordinate
-        for (int x = 0; x < numCols; x++) {
-            int num = 0;                // Number of pixels in accumulation
-            float sum = 0.0;            // Sum of pixels in accumulation
+    int number = ceil(numRows / (float)SIZE); // Number of steps up the columns
+    psVector *values = psVectorAlloc(numRows, PS_TYPE_F32); // Values below center
+    for (int x = 0; x < numCols; x++) {
+        int maxMask = 0;                // Maximum row to which to mask
+        int numValues = 0;              // Number of values
+        for (int i = 0, min = 0, max = size; i < number; i++, min += size, max += size) {
+
+            if (max > numRows) {
+                max = numRows;
+            }
             for (int y = min; y < max; y++) {
                 if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
                     continue;
                 }
-                sum += image->data.F32[y][x] - bg;
-                num++;
+                values->data.F32[numValues++] = image->data.F32[y][x];
             }
-            colSums->data.F32[x] = sum;
-            colNums->data.S32[x] = num;
-            colAvgs->data.F32[x] = sum / (float)num;
-            colMask->data.PS_TYPE_MASK_DATA[x] = (colNums->data.S32[x] == 0 ? 0xFF : 0);
+            values->n = numValues;
+            if (!psVectorStats(stats, values, NULL, NULL, 0)) {
+                // Can't do anything about it
+                psErrorClear();
+                maxMask = max;
+                continue;
+            }
+            float median = stats->sampleMedian;
+
+            if (median > bgMean + threshold * bgStdev / sqrtf(numValues)) {
+                maxMask = max;
+            }
         }
 
-        if (!psVectorStats(stats, colAvgs, NULL, colMask, 0xFF)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on column accumulations.");
-            psFree(colSums);
-            psFree(colNums);
-            psFree(colAvgs);
-            psFree(colMask);
-            psFree(stats);
-            return false;
+        maxMask += size;
+        if (maxMask > numRows) {
+            maxMask = numRows;
         }
+        for (int y = 0; y < maxMask; y++) {
+            mask->data.PS_TYPE_MASK_DATA[y][x] |= maskRem;
+        }
+        numMasked += maxMask;
 
-        float threshold = stats->robustMedian + THRESHOLD * stats->robustStdev; // Threshold for masking
-
-        max = PS_MIN(max + SIZE, numRows);
-
-        for (int x = 0; x < numCols; x++) {
-            if (colAvgs->data.F32[x] > threshold) {
-                for (int y = 0; y < max; y++) {
-                    mask->data.PS_TYPE_MASK_DATA[y][x] = maskRem;
-                    numMasked++;
-                }
-            }
-        }
     }
-
-    psFree(colSums);
-    psFree(colNums);
-    psFree(colAvgs);
-    psFree(colMask);
+    psFree(values);
     psFree(stats);
 

trunk/psModules/src/detrend/pmRemnance.h

@@ -14 +14 @@
 bool pmRemnance(pmReadout *ro,           ///< Readout with input image
                 psMaskType maskVal,      ///< Value of mask
-                psMaskType maskRem       ///< Value to give remance
+                psMaskType maskRem,       ///< Value to give remance
+                int size,               ///< Size of accumulation patches
+                float threshold         ///< Threshold for masking
     );
 
-
 #endif