Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 15756)
+++ trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 15761)
@@ -358,6 +358,9 @@
             return false;
         }
-
-        matrix->data.F64[i][bgIndex] = sumC;
+        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
+            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
+                matrix->data.F64[index][bgIndex] = sumC * polyValues->data.F64[yOrder][xOrder];
+            }
+        }
     }
 
@@ -582,6 +585,6 @@
 
     // Check inputs
-    long numParams = -1;                // Number of parameters
-    long numParams2 = 0;                // Number of parameters for part solution (DUAL mode)
+    int numParams = -1;                // Number of parameters
+    int numParams2 = 0;                // Number of parameters for part solution (DUAL mode)
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -601,11 +604,10 @@
             PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false);
             PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false);
-            // XXX Not sure I've got this the right way around
-            PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams2, numParams, false);
+            PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, false);
             PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false);
             PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false);
             PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false);
             PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false);
-            PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, false);
+            PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
             PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false);
         } else {
@@ -614,5 +616,5 @@
                 numParams = stamp->vector->n;
             }
-            PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, false);
+            PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
             PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false);
             PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false);
@@ -681,5 +683,5 @@
         }
 
-        psVector *permutation1, *permutation2; // Permutation vectors for LUD
+        psVector *permutation1 = NULL, *permutation2 = NULL; // Permutation vectors for LUD
         psImage *lu1 = psMatrixLUD(NULL, &permutation1, sumMatrix1); // LUD of matrix 1
         psImage *lu2 = psMatrixLUD(NULL, &permutation2, sumMatrix2); // LUD of matrix 2