Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/imcombine/pmSubtraction.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtraction.c (modified) (51 diffs)

branches/simtest_nebulous_branches/psModules/src/imcombine/pmSubtraction.c

@@ -17 +17 @@
 #include <pslib.h>
 
+#include "pmErrorCodes.h"
 #include "pmHDU.h"                      // Required for pmFPA.h
 #include "pmFPA.h"
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionEquation.h"
+#include "pmSubtractionVisual.h"
 #include "pmSubtractionThreads.h"
 
@@ -29 +31 @@
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
 #define MIN_SAMPLE_STATS    7           // Minimum number to use sample statistics; otherwise use quartiles
+#define USE_KERNEL_ERR                  // Use kernel error image?
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -36 +39 @@
 // Generate the kernel to apply to the variance from the normal kernel
 static psKernel *varianceKernel(psKernel *out, // Output kernel
-                                psKernel *normalKernel // Normal kernel
+                                const psKernel *normalKernel // Normal kernel
                                 )
 {
@@ -48 +51 @@
 
     // Take the square of the normal kernel
-    double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels
+    double sumVariance = 0.0; // Sum of the variance kernels
     for (int v = yMin; v <= yMax; v++) {
         for (int u = xMin; u <= xMax; u++) {
-            float value = normalKernel->kernel[v][u]; // Value of interest
-            float value2 = PS_SQR(value); // Value squared
-            sumNormal += value;
-            sumVariance += value2;
-            out->kernel[v][u] = value2;
+            sumVariance += out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
         }
     }
@@ -65 +64 @@
     return out;
 }
+
+// Contribute to an image of the solved kernel component using the preCalculated image
+static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated
+                                const pmSubtractionKernels *kernels, // Kernel basis functions
+                                float value,                         // Normalisation value for basis function
+                                int index                  // Index of basis function of interest
+    )
+{
+    int size = kernels->size;           // Kernel half-size
+    pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values
+#if 0
+    // Iterating over the kernel
+    for (int y = 0, v = -size; v <= size; y++, v++) {
+        float yValue = value * preCalc->yKernel->data.F32[y];
+        for (int x = 0, u = -size; u <= size; x++, u++) {
+            kernel->kernel[v][u] +=  yValue * preCalc->xKernel->data.F32[x];
+        }
+    }
+    // Photometric scaling for even kernels only
+    if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
+        kernel->kernel[0][0] -= value;
+    }
+#else
+    for (int v = -size; v <= size; v++) {
+        for (int u = -size; u <= size; u++) {
+            kernel->kernel[v][u] +=  value * preCalc->kernel->kernel[v][u];
+        }
+    }
+#endif
+
+    return;
+}
+
 
 // Generate an image of the solved kernel
@@ -70 +102 @@
                               const pmSubtractionKernels *kernels, // Kernel basis functions
                               const psImage *polyValues, // Spatial polynomial values
+                              bool normalise,            // Add normalisation?
                               bool wantDual // Want the dual (second) kernel?
                               )
@@ -85 +118 @@
     for (int i = 0; i < numKernels; i++) {
         double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value
+        if (wantDual) {
+            // The model is built with the dual convolution terms added, so to produce zero residual the
+            // equation results in negative coefficients which we must undo.
+            value *= -1.0;
+        }
 
         switch (kernels->type) {
@@ -115 +153 @@
           case PM_SUBTRACTION_KERNEL_GUNK: {
               if (i < kernels->inner) {
-                  // 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] += preCalc->kernel[v][u] * value;
-                          // Photometric scaling is built into the preCalc kernel --- no subtraction!
-                      }
-                  }
+                  solvedKernelPreCalc(kernel, kernels, value, i);
               } else {
                   // Using delta function
@@ -133 +163 @@
               break;
           }
-          case PM_SUBTRACTION_KERNEL_ISIS: {
-              psArray *preCalc = kernels->preCalc->data[i]; // Precalculated values
-              psVector *xKernel = preCalc->data[0]; // Kernel in x
-              psVector *yKernel = preCalc->data[1]; // Kernel in y
-              // Iterating over the kernel
-              for (int y = 0, v = -size; v <= size; y++, v++) {
-                  for (int x = 0, u = -size; u <= size; x++, u++) {
-                      kernel->kernel[v][u] +=  value * xKernel->data.F32[x] * yKernel->data.F32[y];
-                  }
-              }
-              // Photometric scaling for even kernels only
-              if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
-                  kernel->kernel[0][0] -= value;
-              }
+          case PM_SUBTRACTION_KERNEL_ISIS:
+          case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
+          case PM_SUBTRACTION_KERNEL_HERM:
+          case PM_SUBTRACTION_KERNEL_DECONV_HERM: {
+              solvedKernelPreCalc(kernel, kernels, value, i);
               break;
           }
           case PM_SUBTRACTION_KERNEL_RINGS: {
-              psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data
-              psVector *uCoords = preCalc->data[0]; // u coordinates
-              psVector *vCoords = preCalc->data[1]; // v coordinates
-              psVector *poly = preCalc->data[2]; // Polynomial values
-              int num = uCoords->n;     // Number of pixels
+              pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Precalculated kernels
+              int num = preCalc->uCoords->n;     // Number of pixels
 
               for (int j = 0; j < num; j++) {
-                  int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
-                  kernel->kernel[v][u] += poly->data.F32[j] * value;
+                  int u = preCalc->uCoords->data.S32[j];
+                  int v = preCalc->vCoords->data.S32[j]; // Kernel coordinates
+                  kernel->kernel[v][u] += preCalc->poly->data.F32[j] * value;
               }
               // Photometric scaling is built into the kernel --- no subtraction!
@@ -168 +187 @@
     }
 
-    // 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;
@@ -250 +271 @@
 static void convolveVarianceFFT(psImage *target,// Place the result in here
                               psImage *variance, // Variance map to convolve
-                              psImage *sys, // Systematic error image
+                              psImage *kernelErr, // Kernel error image
                               psImage *mask, // Mask image
                               psImageMaskType maskVal, // Value to mask
@@ -262 +283 @@
 
     psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map
-    psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image
+    psImage *subKE = kernelErr ? psImageSubset(kernelErr, border) : NULL; // Kernel error image
     psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask
 
     // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
     psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
-    psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
+    psImage *convKE = subKE ? psImageConvolveFFT(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE
 
     psFree(subVariance);
-    psFree(subSys);
+    psFree(subKE);
     psFree(subMask);
 
     // Now, we have to stick it in where it belongs
     int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region
-    if (convSys) {
+    if (convKE) {
         for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
             for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) {
                 target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] +
-                    convSys->data.F32[ySource][xSource];
+                    convKE->data.F32[ySource][xSource];
             }
         }
@@ -290 +311 @@
 
     psFree(convVariance);
-    psFree(convSys);
+    psFree(convKE);
 
     return;
@@ -326 +347 @@
                                   psImage *image, // Image to convolve
                                   psImage *variance, // Variance map to convolve, or NULL
-                                  psImage *sys, // Systematic error image, or NULL
+                                  psImage *kernelErr, // Kernel error image, or NULL
                                   psImage *subMask, // Subtraction mask
                                   const pmSubtractionKernels *kernels, // Kernels
@@ -339 +360 @@
     )
 {
-    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
+    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, true, wantDual);
     if (variance || subMask) {
         *kernelVariance = varianceKernel(*kernelVariance, *kernelImage);
@@ -363 +384 @@
         convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size);
         if (variance) {
-            convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size);
+            convolveVarianceFFT(convVariance, variance, kernelErr, subMask, subBad, *kernelVariance,
+                                region, kernels->size);
         }
     } else {
@@ -373 +395 @@
     }
 
