Changeset 15399 for branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Oct 28, 2007, 3:40:54 PM (19 years ago)

Author:

eugene

Message:

merging diffs from cvs head

File:

• branches/eam_branch_20071023/psModules/src/imcombine/pmSubtraction.c (modified) (4 diffs)

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

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.65 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-17 02:45:40 $
+ *  @version $Revision: 1.65.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-29 01:40:54 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -907 +907 @@
                     for (int j = 0; j < numKernels; j++) {
                         psKernel *convolution = convolutions->data[j]; // Convolution of reference
+
+                        // XXX Precalculate these values, so that we're not calculating the same thing for
+                        // every pixel.
                         double polynomial = 0.0; // Value of the polynomial
                         for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) {
@@ -928 +931 @@
             deviations->data.F32[i] = devNorm * deviation;
             psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
-                    i, (int)stamp->x, (int)stamp->y, deviations->data.F32[i]);
+                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
             totalSquareDev += PS_SQR(deviations->data.F32[i]);
             numStamps++;
@@ -947 +950 @@
     }
 
+    if (!isfinite(sigmaRej) || sigmaRej <= 0.0) {
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "RMS deviation: %f", sqrt(totalSquareDev / numStamps));
+        psFree(deviations);
+        return 0;
+    }
+
+    int numRejected = 0;                // Number of stamps rejected
+    int numGood = 0;                    // Number of good stamps
+    double newSquareDev = 0.0;          // New square deviation
+
     float limit = sigmaRej * sqrt(totalSquareDev / numStamps); // Limit on maximum deviation
     psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
 
-    int numRejected = 0;
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > 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,
-                    (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
-            numRejected++;
-            for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
-                for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
-                    subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
+        if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
+            if (deviations->data.F32[i] > 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,
+                        (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+                numRejected++;
+                for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
+                    for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
+                        subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
+                    }
                 }
-            }
-
-            // Set stamp for replacement
-            stamp->x = 0;
-            stamp->y = 0;
-            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            // Recalculate convolutions
-            psFree(stamp->convolutions);
-            stamp->convolutions = NULL;
-        }
-    }
+
+                // Set stamp for replacement
+                stamp->x = 0;
+                stamp->y = 0;
+                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+                // Recalculate convolutions
+                psFree(stamp->convolutions);
+                stamp->convolutions = NULL;
+            } else {
+                numGood++;
+                newSquareDev += PS_SQR(deviations->data.F32[i]);
+            }
+        }
+    }
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n",
+             numGood, numRejected, sqrt(totalSquareDev / numStamps), sqrt(newSquareDev / numGood));
 
     psFree(deviations);