Index: branches/pap/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- branches/pap/psModules/src/imcombine/pmSubtraction.c	(revision 25834)
+++ branches/pap/psModules/src/imcombine/pmSubtraction.c	(revision 25841)
@@ -70,4 +70,5 @@
                               const pmSubtractionKernels *kernels, // Kernel basis functions
                               const psImage *polyValues, // Spatial polynomial values
+                              bool normalise,            // Add normalisation?
                               bool wantDual // Want the dual (second) kernel?
                               )
@@ -178,6 +179,8 @@
     }
 
-    // Put in the normalisation component
-    kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels));
+    if (normalise) {
+        // Put in the normalisation component
+        kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels));
+    }
 
     return kernel;
@@ -349,5 +352,5 @@
     )
 {
-    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
+    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, true, wantDual);
     if (variance || subMask) {
         *kernelVariance = varianceKernel(*kernelVariance, *kernelImage);
@@ -466,5 +469,5 @@
     psKernel *kernel;                   // Kernel to use
     if (!preKernel) {
-        kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
+        kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual);
     } else {
         kernel = psMemIncrRefCounter(preKernel);
@@ -916,5 +919,5 @@
 
     psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial
-    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel
+    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual); // The appropriate kernel
     psFree(polyValues);
 
@@ -947,5 +950,5 @@
     psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y);
 
-    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel
+    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual); // The appropriate kernel
     psFree(polyValues);
 
@@ -977,12 +980,34 @@
     int num = wantDual ? solution->n - 1 : solution->n; // Number of kernel basis functions
 
-    psArray *images = psArrayAlloc(num); // Images of each kernel to return
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial
+    psArray *images = psArrayAlloc(num + 1); // Images of each kernel to return
+
+    // The whole kernel
+    {
+        psKernel *kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual); // The appropriate kernel
+        images->data[0] = psMemIncrRefCounter(kernel->image);
+        psFree(kernel);
+    }
+
+    // The parts
     psVectorInit(solution, 0.0);
-
     for (int i = 0; i < num; i++) {
         solution->data.F64[i] = backup->data.F64[i];
-        images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual);
+        psKernel *kernel = solvedKernel(NULL, kernels, polyValues, false, wantDual); // The appropriate kernel
+#if 0
+        int size = kernels->size;
+        double sum = 0.0;
+        for (int v = -size; v <= size; v++) {
+            for (int u = -size; u <= size; u++) {
+                sum += kernel->kernel[v][u];
+            }
+        }
+        fprintf(stderr, "Kernel %d: %lf\n", i, sum);
+#endif
+        images->data[i + 1] = psMemIncrRefCounter(kernel->image);
+        psFree(kernel);
         solution->data.F64[i] = 0.0;
     }
+    psFree(polyValues);
     psVectorCopy(solution, backup, PS_TYPE_F64);
     psFree(backup);