Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 18043)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 18163)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.93 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-06-10 02:43:01 $
+ *  @version $Revision: 1.94 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-06-17 22:16:38 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -627,4 +627,30 @@
 }
 
+
+float pmSubtractionVarianceFactor(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
+{
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NAN);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NAN);
+
+    // Precalulate polynomial values
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y);
+
+    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel
+    psFree(polyValues);
+
+    double sumKernel2 = 0.0;            // Sum of the kernel squared
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+
+    psFree(kernel);
+
+    return sumKernel2;
+}
+
 #if 0
 psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)