Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 14420)
+++ trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 14455)
@@ -19,6 +19,4 @@
     psFree(kernels->uStop);
     psFree(kernels->vStop);
-    psFree(kernels->xOrder);
-    psFree(kernels->yOrder);
     psFree(kernels->preCalc);
 }
@@ -50,4 +48,5 @@
 
     kernels->type = type;
+    kernels->num = numBasisFunctions;
     kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
     kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
@@ -55,6 +54,4 @@
     kernels->uStop = NULL;
     kernels->vStop = NULL;
-    kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
-    kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
     kernels->subIndex = 0;
     kernels->preCalc = NULL;
@@ -62,4 +59,5 @@
     kernels->inner = 0;
     kernels->spatialOrder = spatialOrder;
+    kernels->bgOrder = 0;
 
     return kernels;
@@ -71,5 +69,5 @@
     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
 
-    int num = PS_SQR(2 * size + 1) * (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions
+    int num = PS_SQR(2 * size + 1); // Number of basis functions
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS,
@@ -81,26 +79,15 @@
     // Generate a set of kernels for each (u,v)
     for (int v = - size, index = 0; v <= size; v++) {
-        for (int u = - size; u <= size; u++) {
-            // Iterate over spatial order.  This loop creates the terms for
-            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
-            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
-                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
-                    kernels->u->data.S32[index] = u;
-                    kernels->v->data.S32[index] = v;
-                    kernels->xOrder->data.S32[index] = xOrder;
-                    kernels->yOrder->data.S32[index] = yOrder;
-
-                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
-                            u, v, xOrder, yOrder);
-                }
-            }
-        }
-    }
-
-    kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2;
+        for (int u = - size; u <= size; u++, index++) {
+            kernels->u->data.S32[index] = u;
+            kernels->v->data.S32[index] = v;
+
+            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
+        }
+    }
+
+    kernels->subIndex = num / 2;
     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);
+           kernels->v->data.S32[kernels->subIndex] == 0);
 
     return kernels;
@@ -127,5 +114,4 @@
         num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
     }
-    num *= (spatialOrder + 1) * (spatialOrder + 2) / 2;
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS,
@@ -138,10 +124,12 @@
     kernels->widths = psVectorAlloc(num, PS_TYPE_F32);
     kernels->preCalc = psArrayAlloc(num);
+    psKernel *subtract = NULL;          // Kernel to subtract to maintain flux scaling
 
     // Set the kernel parameters
     for (int i = 0, index = 0; i < numGaussians; i++) {
+        float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian
         // Iterate over (u,v) order
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
-            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
 
                 // Set the pre-calculated kernel
@@ -150,37 +138,36 @@
                 for (int v = -size; v <= size; v++) {
                     for (int u = -size; u <= size; u++) {
-                        sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
+                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
                     }
                 }
-                // Normalise sum of kernel component to unity
-                psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32));
-
-                // Iterate over spatial order.  This loop creates the terms for
-                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
-                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
-                    for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
-                        kernels->widths->data.F32[index] = sigmas->data.F32[i];
-                        kernels->u->data.S32[index] = uOrder;
-                        kernels->v->data.S32[index] = vOrder;
-                        kernels->xOrder->data.S32[index] = xOrder;
-                        kernels->yOrder->data.S32[index] = yOrder;
-                        kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc);
-
-                        psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index,
-                                sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder);
+                if (index == 0) {
+                    subtract = preCalc;
+                    for (int v = -size; v <= size; v++) {
+                        for (int u = -size; u <= size; u++) {
+                            preCalc->kernel[v][u] /= sum;
+                        }
+                    }
+                } else if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+                    // Normalise sum of kernel component to unity for even functions
+                    for (int v = -size; v <= size; v++) {
+                        for (int u = -size; u <= size; u++) {
+                            preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u];
+                        }
                     }
                 }
 
-                psFree(preCalc);        // Drop reference
+                kernels->widths->data.F32[index] = sigmas->data.F32[i];
+                kernels->u->data.S32[index] = uOrder;
+                kernels->v->data.S32[index] = vOrder;
+                kernels->preCalc->data[index] = preCalc;
+
+                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
+                        sigmas->data.F32[i], uOrder, vOrder);
             }
         }
     }
 
-    kernels->subIndex = 0;
-    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);
+    kernels->subIndex = -1;
 
     if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
@@ -220,6 +207,5 @@
     psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
 
