Index: /trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- /trunk/psModules/src/imcombine/pmStack.c	(revision 18181)
+++ /trunk/psModules/src/imcombine/pmStack.c	(revision 18182)
@@ -8,6 +8,6 @@
  *  @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,4 +216,5 @@
                           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,4 +234,5 @@
     assert(numIter >= 0);
     assert(buffer);
+    assert(varFactors);
     assert((useVariance && variance) || !useVariance);
 
@@ -339,7 +341,10 @@
               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,4 +369,11 @@
               }
 
+// 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,9 +381,11 @@
                   }
                   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,4 +575,5 @@
     }
 
+    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,4 +586,11 @@
         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,6 +640,6 @@
                 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,6 +653,6 @@
                 }
                 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,6 +680,6 @@
         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);
             }
         }