-    // Convolve the mask for bad pixels
+    // Convolve the mask for bad/poor pixels
     if (subMask && convMask) {
         int box = p_pmSubtractionBadRadius(*kernelImage, kernels, polyValues,
                                            wantDual, poorFrac); // Size of bad box
+        psAssert(box >= 0, "Bad radius must be >= 0");
+
+        int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
+        psImage *convolved = NULL; // Convolved subtraction mask
         if (box > 0) {
-            int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
-            psRegion region = psRegionSet(colMin - box, colMax + box,
-                                          rowMin - box, rowMax + box); // Region to convolve
-
-            psImage *image = subMask ? psImageSubset(subMask, region) : NULL; // Mask to convolve
-
-            psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad,
-                                                     -box, box, -box, box); // Convolved subtraction mask
-
+            psRegion maskRegion = psRegionSet(colMin - box, colMax + box,
+                                              rowMin - box, rowMax + box); // Region to convolve
+            psImage *image = subMask ? psImageSubset(subMask, maskRegion) : NULL; // Mask to convolve
+            convolved = psImageConvolveMask(NULL, image, subBad, subConvBad, -box, box, -box, box);
             psFree(image);
-
-            psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns");
-            psAssert(convolved->numRows - 2 * box == rowMax - rowMin, "Bad number of rows");
-
-            for (int yTarget = rowMin, ySource = box; yTarget < rowMax; yTarget++, ySource++) {
-                // Dereference images
-                psImageMaskType *target = &convMask->data.PS_TYPE_IMAGE_MASK_DATA[yTarget][colMin]; // Target values
-                psImageMaskType *source = &convolved->data.PS_TYPE_IMAGE_MASK_DATA[ySource][box]; // Source values
-                for (int xTarget = colMin; xTarget < colMax; xTarget++, target++, source++) {
-                    if (*source & subConvBad) {
-                        *target |= maskBad;
-                    } else if (*source & subConvPoor) {
-                        *target |= maskPoor;
-                    }
+        } else {
+            convolved = psImageSubset(subMask, region);
+        }
+
+        psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns");
+        psAssert(convolved->numRows - 2 * box == rowMax - rowMin, "Bad number of rows");
+
+        for (int yTarget = rowMin, ySource = box; yTarget < rowMax; yTarget++, ySource++) {
+            // Dereference images
+            psImageMaskType *target = &convMask->data.PS_TYPE_IMAGE_MASK_DATA[yTarget][colMin]; // Target values
+            psImageMaskType *source = &convolved->data.PS_TYPE_IMAGE_MASK_DATA[ySource][box]; // Source values
+            for (int xTarget = colMin; xTarget < colMax; xTarget++, target++, source++) {
+                if (*source & subConvBad) {
+                    *target |= maskBad;
+                } else if (*source & subConvPoor) {
+                    *target &= ~maskBad;
+                    *target |= maskPoor;
+                } else {
+                    *target &= ~maskBad & ~maskPoor;
                 }
             }
-
-            // No need to lock: we own this
-            psFree(convolved);
-        }
+        }
+
+        psFree(convolved);
     }
 
@@ -413 +438 @@
 }
 
