Changeset 14455 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Aug 9, 2007, 10:25:52 AM (19 years ago)

Author:

Paul Price

Message:

Extensive changes to the way spatial variation of the kernel is implemented. I was doing it the stupid way, rather than the smart way outlined in Alard 2000 (A&ASS, 144, 363). Fixed up the normalisation of the kernels (essential for proper flux conservation with spatial variation). Put in support for dividing an image up into regions.

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (33 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-07 19:02:55 $
+ *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-09 20:25:52 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -30 +30 @@
 // Private (file-static) functions
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// Return the number of polynomial terms there are for a given order
+static inline long polyTerms(int order  // Order of polynomial
+                             )
+{
+    return (order + 1) * (order + 2) / 2;
+}
 
 // Given a stamp coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
@@ -160 +167 @@
           case PM_SUBTRACTION_KERNEL_ISIS: { \
               psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \
-              psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \
               for (int v = -size; v <= size; v++) { \
                   for (int u = -size; u <= size; u++) { \
                       (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \
-                      /* Subtract (0,0) kernel from other kernels to preserve photometric scaling */ \
-                      if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
-                          (TARGET)->kernel[v][u] += subValue * FUNC(subKernel->kernel[v][u]); \
-                      } \
+                      /* Photometric scaling is already built in to the precalculated kernel */ \
                   } \
               } \
@@ -207 +210 @@
     )
 {
-    int numKernels = solution->n - 1;   // Number of kernel basis functions
+    int numKernels = kernels->num;      // Number of kernel basis functions
     assert(kernels->u->n == numKernels);
     assert(kernels->v->n == numKernels);
-    assert(kernels->xOrder->n == numKernels);
-    assert(kernels->yOrder->n == numKernels);
+    int spatialOrder = kernels->spatialOrder; // Maximum spatial polynomial order
+    assert(solution->n == numKernels * polyTerms(spatialOrder) + polyTerms(kernels->bgOrder));
 
     // Ensure the subIndex for POIS kernels is what is expected
     assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS &&
             kernels->type != PM_SUBTRACTION_KERNEL_SPAM &&
-            kernels->type != PM_SUBTRACTION_KERNEL_FRIES) ||
-           (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->type != PM_SUBTRACTION_KERNEL_FRIES &&
+            kernels->type != PM_SUBTRACTION_KERNEL_GUNK &&
+            kernels->type != PM_SUBTRACTION_KERNEL_RINGS) ||
+           (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0));
 
     int size = kernels->size;           // Kernel half-size
@@ -228 +231 @@
 
     for (int i = 0; i < numKernels; i++) {
-        int xOrder = kernels->xOrder->data.S32[i]; // Polynomial order in x
-        int yOrder = kernels->yOrder->data.S32[i]; // Polynomial order in y
-        float value = weightFunc(polyValues->data.F64[yOrder][xOrder] *
-                                 solution->data.F64[i]); // Value to sum
-        float subValue = weightFunc(-solution->data.F64[i]); // Value to subtract (actually add, because of -)
+        double value = 0.0;              // Value of kernel coefficient
+        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
+            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
+                value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
+            }
+        }
+
         switch (kernels->type) {
           case PM_SUBTRACTION_KERNEL_POIS: {
               int u = kernels->u->data.S32[i]; // Offset in x
               int v = kernels->v->data.S32[i]; // Offset in y
-              kernel->kernel[v][u] += value;
+              kernel->kernel[v][u] += weightFunc(value);
               if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
                   // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                  kernel->kernel[0][0] += subValue;
+                  kernel->kernel[0][0] += weightFunc(- value);
               }
               break;
@@ -253 +258 @@
 
               // Normalising sum of kernel component to unity
-              value /= weightFunc((uStop - uStart) * (vStop - vStart));
+              float norm = 1.0 / weightFunc((uStop - uStart) * (vStop - vStart));
 
               for (int v = vStart; v <= vStop; v++) {
                   for (int u = uStart; u <= uStop; u++) {
-                      kernel->kernel[v][u] += value;
+                      kernel->kernel[v][u] += norm * weightFunc(value);
+                      if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
+                          // The (0,0) element is subtracted from most kernels to preserve photometric scaling
+                          kernel->kernel[0][0] += weightFunc(- value);
+                      }
                   }
-              }
-              if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
-                  /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */
-                  kernel->kernel[0][0] += subValue;
               }
               break;
