Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c	(revision 30304)
+++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c	(revision 30305)
@@ -18,5 +18,5 @@
 
 //# define TESTING                         // TESTING output for debugging; may not work with threads!
-//# define USE_WEIGHT                      // Include weight (1/variance) in equation?
+# define USE_WEIGHT                      // Include weight (1/variance) in equation?
 # define USE_WINDOW                      // window to avoid neighbor contamination
 
@@ -1589,11 +1589,41 @@
     float momentValue = stats->sampleMean;
 
-    // XXX this is probably wrong : I want to score against higher moments **compared with the raw momements**
-    // score : chisq * dMoments * modeFactor
-    float modeFactor = kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 2.0 : 1.0;
-    float orderFactor = 0.5*(kernels->spatialOrder + 1) * (kernels->spatialOrder + 2);
-    float score = chisqValue * momentValue * modeFactor * orderFactor;
-
-    fprintf (stderr, "chisq: %f, moment: %f, score: %f\n", chisqValue, momentValue, score);
+    double sumKernel1 = 0.0, sumKernel2 = 0.0; // Sum of the kernel
+
+    // calculate the variance contribution from this smoothing kernel
+    psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, false);
+    for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) {
+        for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) {
+            if (!isfinite(modelKernel->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+                return NULL;
+            }
+            sumKernel1 += PS_SQR(modelKernel->kernel[y][x]);
+        }
+    }
+    psFree (modelKernel);
+
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+	psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, true);
+	for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) {
+	    for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) {
+		if (!isfinite(modelKernel->kernel[y][x])) {
+		    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+		    return NULL;
+		}
+		sumKernel2 += PS_SQR(modelKernel->kernel[y][x]);
+	    }
+	}
+	psFree (modelKernel);
+    } else {
+	sumKernel2 = 1.0;
+    }
+
+    // if we modify the chisq value by the (sumKernel1 + sumKernel2), we account for the
+    // smoothing coming from larger kernels adding additional spatial fit terms should be
+    // penalized by increasing the score somewhat.  the 0.01 value is not well-chosen.
+    float orderFactor = 0.01 * kernels->spatialOrder;
+    float score = 2.0 * chisqValue / (sumKernel1 + sumKernel2) + orderFactor;
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqValue, momentValue, sumKernel1, sumKernel2, score);
 
     // save this result if it is the first or the best (skip if bestMatch is NULL)
@@ -1614,5 +1644,5 @@
 	}
 	if (keep) {
-	    fprintf (stderr, "keeping order: %d, mode: %d, score: %f\n", kernels->spatialOrder, kernels->mode, score);
+	    psLogMsg("psModules.imcombine", PS_LOG_INFO, "keeping order: %d, mode: %d, score: %f\n", kernels->spatialOrder, kernels->mode, score);
 	    match->score        = score;
 	    match->spatialOrder = kernels->spatialOrder;