-// Generate an image that can be used to track systematic errors
-static psImage *subtractionSysErrImage(const psImage *image, // Image from which to make sys err image
-                                       float sysError // Relative systematic error
-                                       )
-{
-    if (!isfinite(sysError) || sysError == 0.0) {
+#ifdef USE_KERNEL_ERR
+// Generate an image that can be used to track systematic errors in the kernel
+static psImage *subtractionKernelErrImage(const psImage *image, // Image from which to make kernel error image
+                                          float kernelError // Relative systematic error in kernel
+    )
+{
+    if (!isfinite(kernelError) || kernelError == 0.0) {
         return NULL;
     }
 
     int numCols = image->numCols, numRows = image->numRows; // Size of image
-    psImage *sys = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Systematic error image
-
-    float sysError2 = PS_SQR(sysError); // Square of the systematic error
+    psImage *kernelErr = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Kernel error image
+
+    float kernelError2 = PS_SQR(kernelError); // Square of the kernel error
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            sys->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * sysError2;
-        }
-    }
-
-    return sys;
+            kernelErr->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * kernelError2;
+        }
+    }
+
+    return kernelErr;
+}
+#endif
+
+// Convolve a stamp using a pre-calculated kernel basis function
+static psKernel *convolveStampPreCalc(const psKernel *image, // Image to convolve
+                                      const pmSubtractionKernels *kernels, // Kernel basis functions
+                                      int index,                            // Index of basis function of interest
+                                      int footprint                         // Half-size of stamp
+    )
+{
+    pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data
+#if 0
+    // Convolving using separable convolution
+    int size = kernels->size;     // Size of kernel
+
+    // Convolve in x
+    // Need to convolve a bit more than the footprint, for the y convolution
+    int yMin = -size - footprint, yMax = size + footprint; // Range for y
+    psKernel *temp = psKernelAlloc(yMin, yMax,
+                                   -footprint, footprint); // Temporary convolution; NOTE: wrong way!
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            float value = 0.0;    // Value of convolved pixel
+            int uMin = x - size, uMax = x + size; // Range for u
+            psF32 *xKernelData = &preCalc->xKernel->data.F32[xKernel->n - 1]; // Kernel values
+            psF32 *imageData = &image->kernel[y][uMin]; // Image values
+            for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
+                value += *xKernelData * *imageData;
+            }
+            temp->kernel[x][y] = value; // NOTE: putting in wrong way!
+        }
+    }
+
+    // Convolve in y
+    psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint);// Convolved image
+    for (int x = -footprint; x <= footprint; x++) {
+        for (int y = -footprint; y <= footprint; y++) {
+            float value = 0.0;    // Value of convolved pixel
+            int vMin = y - size, vMax = y + size; // Range for v
+            psF32 *yKernelData = &preCalc->yKernel->data.F32[yKernel->n - 1]; // Kernel values
+            psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
+            for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
+                value += *yKernelData * *imageData;
+            }
+            convolved->kernel[y][x] = value;
+        }
+    }
+    psFree(temp);
+
+    // Photometric scaling for even kernels only
+    if (kernels->u->data.S32[index] % 2 == 0 && kernels->v->data.S32[index] % 2 == 0) {
+        convolveSub(convolved, image, footprint);
+    }
+    return convolved;
+#else
+    // Convolving using precalculated kernel
+    return p_pmSubtractionConvolveStampPrecalc(image, preCalc->kernel);
+#endif
 }
 
@@ -447 +531 @@
     int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
     psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
+
+    // pmSubtractionVisualShowSubtraction(image->image, kernel->image, conv);
+
     psFree(conv);
     return convolved;
@@ -456 +543 @@
     psKernel *kernel;                   // Kernel to use
     if (!preKernel) {
-        kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
+        kernel = solvedKernel(NULL, kernels, polyValues, true, wantDual);
     } else {
         kernel = psMemIncrRefCounter(preKernel);
@@ -491 +578 @@
 }
 
+void p_pmSubtractionPolynomialNormCoords(float *xOut, float *yOut, float xIn, float yIn,
+                                         int xMin, int xMax, int yMin, int yMax)
+{
+    float xNormSize = xMax - xMin, yNormSize = yMax - yMin; // Size to use for normalisation
+    *xOut = 2.0 * (float)(xIn - xMin - xNormSize/2.0) / xNormSize;
+    *yOut = 2.0 * (float)(yIn - yMin - yNormSize/2.0) / yNormSize;
+    return;
+}
+
 psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels,
-                                             int numCols, int numRows, int x, int y)
+                                             int x, int y)
 {
     assert(kernels);
-    assert(numCols > 0 && numRows > 0);
-
-    // Size to use when calculating normalised coordinates (different from actual size when convolving
-    // subimage)
-    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
-    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
-
-    // Normalised coordinates
-    float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize;
-    float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize;
-
+
+    float xNorm, yNorm;                 // Normalised coordinates
+    p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, x, y,
+                                        kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax);
     return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm);
 }
@@ -586 +675 @@
           if (index < kernels->inner) {
               // Photometric scaling is already built in to the precalculated kernel
-              return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]);
+              return convolveStampPreCalc(image, kernels, index, footprint);
           }
           // Using delta function
@@ -595 +684 @@
           return convolved;
       }