@@ -270 +275 @@
                   // Using pre-calculated function
                   psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values
+                  // Iterating over the kernel
                   for (int v = -size; v <= size; v++) {
                       for (int u = -size; u <= size; u++) {
-                          kernel->kernel[v][u] += value * weightFunc(preCalc->kernel[v][u]);
+                          kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value);
+                          // Photometric scaling is built into the preCalc kernel --- no subtraction!
                       }
                   }
               } else {
                   // Using delta function
+                  bool subtract = (kernels->spatialOrder > 0 && i != kernels->subIndex); // Subtract (0,0)?
                   int u = kernels->u->data.S32[i]; // Offset in x
                   int v = kernels->v->data.S32[i]; // Offset in y
-                  kernel->kernel[v][u] += value;
-              }
-              // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
-              if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
-                  kernel->kernel[0][0] += subValue;
-              }
-              break;
-          }
-          case PM_SUBTRACTION_KERNEL_ISIS: {
-              psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values
-              psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract
-              for (int v = -size; v <= size; v++) {
-                  for (int u = -size; u <= size; u++) {
-                      kernel->kernel[v][u] += value * weightFunc(preCalc->kernel[v][u]);
-                      // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
-                      if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
-                          kernel->kernel[v][u] += subValue * weightFunc(subKernel->kernel[v][u]);
-                      }
+                  kernel->kernel[v][u] += weightFunc(value);
+                  if (subtract) {
+                      // The (0,0) element is subtracted from most kernels to preserve photometric scaling
+                      kernel->kernel[0][0] += weightFunc(- value);
                   }
               }
               break;
           }
+          case PM_SUBTRACTION_KERNEL_ISIS: {
+              psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values
+              // Iterating over the kernel
+              for (int v = -size; v <= size; v++) {
+                  for (int u = -size; u <= size; u++) {
+                      kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value);
+                      // Photometric scaling is built into the preCalc kernel --- no subtraction!
+                  }
+              }
+              break;
+          }
           case PM_SUBTRACTION_KERNEL_RINGS: {
               if (i == kernels->subIndex) {
-                  kernel->kernel[0][0] += value;
+                  kernel->kernel[0][0] += weightFunc(value);
                   break;
               }
@@ -311 +316 @@
               psVector *poly = preCalc->data[2]; // Polynomial values
               int num = uCoords->n;     // Number of pixels
-              value /= weightFunc(num); // Normalising sum of kernel component to unity
+
               for (int j = 0; j < num; j++) {
                   int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
-                  kernel->kernel[v][u] += value * weightFunc(poly->data.F32[j]);
+                  kernel->kernel[v][u] += weightFunc(value * poly->data.F32[j]);
               }
-
-              // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
-              if (kernels->spatialOrder > 0) {
-                  kernel->kernel[0][0] += subValue;
-              }
+              kernel->kernel[0][0] += weightFunc(- value * num);
               break;
           }
@@ -336 +337 @@
                                    int index, // Kernel basis function index
                                    int x, int y, // Pixel around which to convolve
-                                   const psImage *image, // Image to convolve
-                                   const psImage *polyValues // Spatial polynomial values
+                                   const psImage *image // Image to convolve
                                    )
 {
-    int xOrder = kernels->xOrder->data.S32[index]; // Polynomial order in x
-    int yOrder = kernels->yOrder->data.S32[index]; // Polynomial order in y
-    double polyValue = polyValues->data.F64[yOrder][xOrder]; // Value of spatial polynomial
-
     switch (kernels->type) {
       case PM_SUBTRACTION_KERNEL_POIS: {
@@ -349 +345 @@
           int u = kernels->u->data.S32[index]; // Offset in x
           int v = kernels->v->data.S32[index]; // Offset in y
-          double value = polyValue * image->data.F32[y + v][x + u]; // Value of convolution
+          float value = image->data.F32[y + v][x + u]; // Value of convolution
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
               // The (0,0) element is subtracted from most kernels to preserve photometric scaling
@@ -366 +362 @@
           for (int v = vStart; v <= vStop; v++) {
               for (int u = uStart; u <= uStop; u++) {
-                  sum += polyValue * image->data.F32[y + v][x + u];
+                  sum += image->data.F32[y + v][x + u];
               }
           }
@@ -377 +373 @@
       }
       case PM_SUBTRACTION_KERNEL_GUNK: {
-          double value;                 // The value to return
           if (index < kernels->inner) {
               // Using pre-calculated function
@@ -388 +383 @@
                   }
               }
-              value = polyValue * sum;
-          } else {
-              // Using delta function
-              int u = kernels->u->data.S32[index]; // Offset in x
-              int v = kernels->v->data.S32[index]; // Offset in y
-              value = polyValue * image->data.F32[y + v][x + u]; // Value of convolution
-          }
+              return sum;
+          }
+          // Using delta function
+          int u = kernels->u->data.S32[index]; // Offset in x
+          int v = kernels->v->data.S32[index]; // Offset in y
+          float value = image->data.F32[y + v][x + u]; // Value of convolution
           // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
