Changeset 20659

Timestamp:

Nov 10, 2008, 3:31:25 PM (17 years ago)

Author:

Paul Price

Message:

Adding rejection of images based on high match chi²

Location:

trunk/ppStack/src

Files:

• ppStack.h (modified) (1 diff)
• ppStackLoop.c (modified) (2 diffs)

trunk/ppStack/src/ppStack.h

@@ -11 +11 @@
 typedef enum {
     PPSTACK_MASK_MATCH  = 0x01,         // PSF-matching failed
-    PPSTACK_MASK_REJECT = 0x02,         // Rejection failed
-    PPSTACK_MASK_BAD    = 0x04,         // Bad image (too many pixels rejected)
+    PPSTACK_MASK_CHI2   = 0x02,         // Chi^2 too deviant
+    PPSTACK_MASK_REJECT = 0x04,         // Rejection failed
+    PPSTACK_MASK_BAD    = 0x08,         // Bad image (too many pixels rejected)
     PPSTACK_MASK_ALL    = 0xff          // All errors
 } ppStackMask;

trunk/ppStack/src/ppStackLoop.c

@@ -532 +532 @@
         psFree(cells);
         psFree(inputMask);
+        psFree(matchChi2);
+        return false;
+    }
+
+    // Reject images out-of-hand on the basis of their match chi^2
+    {
+        psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort
+        for (int i = 0; i < num; i++) {
+            if (inputMask->data.PS_TYPE_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+                continue;
+            }
+            values->data.F32[values->n++] = matchChi2->data.F32[i];
+        }
+        assert(values->n == numGood);
+        if (!psVectorSortInPlace(values)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to sort vector.");
+            psFree(subKernels);
+            psFree(subRegions);
+            psFree(cells);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(values);
+            return false;
+        }
+        float median = numGood % 2 ? values->data.F32[numGood / 2] :
+            0.5 * (values->data.F32[numGood / 2 - 1] + values->data.F32[numGood / 2]);
+
+        float rms = 0.74 * (values->data.F32[numGood * 3 / 4] -
+                            values->data.F32[numGood / 4]); // Estimated RMS from interquartile range
+
+        psFree(values);
+
+        float rej = psMetadataLookupF32(NULL, recipe, "MATCH.REJ"); // Rejection threshold (stdevs)
+        if (isfinite(rej)) {
+            float thresh = median + rej * rms; // Threshold for rejection
+            psLogMsg("ppStack", PS_LOG_INFO, "chi^2 rejection threshold = %f + %f * %f = %f",
+                     median, rej, rms, thresh);
+
+            int numRej = 0;                 // Number rejected
+            numGood = 0;                    // Number of good images
+            for (int i = 0; i < num; i++) {
+                if (inputMask->data.PS_TYPE_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+                    continue;
+                }
+                if (matchChi2->data.F32[i] > thresh) {
+                    numRej++;
+                    inputMask->data.PS_TYPE_MASK_DATA[i] |= PPSTACK_MASK_CHI2;
+                    psLogMsg("ppStack", PS_LOG_INFO, "Rejecting image %d because of large matching chi^2: %f",
+                             i, matchChi2->data.F32[i]);
+                } else {
+                    numGood++;
+                }
+            }
+        }
+    }
+
+    if (numGood == 0) {
+        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
+        psFree(subKernels);
+        psFree(subRegions);
+        psFree(cells);
         psFree(inputMask);
         psFree(matchChi2);
@@ -776 +837 @@
         // Reject bad pixels
         for (int i = 0; i < num; i++) {
+            if (inputMask->data.U8[i]) {
+                continue;
+            }
             psTimerStart("PPSTACK_REJECT");