-      case PM_SUBTRACTION_KERNEL_ISIS: {
-          psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values
-          psVector *xKernel = preCalc->data[0]; // Kernel in x
-          psVector *yKernel = preCalc->data[1]; // Kernel in y
-          int size = kernels->size;     // Size of kernel
-
-          // Convolve in x
-          // Need to convolve a bit more than the footprint, for the y convolution
-          int yMin = -size - footprint, yMax = size + footprint; // Range for y
-          psKernel *temp = psKernelAlloc(yMin, yMax,
-                                         -footprint, footprint); // Temporary convolution; NOTE: wrong way!
-          for (int y = yMin; y <= yMax; y++) {
-              for (int x = -footprint; x <= footprint; x++) {
-                  float value = 0.0;    // Value of convolved pixel
-                  int uMin = x - size, uMax = x + size; // Range for u
-                  psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values
-                  psF32 *imageData = &image->kernel[y][uMin]; // Image values
-                  for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
-                      value += *xKernelData * *imageData;
-                  }
-                  temp->kernel[x][y] = value; // NOTE: putting in wrong way!
-              }
-          }
-
-          // Convolve in y
-          psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint);// Convolved image
-          for (int x = -footprint; x <= footprint; x++) {
-              for (int y = -footprint; y <= footprint; y++) {
-                  float value = 0.0;    // Value of convolved pixel
-                  int vMin = y - size, vMax = y + size; // Range for v
-                  psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values
-                  psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
-                  for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
-                      value += *yKernelData * *imageData;
-                  }
-                  convolved->kernel[y][x] = value;
-              }
-          }
-          psFree(temp);
-
-          // Photometric scaling for even kernels only
-          if (kernels->u->data.S32[index] % 2 == 0 && kernels->v->data.S32[index] % 2 == 0) {
-              convolveSub(convolved, image, footprint);
-          }
-          return convolved;
-      }
+      case PM_SUBTRACTION_KERNEL_ISIS:
+      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
+      case PM_SUBTRACTION_KERNEL_HERM:
+      case PM_SUBTRACTION_KERNEL_DECONV_HERM: {
+            return convolveStampPreCalc(image, kernels, index, footprint);
+        }
       case PM_SUBTRACTION_KERNEL_RINGS: {
-          psKernel *convolved = psKernelAlloc(-footprint, footprint,
-                                              -footprint, footprint); // Convolved image
-          psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
-          psVector *uCoords = preCalc->data[0]; // u coordinates
-          psVector *vCoords = preCalc->data[1]; // v coordinates
-          psVector *poly = preCalc->data[2]; // Polynomial values
-          int num = uCoords->n;         // Number of pixels
-          psS32 *uData = uCoords->data.S32, *vData = vCoords->data.S32; // Dereference u,v coordinates
-          psF32 *polyData = poly->data.F32; // Dereference polynomial values
+          psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
+          pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data
+
+          int num = preCalc->uCoords->n;         // Number of pixels
+          psS32 *uData = preCalc->uCoords->data.S32; // Dereference v coordinate
+          psS32 *vData = preCalc->vCoords->data.S32; // Dereference u coordinate
+          psF32 *polyData = preCalc->poly->data.F32; // Dereference polynomial values
           psF32 **imageData = image->kernel;  // Dereference image
           psF32 **convData = convolved->kernel; // Dereference convolved image
@@ -706 +753 @@
 
     if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Stamp not marked for calculation.");
+        psError(PM_ERR_PROG, true, "Stamp not marked for calculation.");
         return false;
     }
@@ -730 +777 @@
 
 
-
-
 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
