Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 13369)
+++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 13370)
@@ -25,4 +25,7 @@
     )
 {
+    if (exp == 0) {
+        return 1.0;
+    }
     float result = value;               // Result to return
     for (int i = 2; i < exp; i++) {
@@ -82,5 +85,5 @@
     }
 
-    kernels->subIndex = PS_SQR(size + 1);
+    kernels->subIndex = num / 2;
     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
            kernels->v->data.S32[kernels->subIndex] == 0 &&
@@ -106,5 +109,6 @@
     int num = 0;                        // Number of basis functions
     for (int i = 0; i < numGaussians; i++) {
-        num += sigmas->data.F32[i] * orders->data.F32[i];
+        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
+        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
     }
     num *= (spatialOrder + 1) * (spatialOrder + 2) / 2;
@@ -124,7 +128,7 @@
                 psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
                 for (int v = -size; v <= size; v++) {
-                    for (int u = -size; u < size; u++) {
+                    for (int u = -size; u <= size; u++) {
                         preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
-                            expf(0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
+                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
                     }
                 }
@@ -147,4 +151,10 @@
     }
 
+    kernels->subIndex = (spatialOrder + 1) * (spatialOrder + 2) / 4;
+    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
+           kernels->v->data.S32[kernels->subIndex] == 0 &&
+           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
+           kernels->yOrder->data.S32[kernels->subIndex] == 0);
+
     return kernels;
 }