@@ -403 +397 @@
       case PM_SUBTRACTION_KERNEL_ISIS: {
           psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel
-          psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract
           int size = kernels->size;     // Kernel half-size
           double sum = 0.0;             // Accumulated sum from convolution
-          double sub = 0.0;             // Accumulated sum to subtract
           for (int v = -size; v <= size; v++) {
               for (int u = -size; u <= size; u++) {
                   sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u];
-                  // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
-                  if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-                      sub += subKernel->kernel[v][u] * image->data.F32[y + v][x + u];
-                  }
+                  // Photometric scaling is already built in to the precalculated kernel
               }
           }
-          return polyValue * sum - sub;
+          return sum;
       }
       case PM_SUBTRACTION_KERNEL_RINGS: {
           if (index == kernels->subIndex) {
-              return image->data.F32[0][0];
+              return image->data.F32[y][x];
           }
           psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
@@ -432 +421 @@
               sum += image->data.F32[y + v][x + u] * poly->data.F32[j];
           }
-          sum /= (double)num;           // Normalising sum of kernel component to unity
           // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
-          return polyValue * sum - image->data.F32[0][0];
+          return sum - num * image->data.F32[y][x];
       }
       default:
@@ -592 +580 @@
     PS_ASSERT_PTR_NON_NULL(kernels, false);
     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false);
     if (weight) {
         PS_ASSERT_IMAGE_NON_NULL(weight, false);
@@ -599 +585 @@
         PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
     }
-    if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->preCalc, false);
-    }
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
 
@@ -609 +592 @@
 
     int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation
-    int numKernels = kernels->u->n;     // Number of kernel basis functions
-    int numParams = numKernels + 1;     // Total number of parameters to solve for: coefficient of each kernel
-                                        // basis function, and a constant background offset.
-    int bgIndex = numKernels;           // Index in matrix for the background
-    psVector *convolutions = psVectorAlloc(numKernels, PS_TYPE_F64); // Convolutions for each kernel
+    int numKernels = kernels->num;      // Number of kernel basis functions
+    int numSpatial = polyTerms(spatialOrder); // Number of spatial variations
+    int numBackground = polyTerms(kernels->bgOrder); // Number of background terms
+
+    // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the
+    // number of coefficients for the spatial polynomial, and a constant background offset.
+    int numParams = numKernels * numSpatial + numBackground;
+    int bgIndex = numParams - numBackground; // Index in matrix for the background
+    psVector *convolutions = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); // Convolutions
 
     // We iterate over each stamp, allocate the matrix and vectors if
@@ -637 +624 @@
             psVectorInit(stampVector, 0.0);
 
-            // Pre-evaluated spatial polynomial values
-            psImage *polyValues = spatialPolyValues(spatialOrder,
-                                                    2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols,
-                                                    2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows);
+            float xNorm = 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols; // Normalised x coord
+            float yNorm = 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows; // Normalised y coord
+            psImage *polyValues = spatialPolyValues(spatialOrder, xNorm, yNorm); // Spatial polynomial terms
 
             for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
@@ -648 +634 @@
 
                     // Generate the convolutions