-                              const psVector *deviations, psImage *subMask, float sigmaRej, int footprint)
+                              const psVector *deviations, psImage *subMask, float sigmaRej)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
@@ -764 +809 @@
 
     if (numStamps == 0) {
-        psError(PS_ERR_UNKNOWN, true, "No good stamps found.");
+        psError(PM_ERR_STAMPS, true, "No good stamps found.");
         psFree(mask);
         return -1;
@@ -772 +817 @@
                                   PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics for deviatns
     if (!psVectorStats(stats, deviations, NULL, mask, 0xff)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for deviations.");
+        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
         psFree(stats);
         psFree(mask);
@@ -821 +866 @@
     ds9num++;
 
+    int footprint = stamps->footprint;  // Half-size of stamp region of interest
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
     double newMean = 0.0;               // New mean
+    psString log = NULL;                // Log message
+    psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit);
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -836 +884 @@
                 // Mask out the stamp in the image so you it's not found again
                 psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
-                        (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+                        (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
+                psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i,
+                               (int)(stamp->x - 0.5), (int)(stamp->y - 0.5),
+                               fabsf(deviations->data.F32[i] - mean));
                 numRejected++;
                 for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
@@ -857 +908 @@
                 psFree(stamp->image1);
                 psFree(stamp->image2);
-                psFree(stamp->variance);
-                stamp->image1 = stamp->image2 = stamp->variance = NULL;
-                psFree(stamp->matrix1);
-                psFree(stamp->matrix2);
-                psFree(stamp->matrixX);
-                stamp->matrix1 = stamp->matrix2 = stamp->matrixX = NULL;
-                psFree(stamp->vector1);
-                psFree(stamp->vector2);
-                stamp->vector1 = stamp->vector2 = NULL;
+                psFree(stamp->weight);
+                stamp->image1 = stamp->image2 = stamp->weight = NULL;
+                psFree(stamp->matrix);
+                stamp->matrix = NULL;
+                psFree(stamp->vector);
+                stamp->vector = NULL;
             } else {
                 numGood++;
@@ -874 +922 @@
     }
     newMean /= numGood;
+
+    if (numRejected == 0) {
+        psStringAppend(&log, "<none>\n");
+    }
+    psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log);
+    psFree(log);
 
     if (ds9) {
@@ -900 +954 @@
 
     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);
 
@@ -931 +985 @@
     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);
 
@@ -948 +1002 @@
 }
 
-#if 0
+#if 1
 psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
 {
@@ -956 +1010 @@
     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 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(kernels, x, y, wantDual);
-        fakeSolution->data.F64[i] = 0.0;
-    }
-
-    psFree(fakeSolution);
+    psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution of interest
+    psVector *backup = psVectorCopy(NULL, solution, PS_TYPE_F64);  // Backup version
+
+    int num = kernels->num;             // Number of kernel basis functions
+
+    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];
+        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);
 
     return images;
@@ -979 +1059 @@
                                      psImage *convMask, // Output convolved mask
                                      const pmReadout *ro1, const pmReadout *ro2, // Input readouts
-                                     psImage *sys1, psImage *sys2, // Systematic error images
+                                     psImage *kernelErr1, psImage *kernelErr2, // Kernel error images
                                      psImage *subMask, // Input subtraction mask
                                      psImageMaskType maskBad, // Mask value to give bad pixels
@@ -998 +1078 @@
     // Only generate polynomial values every kernel footprint, since we have already assumed
     // (with the stamps) that it does not vary rapidly on this scale.
-    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
-                                                              xMin + x0 + size + 1,
-                                                              yMin + y0 + size + 1);
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + x0 + size + 1,
+                                                              yMin + y0 + size + 1);        // Polynomial
     float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term
 
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance,
-                       ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region,
-                       maskBad, maskPoor, poorFrac, useFFT, false);
+        convolveRegion(out1->image, out1->variance, out1->mask, &kernelImage, &kernelVariance,
+                       ro1->image, ro1->variance, kernelErr1, subMask, kernels, polyValues, background,
+                       *region, maskBad, maskPoor, poorFrac, useFFT, false);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance,
-                       ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region,
-                       maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
+        convolveRegion(out2->image, out2->variance, out2->mask, &kernelImage, &kernelVariance,
+                       ro2->image, ro2->variance, kernelErr2, subMask, kernels, polyValues, background,
+                       *region, maskBad, maskPoor, poorFrac, useFFT,
+                       kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
 
@@ -1019 +1099 @@
 
     if ((kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) && ro1->mask) {
-        psImageMaskType **target = convMask->data.PS_TYPE_IMAGE_MASK_DATA; // Target mask
+        psImageMaskType **target = out1->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Target mask
         psImageMaskType **source = ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Source mask
 
@@ -1029 +1109 @@
     }
     if ((kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) && ro2->mask) {
-        psImageMaskType **target = convMask->data.PS_TYPE_IMAGE_MASK_DATA; // Target mask
+        psImageMaskType **target = out2->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Target mask
         psImageMaskType **source = ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Source mask
 
@@ -1056 +1136 @@
     const pmReadout *ro1 = args->data[7]; // Input readout 1
     const pmReadout *ro2 = args->data[8]; // Input readout 2
-    psImage *sys1 = args->data[9]; // Systematic error image 1
-    psImage *sys2 = args->data[10]; // Systematic error image 2
+    psImage *kernelErr1 = args->data[9]; // Kernel error image 1
+    psImage *kernelErr2 = args->data[10]; // Kernel error image 2
     psImage *subMask = args->data[11]; // Subtraction mask
     psImageMaskType maskBad = PS_SCALAR_VALUE(args->data[12], PS_TYPE_IMAGE_MASK_DATA); // Output mask value for bad pixels
@@ -1067 +1147 @@
     bool useFFT = PS_SCALAR_VALUE(args->data[18], U8); // Use FFT for convolution?
 
-    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
-                                    subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
+    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, kernelErr1,
+                                    kernelErr2, subMask, maskBad, maskPoor, poorFrac, region, kernels,
+                                    doBG, useFFT);
 }
 
 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2,
                            psImage *subMask, int stride, psImageMaskType maskBad, psImageMaskType maskPoor,
