Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14325)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14329)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-19 21:42:12 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-20 01:53:41 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -41,14 +41,23 @@
 
     psImage *polyValues = psImageAlloc(spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); // Matrix to return
-    psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, spatialOrder, spatialOrder); // Polynomial
-
-    for (int j = 0; j <= spatialOrder; j++) {
-        for (int i = 0; i <= spatialOrder - j; i++) {
-            poly->coeff[i][j] = 1.0;
-            polyValues->data.F64[j][i] = psPolynomial2DEval(poly, x, y);
-            poly->coeff[i][j] = 0.0;
-        }
-    }
-    psFree(poly);
+    polyValues->data.F64[0][0] = 1.0;
+
+    double value = 1.0;
+    for (int i = 1; i <= spatialOrder; i++) {
+        value *= x;
+        polyValues->data.F64[0][i] = value;
+    }
+
+    value = 1.0;
+    for (int j = 1; j <= spatialOrder; j++) {
+        value *= y;
+        polyValues->data.F64[j][0] = value;
+    }
+
+    for (int j = 1; j <= spatialOrder; j++) {
+        for (int i = 1; i <= spatialOrder - j; i++) {
+            polyValues->data.F64[j][i] = polyValues->data.F64[j][0] * polyValues->data.F64[0][i];
+        }
+    }
 
     return polyValues;
@@ -875,4 +884,37 @@
 }
 
+psArray *pmSubtractionKernelSolutions(const psVector *solution, const pmSubtractionKernels *kernels,
+                                      float x, float y
+                                      )
+{
+    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
+    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
+    PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL);
+    if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
+        PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, NULL);
+        PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, NULL);
+    }
+
+    psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return
+    psVector *fakeSolution = psVectorAlloc(solution->n, PS_TYPE_F64); // Fake solution vector
+    psVectorInit(fakeSolution, 0.0);
+
+    for (int i = 0; i < solution->n - 1; i++) {
+        fakeSolution->data.F64[i] = solution->data.F64[i];
+        images->data[i] = pmSubtractionKernelImage(fakeSolution, kernels, x, y);
+        fakeSolution->data.F64[i] = 0.0;
+    }
+
+    psFree(fakeSolution);
+
+    return images;
+}
+
+
 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
                            const psImage *inImage, const psImage *inWeight, const psImage *subMask,
@@ -953,6 +995,6 @@
     psKernel *kernelWeight = NULL;      // Kernel for the weight map
 
-    for (int j = size; j < inImage->numRows - fullSize - size; j += fullSize) {
-        for (int i = size; i < inImage->numCols - fullSize - size; i += fullSize) {
+    for (int j = size; j < inImage->numRows - size; j += fullSize) {
+        for (int i = size; i < inImage->numCols - size; i += fullSize) {
 
             // Only generate polynomial values every kernel footprint, since we have already assumed
@@ -960,6 +1002,6 @@
             psFree(polyValues);
             polyValues = spatialPolyValues(kernels->spatialOrder,
-                                           2.0 * (float)(i + size - numCols/2.0) / (float)numCols,
-                                           2.0 * (float)(j + size - numRows/2.0) / (float)numRows);
+                                           2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols,
+                                           2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows);
 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
             kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
@@ -975,6 +1017,6 @@
 #endif
 
-            for (int y = j; y < j + fullSize; y++) {
-                for (int x = i; x < i + fullSize; x++) {
+            for (int y = j; y < PS_MIN(j + fullSize, numRows - size); y++) {
+                for (int x = i; x < PS_MIN(i + fullSize, numCols - size); x++) {
                     // Check and propagate the kernel footprint, if required
                     if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &