Changeset 20497 for trunk/psModules/src/imcombine/pmStack.c

Timestamp:

Oct 31, 2008, 5:00:04 PM (18 years ago)

Author:

Paul Price

Message:

Adding relative systematic error for combination when rejecting.

File:

• trunk/psModules/src/imcombine/pmStack.c (modified) (8 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-31 21:50:59 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-11-01 02:59:33 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -219 +219 @@
                           int numIter, // Number of rejection iterations
                           float rej, // Number of standard deviations at which to reject
+                          float sys,    // Relative systematic error
                           bool useVariance, // Use variance for rejection when combining?
                           bool safe,    // Combine safely?
@@ -311 +312 @@
           if (useVariance && numIter > 0) {
               // Use variance to check that the two are consistent
-              float diff = pixelData->data.F32[0] - pixelData->data.F32[1];
+              float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
+              float var1 = pixelVariances->data.F32[0]; // Variance of first
+              float var2 = pixelVariances->data.F32[1]; // Variance of second
 #ifdef VARIANCE_FACTORS
-              float sigma2 = pixelVariances->data.F32[0] * varFactors->data.F32[pixelSources->data.U16[0]] +
-                  pixelVariances->data.F32[1] * varFactors->data.F32[pixelSources->data.U16[1]];
-#else
-              float sigma2 = pixelVariances->data.F32[0] + pixelVariances->data.F32[1];
+              // Variance factor contributes to the rejection level
+              var1 *= varFactors->data.F32[pixelSources->data.U16[0]];
+              var2 *= varFactors->data.F32[pixelSources->data.U16[1]];
 #endif
+              // Systematic error contributes to the rejection level
+              var1 += PS_SQR(sys * pixelData->data.F32[0]);
+              var2 += PS_SQR(sys * pixelData->data.F32[1]);
+
+              float sigma2 = var1 + var2; //
               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
                   // Not consistent: mark both for inspection
@@ -346 +353 @@
                   float rej2 = PS_SQR(rej); // Rejection level squared
                   for (int i = 0; i < num; i++) {
+                      pixelVariances->data.F32[i] *= rej2;
 #ifdef VARIANCE_FACTORS
-                      pixelVariances->data.F32[i] *= rej2 * varFactors->data.F32[pixelSources->data.U16[i]];
-#else
-                      pixelVariances->data.F32[i] *= rej2;
+                      // Variance factor contributes to the rejection level
+                      pixelVariances->data.F32[i] *= varFactors->data.F32[pixelSources->data.U16[i]];
 #endif
+                      // Systematic error contributes to the rejection level
+                      pixelVariances->data.F32[i] += PS_SQR(sys * pixelData->data.F32[i]);
                   }
               }
@@ -554 +563 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined, psArray *input, psMaskType maskVal, psMaskType bad,
-                    int kernelSize, int numIter, float rej, bool entire, bool useVariance, bool safe)
+                    int kernelSize, int numIter, float rej, float sys,
+                    bool entire, bool useVariance, bool safe)
 {
     PS_ASSERT_PTR_NON_NULL(combined, false);
@@ -660 +670 @@
                                                      x, y); // Reject these images
                     combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                                  reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                                  reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
                 }
             }
@@ -674 +684 @@
                                                  x, y); // Reject these images
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                              reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                              reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
             }
         }
@@ -701 +711 @@
             for (int x = minInputCols; x < maxInputCols; x++) {
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                              NULL, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                              NULL, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
             }
         }