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

Timestamp:

Jun 18, 2008, 3:59:59 PM (18 years ago)

Author:

Paul Price

Message:

Incorporating "variance factors" --- accounting for pixel-scale variance rather than smoothed.

File:

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

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-05-14 22:22:28 $
+ *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-06-19 01:59:59 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -216 +216 @@
                           const psArray *inputs, // Stack data
                           const psVector *weights, // Global (single value) weights for data, or NULL
+                          const psVector *varFactors, // Variance factors for data
                           const psVector *reject, // Indices of pixels to reject, or NULL
                           int x, int y, // Coordinates of interest; frame of output image
@@ -233 +234 @@
     assert(numIter >= 0);
     assert(buffer);
+    assert(varFactors);
     assert((useVariance && variance) || !useVariance);
 
@@ -339 +341 @@
               psVectorInit(pixelMasks, 0);
               if (useVariance) {
-                  // Convert to rejection limits
+                  // Convert to rejection limits --- saves doing it later.
+                  // Using squared rejection limit because it's cheaper than sqrts
+                  rej *= rej;
                   for (int i = 0; i < num; i++) {
-                      pixelVariances->data.F32[i] = rej * sqrtf(pixelVariances->data.F32[i]);
+                      pixelVariances->data.F32[i] = rej * pixelVariances->data.F32[i] *
+                          varFactors->data.F32[i];
                   }
               }
@@ -364 +369 @@
               }
 
+// Mask a pixel for inspection
+#define MASK_PIXEL_FOR_INSPECTION() \
+    pixelMasks->data.PS_TYPE_MASK_DATA[j] = 0xff; \
+    combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
+    numClipped++; \
+    totalClipped++;
+
               for (int j = 0; j < num; j++) {
                   if (pixelMasks->data.PS_TYPE_MASK_DATA[j]) {
@@ -369 +381 @@
                   }
                   float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected
-                  if (diff > (useVariance ? pixelVariances->data.F32[i] : limit)) {
-                      pixelMasks->data.PS_TYPE_MASK_DATA[j] = 0xff;
-                      combineInspect(inputs, x, y, pixelSources->data.U16[j]);
-                      numClipped++;
-                      totalClipped++;
+                  if (useVariance) {
+                      // Comparing squares --- cheaper than lots of sqrts
+                      if (PS_SQR(diff) > pixelVariances->data.F32[i]) {
+                          MASK_PIXEL_FOR_INSPECTION();
+                      }
+                  } else if (diff > limit) {
+                      MASK_PIXEL_FOR_INSPECTION();
                   }
               }
@@ -561 +575 @@
     }
 
+    psVector *varFactors = psVectorAlloc(num, PS_TYPE_F32); // Variance factors for each image
     psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image
     psArray *stack = psArrayAlloc(num); // Stack of readouts
@@ -571 +586 @@
         weights->data.F32[i] = data->weight;
         stack->data[i] = psMemIncrRefCounter(data->readout);
+        // Variance factor
+        float vf = psMetadataLookupF32(NULL, data->readout->parent->concepts, "CELL.VARFACTOR"); // Var factor
+        if (!isfinite(vf)) {
+            psWarning("Non-finite CELL.VARFACTOR for image %d --- setting to unity.", i);
+            vf = 1.0;
+        }
+        varFactors->data.F32[i] = vf;
     }
 
@@ -618 +640 @@
                 for (int x = minInputCols; x < maxInputCols; x++) {
                     psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
-                    combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, reject, x, y,
-                                  maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                    combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
+                                  reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
                 }
             }
@@ -631 +653 @@
                 }
                 psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, reject, x, y,
-                              maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
+                              reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
             }
         }
@@ -658 +680 @@
         for (int y = minInputRows; y < maxInputRows; y++) {
             for (int x = minInputCols; x < maxInputCols; x++) {
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, NULL, x, y,
-                              maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
+                              NULL, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
             }
         }