-                    for (int i = 0; i < numKernels; i++) {
-                        convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues);
+                    for (int j = 0; j < numKernels; j++) {
+                        double value = convolvePixel(kernels, j, x, y, reference); // Value from convolution
+                        // Generate the pseudo-convolutions from the spatial polynomial terms
+                        for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) {
+                            for (int xOrder = 0; xOrder <= spatialOrder - yOrder;
+                                 xOrder++, index += numKernels) {
+                                convolutions->data.F64[index] = value * polyValues->data.F64[yOrder][xOrder];
+                            }
+                        }
                     }
 
                     // Generate the least-squares matrix and vector
                     // Upper diagonal only
-                    for (int i = 0; i < numKernels; i++) {
-                        for (int j = i; j < numKernels; j++) {
+                    for (int i = 0; i < bgIndex; i++) {
+                        for (int j = i; j < bgIndex; j++) {
                             stampMatrix->data.F64[i][j] += convolutions->data.F64[i] *
                                 convolutions->data.F64[j] * invNoise2;
@@ -674 +667 @@
             // Fill in lower diagonal of symmetric matrix, while checking for bad values
             bool bad = false;           // Are there bad values?
-            for (int i = 0; i < numKernels; i++) {
+            for (int i = 0; i < bgIndex; i++) {
                 for (int j = 0; j < i; j++) {
                     stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i];
@@ -694 +687 @@
             if (bad) {
                 stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-                psTrace("psModules.imcombine", 3, "Rejecting stamp %d because of bad equation\n", i);
+                psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n",
+                        i, stamp->x, stamp->y);
             } else {
                 stamp->status = PM_SUBTRACTION_STAMP_USED;
@@ -807 +801 @@
     PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1);
     PS_ASSERT_PTR_NON_NULL(kernels, -1);
-    PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, -1);
+    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
+                          polyTerms(kernels->bgOrder), -1);
+    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, -1);
     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, -1);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, -1);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, -1);
-    if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
-        PS_ASSERT_ARRAY_NON_EMPTY(kernels->preCalc, -1);
-        PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, -1);
-    }
 
     // Measure deviations
@@ -836 +826 @@
             int xStamp = stamp->x, yStamp = stamp->y; // Coordinates of stamp
 
-            // Precalulate polynomial values
-            psImage *polyValues = spatialPolyValues(kernels->spatialOrder,
-                                                    2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols,
-                                                    2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows);
+            float xNorm = 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols; // Normalised x coord
+            float yNorm = 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows; // Normalised y coord
+            psImage *polyValues = spatialPolyValues(kernels->spatialOrder, xNorm, yNorm); // Polynomial terms
 
 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
@@ -878 +867 @@
             }
             deviations->data.F32[i] = sqrt(stats->sampleMean / 2.0);
-            psTrace("psModules.imcombine", 1, "Deviation for stamp %d: %f\n", i, deviations->data.F32[i]);
+            psTrace("psModules.imcombine", 1, "Deviation for stamp %d (%d,%d): %f\n",
+                    i, stamp->x, stamp->y, deviations->data.F32[i]);
             totalSquareDev += PS_SQR(deviations->data.F32[i]);
             numStamps++;
@@ -896 +886 @@
         if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) {
             // Mask out the stamp in the image so you it's not found again
-            psTrace("psModules.imcombine", 3, "Rejecting stamp %d\n", i);
+            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, stamp->x, stamp->y);
             numRejected++;
             for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
@@ -922 +912 @@
     PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
     PS_ASSERT_PTR_NON_NULL(kernels, NULL);
+    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
+                          polyTerms(kernels->bgOrder), 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_VECTOR_SIZE(kernels->u, kernels->num, 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);
-    }
 
     // Precalulate polynomial values
@@ -958 +944 @@
     PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
     PS_ASSERT_PTR_NON_NULL(kernels, NULL);
+    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
+                          polyTerms(kernels->bgOrder), 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_VECTOR_SIZE(kernels->u, kernels->num, 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
@@ -987 +969 @@
 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
                            const psImage *inImage, const psImage *inWeight, const psImage *subMask,
-                           psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels)
+                           psMaskType blank, const psRegion *region, const psVector *solution,
+                           const pmSubtractionKernels *kernels)
 {
     PS_ASSERT_IMAGE_NON_NULL(inImage, false);
@@ -1022 +1005 @@
     PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false);
     PS_ASSERT_PTR_NON_NULL(kernels, false);
