Changeset 18350 for branches/eam_branch_20080620/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Jun 27, 2008, 1:56:41 PM (18 years ago)

Author:

eugene

Message:

accepting updates from HEAD

File:

• branches/eam_branch_20080620/psModules/src/imcombine/pmSubtraction.c (modified) (5 diffs)

branches/eam_branch_20080620/psModules/src/imcombine/pmSubtraction.c

@@ -3 +3 @@
  *  @author Paul Price, IfA
  *  @author GLG, MHPCC
- *
- *  @version $Revision: 1.95 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-06-19 03:03:53 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -512 +509 @@
     PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1);
 
-    double totalSquareDev = 0.0;        // Total square deviation from zero
+    // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but
+    // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like
+    // with large N.  Therefore, let's just treat this as a Gaussian distribution.
+
+    double mean = 0.0;                  // Mean deviation
     int numStamps = 0;                  // Number of used stamps
     for (int i = 0; i < stamps->num; i++) {
@@ -519 +520 @@
             continue;
         }
-        totalSquareDev += PS_SQR(deviations->data.F32[i]);
+        mean += deviations->data.F32[i];
         numStamps++;
     }
-
-    float rms = sqrt(totalSquareDev / (double)numStamps); // Convert to RMS
+    mean /= numStamps;
+
+    double rms = 0.0;                   // Standard deviation
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            continue;
+        }
+        rms += PS_SQR(deviations->data.F32[i] - mean);
+    }
+    rms = sqrt(rms / (numStamps - 1));
 
     if (rmsPtr) {
@@ -549 +559 @@
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
-    double newSquareDev = 0.0;          // New square deviation
+    double newMean = 0.0;               // New mean
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
         if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
-            if (deviations->data.F32[i] > limit) {
+            // Should we reject stars with low deviation?  Well, if this is really a Gaussian-like
+            // distribution and they're low, then we have the right to ask why.  Isn't it suspicious that
+            // they're anomalously low, compared to the rest of the population which (we hope) is indicative
+            // of normality?  Besides, the standard deviation is going to be blown up by stars that didn't
+            // subtract well, in which case very few (if any) stars will be legitimately rejected for being
+            // low.
+            if (deviations->data.F32[i] - mean > limit) {
                 // Mask out the stamp in the image so you it's not found again
                 psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
@@ -587 +603 @@
             } else {
                 numGood++;
-                newSquareDev += PS_SQR(deviations->data.F32[i]);
+                newMean += deviations->data.F32[i];
             }
         }
     }
+    newMean /= numGood;
 
     if (numRejected > 0) {
         psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n",
-                 numGood, numRejected, rms,
-                 sqrt(newSquareDev / (double)numGood));
+                 "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n",
+                 numGood, numRejected, mean, newMean);
     } else {
         psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; 0 rejected.\nRMS deviation: %f\n",
-                 numGood, rms);
+                 "%d good stamps; 0 rejected.\nMean deviation: %lf\n",
+                 numGood, mean);
     }