-    int num = PS_SQR(2 * numTotal + 1) *
-        (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions
+    int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
 
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
@@ -264,33 +250,23 @@
         int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel
 
-        for (int j = - numTotal; j <= numTotal; j++) {
+        for (int j = - numTotal; j <= numTotal; j++, index++) {
             int v = locations->data.S32[numTotal + j]; // Location of pixel
             int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel
 
-            // Iterate over spatial order.  This loop creates the terms for
-            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
-            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
-                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
-                    kernels->u->data.S32[index] = u;
-                    kernels->v->data.S32[index] = v;
-                    kernels->uStop->data.S32[index] = uStop;
-                    kernels->vStop->data.S32[index] = vStop;
-                    kernels->xOrder->data.S32[index] = xOrder;
-                    kernels->yOrder->data.S32[index] = yOrder;
-
-                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index,
-                            u, uStop, v, vStop, xOrder, yOrder);
-                }
-            }
-        }
-    }
-
-    kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2;
+            kernels->u->data.S32[index] = u;
+            kernels->v->data.S32[index] = v;
+            kernels->uStop->data.S32[index] = uStop;
+            kernels->vStop->data.S32[index] = vStop;
+
+            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
+                    u, uStop, v, vStop);
+        }
+    }
+
+    kernels->subIndex = num / 2;
     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
            kernels->v->data.S32[kernels->subIndex] == 0 &&
            kernels->uStop->data.S32[kernels->subIndex] == 0 &&
-           kernels->vStop->data.S32[kernels->subIndex] == 0 &&
-           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
-           kernels->yOrder->data.S32[kernels->subIndex] == 0);
+           kernels->vStop->data.S32[kernels->subIndex] == 0);
 
     psFree(locations);
@@ -327,6 +303,5 @@
     psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
 
-    int num = PS_SQR(2 * numTotal + 1) *
-        (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions
+    int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
 
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
@@ -369,33 +344,23 @@
         int u = start->data.S32[numTotal + i]; // Location of pixel
         int uStop = stop->data.S32[numTotal + i]; // Width of pixel
-        for (int j = - numTotal; j <= numTotal; j++) {
+        for (int j = - numTotal; j <= numTotal; j++, index++) {
             int v = start->data.S32[numTotal + j]; // Location of pixel
             int vStop = stop->data.S32[numTotal + j]; // Width of pixel
 
-            // Iterate over spatial order.  This loop creates the terms for
-            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
-            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
-                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
-                    kernels->u->data.S32[index] = u;
-                    kernels->v->data.S32[index] = v;
-                    kernels->uStop->data.S32[index] = uStop;
-                    kernels->vStop->data.S32[index] = vStop;
-                    kernels->xOrder->data.S32[index] = xOrder;
-                    kernels->yOrder->data.S32[index] = yOrder;
-
-                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index,
-                            u, uStop, v, vStop, xOrder, yOrder);
-                }
-            }
-        }
-    }
-
-    kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2;
+            kernels->u->data.S32[index] = u;
+            kernels->v->data.S32[index] = v;
+            kernels->uStop->data.S32[index] = uStop;
+            kernels->vStop->data.S32[index] = vStop;
+
+            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
+                    u, uStop, v, vStop);
+        }
+    }
+
+    kernels->subIndex = num / 2;
     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
            kernels->v->data.S32[kernels->subIndex] == 0 &&
            kernels->uStop->data.S32[kernels->subIndex] == 0 &&
-           kernels->vStop->data.S32[kernels->subIndex] == 0 &&
-           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
-           kernels->yOrder->data.S32[kernels->subIndex] == 0);
+           kernels->vStop->data.S32[kernels->subIndex] == 0);
 
     psFree(start);
@@ -429,7 +394,6 @@
 
     int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements
-    int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel
-
-    int num = (numGaussianVars + numInner) * numSpatial; // Total number of basis functions
+
+    int num = numGaussianVars + numInner; // Total number of basis functions
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK,
@@ -440,15 +404,15 @@
     psFree(params);
 
-    kernels->widths = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32);
-    kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial);
-    kernels->inner = numGaussianVars * numSpatial;
+    kernels->widths = psVectorAlloc(numGaussianVars, PS_TYPE_F32);
+    kernels->preCalc = psArrayAlloc(numGaussianVars);
+    kernels->inner = numGaussianVars;
+    psKernel *subtract = NULL;          // Kernel to subtract to maintain flux scaling
 
     // Set the Gaussian kernel parameters
     for (int i = 0, index = 0; i < numGaussians; i++) {
+        float norm = 1.0 / (M_2_PI * sqrtf(sigmas->data.F32[i])); // Normalisation for Gaussian
         // Iterate over (u,v) order
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
-            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) {
-
-
+            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
                 // Set the pre-calculated kernel
                 psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
@@ -456,57 +420,48 @@
                 for (int v = -size; v <= size; v++) {
                     for (int u = -size; u <= size; u++) {
-                        sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
+                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
                     }
                 }
-                // Normalise sum of kernel component to unity
-                psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32));
-
-                // Iterate over spatial order.  This loop creates the terms for
-                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
-                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
-                    for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
-                        kernels->widths->data.F32[index] = sigmas->data.F32[i];
-                        kernels->u->data.S32[index] = uOrder;
-                        kernels->v->data.S32[index] = vOrder;
-                        kernels->xOrder->data.S32[index] = xOrder;
-                        kernels->yOrder->data.S32[index] = yOrder;
-                        kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc);
-
-                        psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index,
-                                sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder);
+                if (index == 0) {
+                    subtract = preCalc;
+                    for (int v = -size; v <= size; v++) {
+                        for (int u = -size; u <= size; u++) {
+                            preCalc->kernel[v][u] /= sum;
+                        }
+                    }
+                } else if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+                    // Normalise sum of kernel component to unity for even functions
+                    for (int v = -size; v <= size; v++) {
+                        for (int u = -size; u <= size; u++) {
+                            preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum - subtract->kernel[v][u];
+                        }
                     }
                 }
 