-    PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false);
+    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
+                          polyTerms(kernels->bgOrder), -1);
+    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false);
     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, false);
-    if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
-        PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false);
-        PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, false);
-    }
-
-    int size = kernels->size;           // Half-size of kernel
+    if (region) {
+        if (psRegionIsNaN(*region)) {
+            psString string = psRegionToString(*region);
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
+            psFree(string);
+            return false;
+        }
+        if (region->x0 < 0 || region->x1 > inImage->numCols ||
+            region->y0 < 0 || region->y1 > inImage->numRows) {
+            psString string = psRegionToString(*region);
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
+                    string, inImage->numCols, inImage->numRows);
+            psFree(string);
+            return false;
+        }
+    }
+
     float background = solution->data.F64[solution->n-1]; // The difference in background
     int numCols = inImage->numCols, numRows = inImage->numRows; // Image dimensions
@@ -1057 +1051 @@
     }
 
+    int size = kernels->size;           // Half-size of kernel
+    int fullSize = 2 * size + 1;        // Full size of kernel
     psImage *polyValues = NULL;         // Pre-calculated polynomial values
-    int fullSize = 2 * size + 1;        // Full size of kernel
 
     psKernel *kernelImage = NULL;       // Kernel for the image
     psKernel *kernelWeight = NULL;      // Kernel for the weight map
 
-    for (int j = size; j < inImage->numRows - size; j += fullSize) {
-        for (int i = size; i < inImage->numCols - size; i += fullSize) {
+    // Get region for convolution: [xMin:xMax,yMin:yMax]
+    int xMin = size, xMax = numCols - size;
+    int yMin = size, yMax = numRows - size;
+    if (region) {
+        xMin = PS_MAX(region->x0, xMin);
+        xMax = PS_MIN(region->x1, xMax);
+        yMin = PS_MAX(region->y0, yMin);
+        yMax = PS_MIN(region->y1, yMax);
+    }
+
+    for (int j = yMin; j < yMax; j += fullSize) {
+        for (int i = xMin; i < xMax; i += fullSize) {
 
             // Only generate polynomial values every kernel footprint, since we have already assumed
@@ -1085 +1090 @@
 #endif
 
-            for (int y = j; y < PS_MIN(j + fullSize, numRows - size); y++) {
-                for (int x = i; x < PS_MIN(i + fullSize, numCols - size); x++) {
+            for (int y = j; y < PS_MIN(j + fullSize, yMax); y++) {
+                for (int x = i; x < PS_MIN(i + fullSize, xMax); x++) {
                     // Check and propagate the kernel footprint, if required
                     if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &
@@ -1123 +1128 @@
     psFree(polyValues);
 
-    // Mark the non-convolved part as blank
-    for (int y = size; y < inImage->numRows - size; y++) {
+    return true;
+}
+
+bool pmSubtractionBorder(psImage *image, psImage *weight, psImage *mask,
+                         const pmSubtractionKernels *kernels, psMaskType blank)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, false);
+    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
+    if (mask) {
+        PS_ASSERT_IMAGE_NON_NULL(mask, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false);
+        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false);
+    }
+    if (weight) {
+        PS_ASSERT_IMAGE_NON_NULL(weight, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image, false);
+        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
+    }
+
+    int numCols = image->numCols, numRows = image->numRows; // Image dimensions
+
+    int size = kernels->size;           // Half-size of kernel
+    for (int y = size; y < numRows - size; y++) {
         for (int x = 0; x < size; x++) {
-            markBlank(convImage, convMask, convWeight, x, y, blank);
+            markBlank(image, mask, weight, x, y, blank);
         }
         for (int x = numCols - size; x < numCols; x++) {
-            markBlank(convImage, convMask, convWeight, x, y, blank);
+            markBlank(image, mask, weight, x, y, blank);
         }
     }
     for (int y = 0; y < size; y++) {
         for (int x = 0; x < numCols; x++) {
-            markBlank(convImage, convMask, convWeight, x, y, blank);
+            markBlank(image, mask, weight, x, y, blank);
         }
     }
     for (int y = numRows - size; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            markBlank(convImage, convMask, convWeight, x, y, blank);
+            markBlank(image, mask, weight, x, y, blank);
         }
     }