-                           float poorFrac, float sysError, const psRegion *region,
+                           float poorFrac, float kernelError, float covarFrac, const psRegion *region,
                            const pmSubtractionKernels *kernels, bool doBG, bool useFFT)
 {
@@ -1082 +1163 @@
         PM_ASSERT_READOUT_NON_NULL(ro1, false);
         PM_ASSERT_READOUT_IMAGE(ro1, false);
+        PM_ASSERT_READOUT_IMAGE(out1, false);
         numCols = ro1->image->numCols;
         numRows = ro1->image->numRows;
@@ -1091 +1173 @@
         PM_ASSERT_READOUT_NON_NULL(ro2, false);
         PM_ASSERT_READOUT_IMAGE(ro2, false);
+        PM_ASSERT_READOUT_IMAGE(out2, false);
         if (numCols == 0 && numRows == 0) {
             numCols = ro2->image->numCols;
@@ -1111 +1194 @@
     PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(poorFrac, 0.0, false);
     PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, false);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
-    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(sysError, 1.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(kernelError, 1.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(covarFrac, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN(covarFrac, 1.0, false);
     if (region && psRegionIsNaN(*region)) {
         psString string = psRegionToString(*region);
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
+        psError(PM_ERR_PROG, true, "Input region (%s) contains NAN values", string);
         psFree(string);
         return false;
@@ -1124 +1209 @@
     bool threaded = pmSubtractionThreaded(); // Running threaded?
 
-    // Outputs
-    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        if (!out1->image) {
-            out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            // XXX if (threaded) {
-            // XXX     psMutexInit(out1->image);
-            // XXX }
-        }
-        if (ro1->variance) {
-            if (!out1->variance) {
-                out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-                // XXX if (threaded) {
-                // XXX     psMutexInit(out1->variance);
-                // XXX }
-            }
-            psImageInit(out1->variance, 0.0);
-        }
-    }
-    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        if (!out2->image) {
-            out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            // XXX if (threaded) {
-            // XXX     psMutexInit(out2->image);
-            // XXX }
-        }
-        if (ro2->variance) {
-            if (!out2->variance) {
-                out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-                // XXX if (threaded) {
-                // XXX     psMutexInit(out2->variance);
-                // XXX }
-            }
-            psImageInit(out2->variance, 0.0);
-        }
-    }
     psImage *convMask = NULL;           // Convolved mask image (common to inputs 1 and 2)
     if (subMask) {
-        // XXX if (threaded) {
-        // XXX     psMutexInit(subMask);
-        // XXX }
         if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            if (!out1->mask) {
-                out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            }
             convMask = out1->mask;
         }
         if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            if (convMask) {
-                if (out2->mask) {
-                    psFree(out2->mask);
-                }
-                out2->mask = psMemIncrRefCounter(convMask);
-            } else {
-                if (!out2->mask) {
-                    out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-                }
+            if (!convMask) {
                 convMask = out2->mask;
             }
         }
-        psImageInit(convMask, 0);
-    }
-
-    psImage *sys1 = NULL, *sys2 = NULL; // Systematic error images
+    }
+
+    psImage *kernelErr1 = NULL, *kernelErr2 = NULL; // Kernel error images
+#ifdef USE_KERNEL_ERR
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        sys1 = subtractionSysErrImage(ro1->image, sysError);
-        // XXX if (threaded && sys1) {
-        // XXX     psMutexInit(sys1);
-        // XXX }
+        kernelErr1 = subtractionKernelErrImage(ro1->image, kernelError);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        sys2 = subtractionSysErrImage(ro2->image, sysError);
-        // XXX if (threaded && sys2) {
-        // XXX     psMutexInit(sys2);
-        // XXX }
-    }
+        kernelErr2 = subtractionKernelErrImage(ro2->image, kernelError);
+    }
+#endif
 
     int size = kernels->size;           // Half-size of kernel
 
     // Get region for convolution: [xMin:xMax,yMin:yMax]
