Changeset 5717 for trunk/pois/src/poisRejectStamps.c

Timestamp:

Dec 6, 2005, 6:43:52 PM (21 years ago)

Author:

Paul Price

Message:

Importing current working PAP version (many bug fixes since this branch) straight on top of this version, and tidying up a bit

File:

• trunk/pois/src/poisRejectStamps.c (modified) (2 diffs)

trunk/pois/src/poisRejectStamps.c

@@ -6 +6 @@
 #define SQUARE(x) ((x)*(x))
 
-psList *restrict poisRejectStamps(psArray *stamps, // Array of stamps
-                                  psImage *mask,        // Mask image
-                                  const psVector *deviations, // Vector of deviations for the stamps
-                                  const poisConfig *config // Configuration
+bool poisRejectStamps(psArray *stamps,  // Array of stamps
+                      psImage *mask,    // Mask image
+                      const psVector *deviations, // Vector of deviations for the stamps
+                      const poisConfig *config // Configuration
     )
 {
@@ -21 +21 @@
     assert(deviations->type.type == PS_TYPE_F32);
 
-    psVector *devMask = psVectorAlloc(stamps->n, PS_TYPE_U8); // Mask for statistics
+    bool masked = false;                // Have we masked any stamps?
 
     // Don't want the standard deviation, but rather the deviation from zero
-    double meanDev = 0.0;
+    double sum = 0.0;
     int numDev = 0;
     for (int i = 0; i < deviations->n; i++) {
-        poisStamp *stamp = stamps->data[i];     // The stamp
-        // Only interested in the stamps that we're actually using
-        if (stamp != NULL && stamp->status == POIS_STAMP_USED) {
-            meanDev += SQUARE(deviations->data.F32[i]);
-            numDev++;
-        }
+        poisStamp *stamp = stamps->data[i];     // The stamp
+        // Only interested in the stamps that we're actually using
+        if (stamp->status == POIS_STAMP_USED) {
+            sum += SQUARE(deviations->data.F32[i]);
+            numDev++;
+        }
     }
-    float rmsDev = sqrt(meanDev/numDev);
-    float limit = rmsDev * config->sigmaRej;
-    psTrace("pois.rejectStamps", 2, "RMS deviation: %f\n", rmsDev);
+    float meanDev = sqrtf(sum / (float)numDev);
+    float limit = meanDev * config->sigmaRej;
+    psTrace("pois.rejectStamps", 2, "Mean RMS deviation: %f\n", meanDev);
     psTrace("pois.rejectStamps", 2, "Rejection limit: %f\n", limit);
 
     // Reject stamps
-    psList *rejected = NULL;            // rejected stamps
     for (int s = 0; s < deviations->n; s++) {
-        poisStamp *stamp = stamps->data[s];
-        if (stamp == NULL) {
-            continue;
-        }
-        
-        if (stamp->status == POIS_STAMP_USED && fabsf(deviations->data.F32[s]) > limit) {
-            psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y,
-                    deviations->data.F32[s]);
-            
-            // Mask out the stamp in the image so you don't find it again
-            for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
-                for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
-                    mask->data.U8[y][x] |= POIS_MASK_STAMP;
-                }
-            }
-            
-            // Mark stamp for replacement
-            stamps->data[s] = NULL;
-            
-            if (rejected == NULL) {
-                rejected = psListAlloc(stamp);
-            } else {
-                bool success = psListAdd(rejected, PS_LIST_TAIL, stamp);
-                assert(success);
-            }
-            psFree(stamp);
-        } // Bad stamps
+        poisStamp *stamp = stamps->data[s];
+        if (stamp->status == POIS_STAMP_USED &&
+            (fabsf(deviations->data.F32[s]) > limit || ! isfinite(deviations->data.F32[s]))) {
+            masked = true;
+            psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y,
+                    deviations->data.F32[s]);
+
+            // Mask out the stamp in the image so you don't find it again
+            for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
+                for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
+                    mask->data.U8[y][x] |= POIS_MASK_STAMP;
+                }
+            }
+
+            // Set stamp for replacement
+            stamp->x = 0;
+            stamp->y = 0;
+            stamp->status = POIS_STAMP_RESET;
+
+        } // Bad stamps
     } // Iterating over stamps
 
-    psFree(devMask);
-
-    return rejected;
+    return masked;
 }