Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c	(revision 26552)
+++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c	(revision 26553)
@@ -154,5 +154,5 @@
         for (int u = -size; u <= size; u++) {
             double value = preCalc->kernel->kernel[v][u];
-            moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
+            moment += PS_SQR(value) * PS_SQR((PS_SQR(u) + PS_SQR(v)));
         }
     }
@@ -181,7 +181,5 @@
     // only even terms have non-zero sums
     if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) {
-        moment /= sum;
-    } else {
-        moment = 0.0;
+        moment /= PS_SQR(sum);
     }
 
@@ -231,7 +229,4 @@
     kernels->preCalc->data[index] = preCalc;
     kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
-    if (!isfinite(kernels->penalties->data.F32[index])) {
-        psAbort ("invalid penalty");
-    }
 
     psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
@@ -893,5 +888,5 @@
                                     preCalc->poly->data.F32[j] = polyVal;
                                     norm += polyVal;
-                                    moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v));
+                                    moment += PS_SQR(polyVal) * PS_SQR(PS_SQR(u) + PS_SQR(v));
 
                                     psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1);
@@ -919,5 +914,5 @@
                         psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
                     }
-                    moment /= norm;
+                    moment /= PS_SQR(norm);
                 }