-    int xMin = size, xMax = numCols - size;
-    int yMin = size, yMax = numRows - size;
+    int xMin = kernels->xMin + size, xMax = kernels->xMax - size;
+    int yMin = kernels->yMin + size, yMax = kernels->yMax - size;
     if (region) {
         xMin = PS_MAX(region->x0, xMin);
@@ -1247 +1278 @@
                 psArrayAdd(args, 1, (pmReadout*)ro1); // Casting away const
                 psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const
-                // Since adding to the array can impact the reference count, we need to lock
-                // XXX if (sys1) {
-                // XXX     psMutexLock(sys1);
-                // XXX }
-                psArrayAdd(args, 1, sys1);
-                // XXX if (sys1) {
-                // XXX     psMutexUnlock(sys1);
-                // XXX }
-                // XXX if (sys2) {
-                // XXX     psMutexLock(sys2);
-                // XXX }
-                psArrayAdd(args, 1, sys2);
-                // XXX if (sys2) {
-                // XXX     psMutexUnlock(sys2);
-                // XXX }
-                // XXX if (subMask) {
-                // XXX     psMutexLock(subMask);
-                // XXX }
+                psArrayAdd(args, 1, kernelErr1);
+                psArrayAdd(args, 1, kernelErr2);
                 psArrayAdd(args, 1, subMask);
-                // XXX if (subMask) {
-                // XXX     psMutexUnlock(subMask);
-                // XXX }
                 PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_IMAGE_MASK);
                 PS_ARRAY_ADD_SCALAR(args, maskPoor, PS_TYPE_IMAGE_MASK);
@@ -1283 +1295 @@
                 psFree(job);
             } else {
-                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
-                                         subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG,
-                                         useFFT);
+                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,
+                                         kernelErr1, kernelErr2, subMask, maskBad, maskPoor, poorFrac,
+                                         subRegion, kernels, doBG, useFFT);
             }
             psFree(subRegion);
@@ -1292 +1304 @@
 
     if (!psThreadPoolWait(false)) {
-        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        psError(psErrorCodeLast(), false, "Error waiting for threads.");
         return false;
     }
@@ -1307 +1319 @@
             psFree(job);
         }
-
-        // XXX if (subMask) {
-        // XXX     psMutexDestroy(subMask);
-        // XXX }
-        // XXX if (sys1) {
-        // XXX     psMutexDestroy(sys1);
-        // XXX }
-        // XXX if (sys2) {
-        // XXX     psMutexDestroy(sys2);
-        // XXX }
     }
     psImageConvolveSetThreads(oldThreads);
 
-    psFree(sys1);
-    psFree(sys2);
+    psFree(kernelErr1);
+    psFree(kernelErr2);
 
     // Calculate covariances
-    // This can take a while, so we only do it for a single instance
-    // XXX psImageCovarianceCalculate could be multithreaded
+    // This can be fairly involved, so we only do it for a single instance
+    // Enable threads for covariance calculation, since we're not threading on top of it.
+    oldThreads = psImageCovarianceSetThreads(true);
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+        psKernelTruncate(kernel, covarFrac);
         out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance);
         psFree(kernel);
@@ -1334 +1338 @@
         psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0,
                                                kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel
+        psKernelTruncate(kernel, covarFrac);
         out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance);
         psFree(kernel);
     }
+    psImageCovarianceSetThreads(oldThreads);
 
     // Copy anything that wasn't convolved
@@ -1376 +1382 @@
     }
 
-
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve image: %f sec",
              psTimerClear("pmSubtractionConvolve"));