-                psFree(preCalc);        // Drop reference
+                kernels->widths->data.F32[index] = sigmas->data.F32[i];
+                kernels->u->data.S32[index] = uOrder;
+                kernels->v->data.S32[index] = vOrder;
+                kernels->preCalc->data[index] = preCalc;
+
+                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
+                        sigmas->data.F32[i], uOrder, vOrder);
             }
         }
     }
-
 
     // Generate a grid set of kernels for each (u,v)
     for (int v = - inner, index = kernels->inner; v <= inner; v++) {
-        for (int u = - inner; u <= inner; u++) {
-            // Iterate over spatial order.  This loop creates the terms for
-            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
-            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
-                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
-                    kernels->u->data.S32[index] = u;
-                    kernels->v->data.S32[index] = v;
-                    kernels->xOrder->data.S32[index] = xOrder;
-                    kernels->yOrder->data.S32[index] = yOrder;
-
-                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
-                            u, v, xOrder, yOrder);
-                }
-            }
-        }
-    }
-
-    kernels->subIndex = kernels->inner + (numInner - 1) * numSpatial / 2;
+        for (int u = - inner; u <= inner; u++, index++) {
+            kernels->u->data.S32[index] = u;
+            kernels->v->data.S32[index] = v;
+
+            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
+        }
+    }
+
+    kernels->subIndex = inner + num / 2;
     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);
+           kernels->v->data.S32[kernels->subIndex] == 0);
 
     return kernels;
@@ -543,7 +498,6 @@
     int numRings = numOuter + numInner; // Number of rings (not including the central pixel)
     int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring
-    int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel
-
-    int num = (numRings * numPoly + 1) * numSpatial; // Total number of basis functions
+
+    int num = numRings * numPoly + 1; // Total number of basis functions
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS,
@@ -589,5 +543,5 @@
         // Iterate over (u,v) order
         for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) {
-            for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++) {
+            for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++, index++) {
 
                 psArray *data = psArrayAlloc(3); // Container for data
@@ -604,4 +558,5 @@
                 } else {
                     int j = 0;          // Index for data
+                    double norm = 0.0;  // Normalisation
                     for (int v = -size; v <= size; v++) {
                         int v2 = PS_SQR(v);   // Square of v
@@ -616,5 +571,5 @@
                                 uCoords->data.S32[j] = u;
                                 vCoords->data.S32[j] = v;
-                                poly->data.F32[j] = uPoly * vPoly;
+                                norm += poly->data.F32[j] = uPoly * vPoly;
 
                                 psVectorExtend(uCoords, RINGS_BUFFER, 1);
@@ -629,23 +584,17 @@
                         }
                     }
+                    // Normalise kernel component to unit sum
+                    psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
+
                 }
 
                 psTrace("psModules.imcombine", 8, "%ld pixels in ring\n", uCoords->n);
 
-                // Iterate over spatial order.  This loop creates the terms for
-                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
-                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
-                    for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
-                        kernels->preCalc->data[index] = psMemIncrRefCounter(data);
-                        kernels->u->data.S32[index] = uOrder;
-                        kernels->v->data.S32[index] = vOrder;
-                        kernels->xOrder->data.S32[index] = xOrder;
-                        kernels->yOrder->data.S32[index] = yOrder;
-
-                        psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d\n", index,
-                                i, uOrder, vOrder, xOrder, yOrder);
-                    }
-                }
-                psFree(data);
+                kernels->preCalc->data[index] = data;
+                kernels->u->data.S32[index] = uOrder;
+                kernels->v->data.S32[index] = vOrder;
+
+                psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,
+                        i, uOrder, vOrder);
             }
         }
@@ -653,6 +602,4 @@
 
     kernels->subIndex = 0;
-    assert(kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
-           kernels->yOrder->data.S32[kernels->subIndex] == 0);
 
     return kernels;