Changeset 20621 for trunk/psModules/src/detrend/pmRemnance.c

Timestamp:

Nov 9, 2008, 7:17:43 PM (18 years ago)

Author:

Paul Price

Message:

Changing from drop-dead simple remnance algorithm to something a bit
more complicated: accumulate up the columns and compare each column
with the others.

File:

• trunk/psModules/src/detrend/pmRemnance.c (modified) (2 diffs)

trunk/psModules/src/detrend/pmRemnance.c

@@ -8 +8 @@
 #include "pmFPA.h"
 #include "pmRemnance.h"
+
+#define SIZE 30                         // Size of accumulation patch
 
 #define THRESHOLD 4.0                   // Threshold above background
@@ -34 +36 @@
     }
     psFree(rng);
-    float thresh = stats->robustMedian + THRESHOLD * stats->robustStdev; // Threshold for significant pixels
-    psFree(stats);
+    float bg = stats->robustMedian;     // Background level
 
     // 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
     int numMasked = 0;                  // Number of pixels masked
-    for (int x = 0; x < numCols; x++) {
-        bool significant = true;        // Are pixels significant?
-        for (int y = 0; significant && y < numRows; y++) {
-            if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
-                continue;
+    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
+            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++;
             }
-            if (image->data.F32[y][x] > thresh) {
-                mask->data.PS_TYPE_MASK_DATA[y][x] |= maskRem;
-                numMasked++;
-            } else {
-                significant = false;
+            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);
+        }
+
+        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;
+        }
+
+        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(stats);
+
 
     psMetadataAddS32(ro->analysis, PS_LIST_TAIL, PM_REMNANCE_ANALYSIS_NUM, 0,