Index: branches/pap/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- branches/pap/psModules/src/imcombine/pmSubtractionEquation.c	(revision 25897)
+++ branches/pap/psModules/src/imcombine/pmSubtractionEquation.c	(revision 25899)
@@ -968,6 +968,4 @@
 
         psVector *solution = NULL;                       // Solution to equation!
-        psVector *mask = psVectorAlloc(numParams + numParams2, PS_TYPE_U8); // Mask of parameters
-        psVectorInit(mask, 0);
         {
             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
@@ -976,38 +974,42 @@
                 psFree(sumMatrix);
                 psFree(sumVector);
-                psFree(mask);
                 return NULL;
             }
 
+            int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
+            int numKernels = kernels->num; // Number of kernel basis functions
+
             // Remove a kernel basis for image 1 from the equation
-#define MASK_BASIS_1(INDEX) \
-            { \
-                for (int k = 0; k < numParams2; k++) { \
-                    sumMatrix1->data.F64[k][INDEX] = 0.0; \
-                    sumMatrix1->data.F64[INDEX][k] = 0.0; \
-                    sumMatrixX->data.F64[k][INDEX] = 0.0; \
-                } \
-                sumMatrix1->data.F64[bgIndex][INDEX] = 0.0; \
-                sumMatrix1->data.F64[INDEX][bgIndex] = 0.0; \
-                sumMatrix1->data.F64[normIndex][INDEX] = 0.0; \
-                sumMatrix1->data.F64[INDEX][normIndex] = 0.0; \
-                sumMatrix1->data.F64[INDEX][INDEX] = 1.0; \
-                sumVector1->data.F64[INDEX] = 0.0; \
-                mask->data.U8[INDEX] = 0xFF; \
+#define MASK_BASIS_1(INDEX)                                             \
+            {                                                           \
+                for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
+                    for (int k = 0; k < numParams2; k++) {              \
+                        sumMatrix1->data.F64[k][index] = 0.0;           \
+                        sumMatrix1->data.F64[index][k] = 0.0;           \
+                        sumMatrixX->data.F64[k][index] = 0.0;           \
+                    }                                                   \
+                    sumMatrix1->data.F64[bgIndex][index] = 0.0;         \
+                    sumMatrix1->data.F64[index][bgIndex] = 0.0;         \
+                    sumMatrix1->data.F64[normIndex][index] = 0.0;       \
+                    sumMatrix1->data.F64[index][normIndex] = 0.0;       \
+                    sumMatrix1->data.F64[index][index] = 1.0;           \
+                    sumVector1->data.F64[index] = 0.0;                  \
+                }                                                       \
             }
 
             // Remove a kernel basis for image 2 from the equation
-#define MASK_BASIS_2(INDEX) \
-            { \
-                for (int k = 0; k < numParams2; k++) { \
-                    sumMatrix2->data.F64[k][j] = 0.0; \
-                    sumMatrix2->data.F64[j][k] = 0.0; \
-                    sumMatrixX->data.F64[j][k] = 0.0; \
-                } \
-                sumMatrix2->data.F64[INDEX][INDEX] = 1.0; \
-                sumMatrixX->data.F64[j][normIndex] = 0.0; \
-                sumMatrixX->data.F64[j][bgIndex] = 0.0; \
-                sumVector2->data.F64[j] = 0.0; \
-                mask->data.U8[numParams + j] = 0xFF; \
+#define MASK_BASIS_2(INDEX)                                             \
+            {                                                           \
+                for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
+                    for (int k = 0; k < numParams2; k++) {              \
+                        sumMatrix2->data.F64[k][index] = 0.0;           \
+                        sumMatrix2->data.F64[index][k] = 0.0;           \
+                        sumMatrixX->data.F64[index][k] = 0.0;           \
+                    }                                                   \
+                    sumMatrix2->data.F64[index][index] = 1.0;           \
+                    sumMatrixX->data.F64[index][normIndex] = 0.0;       \
+                    sumMatrixX->data.F64[index][bgIndex] = 0.0;         \
+                    sumVector2->data.F64[index] = 0.0;                  \
+                }                                                       \
             }
 
@@ -1016,23 +1018,30 @@
             double norm = solution->data.F64[normIndex];        // Normalisation
             double thresh = norm * TOL;                         // Threshold for low parameters
-            for (int j = 0; j < numParams2; j++) {
-                double param1 = solution->data.F64[j],
-                    param2 = solution->data.F64[numParams + j]; // Parameters of interest
+            for (int i = 0; i < numKernels; i++) {
+                // Getting 0th order parameter value.  In the presence of spatial variation, the actual value
+                // of the parameter will vary over the image.  We are in effect getting the value in the
+                // centre of the image.  If we use different polynomial functions (e.g., Chebyshev), we may
+                // have to change this to properly determine the value of the parameter at the centre.
+                double param1 = solution->data.F64[i],
+                    param2 = solution->data.F64[numParams + i]; // Parameters of interest
+                bool mask1 = false, mask2 = false;              // Masked the parameter?
                 if (fabs(param1) < thresh) {
-                    psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", j);
-                    MASK_BASIS_1(j);
+                    psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i);
+                    MASK_BASIS_1(i);
+                    mask1 = true;
                 }
                 if (fabs(param2) < thresh) {
-                    psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", j);
-                    MASK_BASIS_2(j);
-                }
-
-                if (!mask->data.U8[j] && !mask->data.U8[numParams + j]) {
+                    psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i);
+                    MASK_BASIS_2(i);
+                    mask2 = true;
+                }
+
+                if (!mask1 && !mask2) {
                     if (fabs(param1) < fabs(param2)) {
-                        psTrace("psModules.imcombine", 7, "Parameter %d: 1 < 2\n", j);
-                        MASK_BASIS_1(j);
+                        psTrace("psModules.imcombine", 7, "Parameter %d: 1 < 2\n", i);
+                        MASK_BASIS_1(i);
                     } else {
-                        psTrace("psModules.imcombine", 7, "Parameter %d: 2 < 1\n", j);
-                        MASK_BASIS_2(j);
+                        psTrace("psModules.imcombine", 7, "Parameter %d: 2 < 1\n", i);
+                        MASK_BASIS_2(i);
                     }
                 }
@@ -1066,16 +1075,6 @@
             psFree(sumMatrix);
             psFree(sumVector);
-            psFree(mask);
             return NULL;
         }
-
-#if 0
-        for (int i = 0; i < num; i++) {
-            if (mask->data.U8[i]) {
-                solution->data.F64[i] = 0.0;
-            }
-        }
-#endif
-        psFree(mask);
 
         psFree(sumMatrix1);