Changeset 26893 for trunk/psModules/src/imcombine

trunk/psModules/src/imcombine/Makefile.am

-              r23242
+              r26893
         pmSubtractionIO.c       \
         pmSubtractionKernels.c  \
+        pmSubtractionHermitian.c        \
+        pmSubtractionDeconvolve.c       \
         pmSubtractionMask.c     \
         pmSubtractionMatch.c    \
 …
         pmSubtractionIO.h       \
         pmSubtractionKernels.h  \
+        pmSubtractionHermitian.h        \
+        pmSubtractionDeconvolve.h       \
         pmSubtractionMask.h     \
         pmSubtractionMatch.h    \

trunk/psModules/src/imcombine/pmImageCombine.c

-              r23989
+              r26893
             psImage *mask = masks->data[i]; // Mask of interest
             pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal);
+        }
         // Set the pixel error data, if necessary
+        }
+        // Set the pixel error data, if necessary
         if (errors) {
             psImage *error = errors->data[i]; // Error image of interest
 …
         // Combine all the pixels, using the specified stat.
         if (!psVectorStats(stats, pixelData, pixelErrors, pixelMasks, 0xff)) {
             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
             return false;
+        }
         if (isnan(stats->sampleMean)) {
+            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+            return false;
+        }
+        if (isnan(stats->sampleMean)) {
             combine->data.F32[y][x] = NAN;
             psFree(buffer);
 …
+        }
         float combinedPixel = stats->sampleMean; // Value of the combination
         if (iter == 0) {
             // Save the value produced with no rejection, since it may be useful later
 …
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     if (!psVectorStats(stats, pixels, NULL, mask, 1)) {
         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+    }
     float median = stats->sampleMedian;
 …
                                               (psPlaneTransform * )outToIn->data[otherImg],
                                               outCoords);
                         psS32 xPix = (int)(inCoords->x + 0.5);
                         psS32 yPix = (int)(inCoords->y + 0.5);
+                        psS32 xPix = (int)(inCoords->x - 0.5);
+                        psS32 yPix = (int)(inCoords->y - 0.5);
                         if ((xPix >= 0) && (xPix <= ((psImage*)(images->data[otherImg]))->numCols - 1) &&
                                 (yPix >= 0) && (yPix <= ((psImage*)(images->data[otherImg]))->numRows - 1)) {

trunk/psModules/src/imcombine/pmPSFEnvelope.c

-              r25754
+              r26893
                      int radius,        // Radius of each PSF
                      const char *modelName,// Name of PSF model to use
+                     int xOrder, int yOrder // Order for PSF variation fit
+                     int xOrder, int yOrder, // Order for PSF variation fit
+                     psImageMaskType maskVal
+                     )
+{
 …
             pmModelClassSetLimits(PM_MODEL_LIMITS_IGNORE);
             pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, PEAK_FLUX); // Test model
+            model->modelSetLimits(PM_MODEL_LIMITS_STRICT);
+            for (int j = 0; j < model->params->n && goodPSF; j++) {
+                if (!model->modelLimits(PS_MINIMIZE_PARAM_MIN, j, model->params->data.F32, NULL) ||
+                    !model->modelLimits(PS_MINIMIZE_PARAM_MAX, j, model->params->data.F32, NULL)) {
+                    goodPSF = false;
+            if (!model) {
+                goodPSF = false;
+            } else {
+                model->modelSetLimits(PM_MODEL_LIMITS_STRICT);
+                for (int j = 0; j < model->params->n && goodPSF; j++) {
+                    if (!model->modelLimits(PS_MINIMIZE_PARAM_MIN, j, model->params->data.F32, NULL) ||
+                        !model->modelLimits(PS_MINIMIZE_PARAM_MAX, j, model->params->data.F32, NULL)) {
+                        goodPSF = false;
+                    }
+                }
+            }
             psFree(model);
+                psFree(model);
+            }
             if (!goodPSF) {
                 psWarning("PSF %d is bad --- not including in envelope calculation.", i);
 …
         // measure the source moments: tophat windowing, no pixel S/N cutoff
+        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) {
+        // XXX probably should be passing the maskVal to this function so we can pass it along here...
+        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0, maskVal)) {
             // Can't do anything about it; limp along as best we can
             psErrorClear();

trunk/psModules/src/imcombine/pmPSFEnvelope.h

-              r15837
+              r26893
                      int radius,        // Radius of each PSF
                      const char *modelName, // Name of PSF model to use
+                     int xOrder, int yOrder // Order for PSF variation
+                     int xOrder, int yOrder, // Order for PSF variation
+                     psImageMaskType maskVal
     );

trunk/psModules/src/imcombine/pmStackReject.c

-              r25468
+              r26893
+{
     int size = kernels->size;           // Half-size of convolution kernel
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
+                                                              xMin + size + 1, yMin + size + 1); // Polynomial
+    int x = PS_MIN(xMin + size + 1, kernels->xMax); // x coordinate of interest
+    int y = PS_MIN(yMin + size + 1, kernels->yMax); // y coordinate of interest
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, x, y); // Polynomial
     int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues, false, poorFrac); // Radius of bad box
     psTrace("psModules.imcombine", 10, "Growing by %d", box);
 …
     pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
     inRO->image = image;
+    convRO->image = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = subRegions->data[i]; // Region of interest
         pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest
         if (!pmSubtractionConvolve(NULL, convRO, NULL, inRO, NULL, stride, 0, 0, 1.0, 0.0,
+        if (!pmSubtractionConvolve(NULL, convRO, NULL, inRO, NULL, stride, 0, 0, 1.0, 0.0, 0.0,
                                    region, kernels, false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
 …
         // Image of the kernel at the centre of the region
+        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) /
+            (float)kernels->numCols;
+        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) /
+            (float)kernels->numRows;
+        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
+        psImage *kernel = pmSubtractionKernelImage(kernels, 0.5, 0.5, false);
         if (!kernel) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");

trunk/psModules/src/imcombine/pmSubtraction.c

-              r26584
+              r26893
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionEquation.h"
+#include "pmSubtractionVisual.h"
 #include "pmSubtractionThreads.h"
 …
+}
 // Contribute to an image of the solved kernel component for ISIS
 static void solvedKernelISIS(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
+// 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
     psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values
+    pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values
 #if 0
-    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++) {
         float yValue = value * yKernel->data.F32[y];
+        float yValue = value * preCalc->yKernel->data.F32[y];
         for (int x = 0, u = -size; u <= size; x++, u++) {
             kernel->kernel[v][u] +=  yValue * xKernel->data.F32[x];
+            kernel->kernel[v][u] +=  yValue * preCalc->xKernel->data.F32[x];
+        }
+    }
 …
+    }
 #else
-    psKernel *k = preCalc->data[2]; // Kernel image
     for (int v = -size; v <= size; v++) {
         for (int u = -size; u <= size; u++) {
             kernel->kernel[v][u] +=  value * k->kernel[v][u];
+            kernel->kernel[v][u] +=  value * preCalc->kernel->kernel[v][u];
+        }
+    }
 …
     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) {
 …
           case PM_SUBTRACTION_KERNEL_GUNK: {
               if (i < kernels->inner) {
                   solvedKernelISIS(kernel, kernels, value, i);
+                  solvedKernelPreCalc(kernel, kernels, value, i);
               } else {
                   // Using delta function
 …
               break;
+          }
+          case PM_SUBTRACTION_KERNEL_ISIS: {
+              solvedKernelISIS(kernel, kernels, value, i);
+          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!
 …
                     *target |= maskBad;
                 } else if (*source & subConvPoor) {
+                    *target &= ~maskBad;
                     *target |= maskPoor;
+                } else {
+                    *target &= ~maskBad & ~maskPoor;
+                }
+            }
 …
 #endif
 // Convolve a stamp using an ISIS kernel basis function
 static psKernel *convolveStampISIS(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
+// 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
+    )
+{
     psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
+    pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data
 #if 0
     // Convolving using separable convolution
-    psVector *xKernel = preCalc->data[0]; // Kernel in x
-    psVector *yKernel = preCalc->data[1]; // Kernel in y
     int size = kernels->size;     // Size of kernel
 …
             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 *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++) {
 …
             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 *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++) {
 …
 #else
     // Convolving using precalculated kernel
     return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]);
+    return p_pmSubtractionConvolveStampPrecalc(image, preCalc->kernel);
 #endif
+}
 …
     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;
 …
+}
+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);
+}
 …
           if (index < kernels->inner) {
               // Photometric scaling is already built in to the precalculated kernel
               return convolveStampISIS(image, kernels, index, footprint);
+              return convolveStampPreCalc(image, kernels, index, footprint);
+          }
           // Using delta function
 …
           return convolved;
+      }
+      case PM_SUBTRACTION_KERNEL_ISIS: {
+          return convolveStampISIS(image, kernels, index, footprint);
+      }
+      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
 …
 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
                               const psVector *deviations, psImage *subMask, float sigmaRej)
 …
     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
 …
                 // 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++) {
 …
                 psFree(stamp->weight);
                 stamp->image1 = stamp->image2 = stamp->weight = 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->matrix);
+                stamp->matrix = NULL;
+                psFree(stamp->vector);
+                stamp->vector = NULL;
             } else {
                 numGood++;
 …
+    }
     newMean /= numGood;
+    if (numRejected == 0) {
+        psStringAppend(&log, "<none>\n");
+    }
+    psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log);
+    psFree(log);
     if (ds9) {
 …
     psVector *backup = psVectorCopy(NULL, solution, PS_TYPE_F64);  // Backup version
     int num = wantDual ? solution->n - 1 : solution->n; // Number of kernel basis functions
+    int num = kernels->num;             // Number of kernel basis functions
     psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial
 …
     // 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
 …
 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2,
                            psImage *subMask, int stride, psImageMaskType maskBad, psImageMaskType maskPoor,
                            float poorFrac, float kernelError, const psRegion *region,
+                           float poorFrac, float kernelError, float covarFrac, const psRegion *region,
                            const pmSubtractionKernels *kernels, bool doBG, bool useFFT)
+{
 …
         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;
 …
         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;
 …
     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);
 …
     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);
+        }
-        if (ro1->variance) {
-            if (!out1->variance) {
-                out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
-            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);
+        }
-        if (ro2->variance) {
-            if (!out2->variance) {
-                out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
-            psImageInit(out2->variance, 0.0);
+        }
+    }
     psImage *convMask = NULL;           // Convolved mask image (common to inputs 1 and 2)
     if (subMask) {
         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);
+            }
-            psImageInit(out1->mask, 0);
             convMask = out1->mask;
+        }
         if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            if (!out2->mask) {
-                out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            }
-            psImageInit(out2->mask, 0);
             if (!convMask) {
                 convMask = out2->mask;
 …
     // 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);
 …
     // 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);
 …
         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
 …
+    }
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve image: %f sec",
              psTimerClear("pmSubtractionConvolve"));

trunk/psModules/src/imcombine/pmSubtraction.h

-              r25279
+              r26893
                            psImageMaskType maskPoor, ///< Mask value to give poor pixels
                            float poorFrac, ///< Fraction for "poor"
+                           float sysError, ///< Relative systematic error
+                           float kernelError, ///< Relative systematic error in kernel
+                           float covarFrac,  ///< Truncation fraction for kernel before covariance calculation
                            const psRegion *region, ///< Region to convolve (or NULL)
                            const pmSubtractionKernels *kernels, ///< Kernel parameters
 …
     );
+/// Return normalised coordinates
+void p_pmSubtractionPolynomialNormCoords(
+    float *xOut, float *yOut,           ///< Normalised coordinates, returned
+    float xIn, float yIn,               ///< Input coordinates
+    int xMin, int xMax, int yMin, int yMax ///< Bounds of validity
+    );
 /// Given (normalised) coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
 psImage *p_pmSubtractionPolynomial(psImage *output, ///< Output matrix, or NULL
 …
 psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, ///< Output matrix, or NULL
                                              const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                             int numCols, int numRows, ///< Size of image of interest
                                              int x, int y ///< Position of interest
     );

trunk/psModules/src/imcombine/pmSubtractionAnalysis.c

-              r25999
+              r26893
+    }
+    // sample difference images
+    {
+        psMetadataAddArray(analysis, PS_LIST_TAIL, "SUBTRACTION.SAMPLE.STAMP.SET", PS_META_DUPLICATE_OK, "Sample Difference Stamps", kernels->sampleStamps);
+    }
 #ifdef TESTING
 …
     // Difference in background
+    {
+        psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
+                                                                  numCols / 2.0, numRows / 2.0); // Polynomial
+        psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0); // Polynomial
         float bg = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background difference
 …
         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
                          kernels->rms);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
+        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
+    }

trunk/psModules/src/imcombine/pmSubtractionAnalysis.h

-              r25060
+              r26893
 #define PM_SUBTRACTION_ANALYSIS_DECONV_MAX   "SUBTRACTION.DECONV.MAX"   // Maximum deconvolution fraction
+#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN  "SUBTRACTION.RES.FSIGMA.MEAN"      // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV "SUBTRACTION.RES.FSIGMA.STDEV"      // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN    "SUBTRACTION.RES.FMIN.MEAN"      // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV   "SUBTRACTION.RES.FMIN.STDEV"      // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN    "SUBTRACTION.RES.FMAX.MEAN"      // RMS stamp deviation
+#define PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV   "SUBTRACTION.RES.FMAX.STDEV"      // RMS stamp deviation
 // Derive QA information about the subtraction
 bool pmSubtractionAnalysis(

trunk/psModules/src/imcombine/pmSubtractionEquation.c

-              r26035
+              r26893
 //#define TESTING                         // TESTING output for debugging; may not work with threads!
+#define USE_WEIGHT                      // Include weight (1/variance) in equation?
+//#define USE_WEIGHT                      // Include weight (1/variance) in equation?
+//#define USE_WINDOW                      // Include weight (1/variance) in equation?
 …
 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
                                   psVector *vector, // Least-squares vector, updated
+                                  double *norm,     // Normalisation, updated
                                   const psKernel *input, // Input image (target)
                                   const psKernel *reference, // Reference image (convolution source)
                                   const psKernel *weight,  // Weight image
+                                  const psKernel *window,  // Window image
                                   const psArray *convolutions,         // Convolutions for each kernel
                                   const pmSubtractionKernels *kernels, // Kernels
                                   const psImage *polyValues, // Spatial polynomial values
+                                  int footprint // (Half-)Size of stamp
+                                  int footprint, // (Half-)Size of stamp
+                                  int normWindow, // Window (half-)size for normalisation measurement
+                                  const pmSubtractionEquationCalculationMode mode
+                                  )
+{
 …
     // Evaluate polynomial-polynomial terms
+    // XXX we can skip this if we are not calculating kernel coeffs
     for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) {
         for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) {
 …
+    }
+    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
+    // choose to calculate
+    psImageInit(matrix, 0.0);
+    psVectorInit(vector, 1.0);
+    for (int i = 0; i < matrix->numCols; i++) {
+        matrix->data.F64[i][i] = 1.0;
+    }
+    // the order of the elements in the matrix and vector is:
+    // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0]
+    // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0]
+    // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m]
+    // normalization
+    // bg 0, bg 1, bg 2 (only 0 is currently used?)
     for (int i = 0; i < numKernels; i++) {
 …
                 for (int x = - footprint; x <= footprint; x++) {
                     double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
+#ifdef USE_WEIGHT
+                    cc *= weight->kernel[y][x];
+#endif
+                    if (weight) {
+                        cc *= weight->kernel[y][x];
+                    }
+                    if (window) {
+                        cc *= window->kernel[y][x];
+                    }
                     sumCC += cc;
+                }
+            }
+            // Spatial variation
+            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                    double value = sumCC * poly2[iTerm][jTerm];
+                    matrix->data.F64[iIndex][jIndex] = value;
+                    matrix->data.F64[jIndex][iIndex] = value;
+            // Spatial variation of kernel coeffs
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                        double value = sumCC * poly2[iTerm][jTerm];
+                        matrix->data.F64[iIndex][jIndex] = value;
+                        matrix->data.F64[jIndex][iIndex] = value;
+                    }
+                }
+            }
 …
                 double rc = ref * conv;
                 double c = conv;
+#ifdef USE_WEIGHT
+                float wtVal = weight->kernel[y][x];
+                ic *= wtVal;
+                rc *= wtVal;
+                c *= wtVal;
+#endif
+                if (weight) {
+                    float wtVal = weight->kernel[y][x];
+                    ic *= wtVal;
+                    rc *= wtVal;
+                    c *= wtVal;
+                }
+                if (window) {
+                    float winVal = window->kernel[y][x];
+                    ic *= winVal;
+                    rc *= winVal;
+                    c  *= winVal;
+                }
                 sumIC += ic;
                 sumRC += rc;
 …
             double normTerm = sumRC * poly[iTerm];
             double bgTerm = sumC * poly[iTerm];
+            matrix->data.F64[iIndex][normIndex] = normTerm;
+            matrix->data.F64[normIndex][iIndex] = normTerm;
+            matrix->data.F64[iIndex][bgIndex] = bgTerm;
+            matrix->data.F64[bgIndex][iIndex] = bgTerm;
+            vector->data.F64[iIndex] = sumIC * poly[iTerm];
+            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
+                matrix->data.F64[iIndex][normIndex] = normTerm;
+                matrix->data.F64[normIndex][iIndex] = normTerm;
+            }
+            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
+                matrix->data.F64[iIndex][bgIndex] = bgTerm;
+                matrix->data.F64[bgIndex][iIndex] = bgTerm;
+            }
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+                vector->data.F64[iIndex] = sumIC * poly[iTerm];
+                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
+                    // subtract norm * sumRC * poly[iTerm]
+                    psAssert (kernels->solution1, "programming error: define solution first!");
+                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
+                    vector->data.F64[iIndex] -= norm * normTerm;
+                }
+            }
+        }
+    }
 …
     double sumR = 0.0;                  // Sum of the reference
     double sumI = 0.0;                  // Sum of the input
+    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
     for (int y = - footprint; y <= footprint; y++) {
         for (int x = - footprint; x <= footprint; x++) {
 …
             double rr = PS_SQR(ref);
             double one = 1.0;
+#ifdef USE_WEIGHT
+            float wtVal = weight->kernel[y][x];
+            rr *= wtVal;
+            ir *= wtVal;
+            in *= wtVal;
+            ref *= wtVal;
+            one *= wtVal;
+#endif
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
+                normI1 += ref;
+                normI2 += in;
+            }
+            if (weight) {
+                float wtVal = weight->kernel[y][x];
+                rr *= wtVal;
+                ir *= wtVal;
+                in *= wtVal;
+                ref *= wtVal;
+                one *= wtVal;
+            }
+            if (window) {
+                float  winVal = window->kernel[y][x];
+                rr      *= winVal;
+                ir      *= winVal;
+                in      *= winVal;
+                ref *= winVal;
+                one *= winVal;
+            }
             sumRR += rr;
             sumIR += ir;
 …
+        }
+    }
+    matrix->data.F64[normIndex][normIndex] = sumRR;
+    matrix->data.F64[bgIndex][bgIndex] = sum1;
+    matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR;
+    vector->data.F64[normIndex] = sumIR;
+    vector->data.F64[bgIndex] = sumI;
+    *norm = normI2 / normI1;
+    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
+        matrix->data.F64[normIndex][normIndex] = sumRR;
+        vector->data.F64[normIndex] = sumIR;
+        // subtract sum over kernels * kernel solution
+    }
+    if (mode & PM_SUBTRACTION_EQUATION_BG) {
+        matrix->data.F64[bgIndex][bgIndex] = sum1;
+        vector->data.F64[bgIndex] = sumI;
+    }
+    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
+        matrix->data.F64[normIndex][bgIndex] = sumR;
+        matrix->data.F64[bgIndex][normIndex] = sumR;
+    }
+    // check for any NAN values in the result, skip if found:
+    for (int iy = 0; iy < matrix->numRows; iy++) {
+        for (int ix = 0; ix < matrix->numCols; ix++) {
+            if (!isfinite(matrix->data.F64[iy][ix])) {
+                fprintf (stderr, "WARNING: NAN in matrix\n");
+                return false;
+            }
+        }
+    }
+    for (int ix = 0; ix < vector->n; ix++) {
+        if (!isfinite(vector->data.F64[ix])) {
+            fprintf (stderr, "WARNING: NAN in vector\n");
+            return false;
+        }
+    }
     return true;
+}
 // Calculate the least-squares matrix and vector for dual convolution
+static bool calculateDualMatrixVector(psImage *matrix1, // Least-squares matrix, updated
+                                      psVector *vector1, // Least-squares vector, updated
+                                      psImage *matrix2,  // Least-squares matrix, updated
+                                      psVector *vector2, // Least-squares vector, updated
+                                      psImage *matrixX,  // Cross-matrix
+static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated
+                                      psVector *vector, // Least-squares vector, updated
+                                      double *norm,     // Normalisation, updated
                                       const psKernel *image1, // Image 1
                                       const psKernel *image2, // Image 2
                                       const psKernel *weight,  // Weight image
+                                      const psKernel *window,  // Window image
                                       const psArray *convolutions1, // Convolutions of image 1 for each kernel
                                       const psArray *convolutions2, // Convolutions of image 2 for each kernel
                                       const pmSubtractionKernels *kernels, // Kernels
                                       const psImage *polyValues, // Spatial polynomial values
+                                      int footprint // (Half-)Size of stamp
+                                      int footprint, // (Half-)Size of stamp
+                                      int normWindow, // Window (half-)size for normalisation measurement
+                                      const pmSubtractionEquationCalculationMode mode
+                                      )
+{
-    // A_ij = A_i A_j
-    // B_ij = B_i B_j
-    // C_ij = A_i B_j
-    // d_i = A_i I_2
-    // e_i = B_i I_2
-    // A_i = I_1 * k_i
-    // B_i = I_2 * k_i
-    // Background: A_i = 1.0
-    // Normalisation: A_i = I_1
     int numKernels = kernels->num;                      // Number of kernels
     int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
 …
     double poly[numPoly];                                 // Polynomial terms
     double poly2[numPoly][numPoly];                       // Polynomial-polynomial values
+    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
+    int numParams = numKernels * numPoly + 1 + numBackground;       // Number of regular parameters
+    int numParams2 = numKernels * numPoly;                          // Number of additional parameters for dual
+    int numDual = numParams + numParams2;                           // Total number of parameters for dual
+    psAssert(matrix &&
+             matrix->type.type == PS_TYPE_F64 &&
+             matrix->numCols == numDual &&
+             matrix->numRows == numDual,
+             "Least-squares matrix is bad.");
+    psAssert(vector &&
+             vector->type.type == PS_TYPE_F64 &&
+             vector->n == numDual,
+             "Least-squares vector is bad.");
     // Evaluate polynomial-polynomial terms
 …
+    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
+    // choose to calculate
+    psImageInit(matrix, 0.0);
+    psVectorInit(vector, 1.0);
+    for (int i = 0; i < matrix->numCols; i++) {
+        matrix->data.F64[i][i] = 1.0;
+    }
     for (int i = 0; i < numKernels; i++) {
         psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i
 …
             psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j
             double sumAA = 0.0;         // Sum of convolution products for matrix A
             double sumBB = 0.0;         // Sum of convolution products for matrix B
             double sumAB = 0.0;         // Sum of convolution products for matrix C
+            double sumAA = 0.0;         // Sum of convolution products between image 1
+            double sumBB = 0.0;         // Sum of convolution products between image 2
+            double sumAB = 0.0;         // Sum of convolution products across images 1 and 2
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
 …
                     double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
+#ifdef USE_WEIGHT
+                    float wtVal = weight->kernel[y][x];
+                    aa *= wtVal;
+                    bb *= wtVal;
+                    ab *= wtVal;
+#endif
+                    if (weight) {
+                        float wtVal = weight->kernel[y][x];
+                        aa *= wtVal;
+                        bb *= wtVal;
+                        ab *= wtVal;
+                    }
+                    if (window) {
+                        float wtVal = window->kernel[y][x];
+                        aa *= wtVal;
+                        bb *= wtVal;
+                        ab *= wtVal;
+                    }
                     sumAA += aa;
                     sumBB += bb;
 …
+            }
+            // Spatial variation
+            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                    double aa = sumAA * poly2[iTerm][jTerm];
+                    double bb = sumBB * poly2[iTerm][jTerm];
+                    double ab = sumAB * poly2[iTerm][jTerm];
+                    matrix1->data.F64[iIndex][jIndex] = aa;
+                    matrix1->data.F64[jIndex][iIndex] = aa;
+                    matrix2->data.F64[iIndex][jIndex] = bb;
+                    matrix2->data.F64[jIndex][iIndex] = bb;
+                    matrixX->data.F64[iIndex][jIndex] = ab;
+            // Spatial variation of kernel coeffs
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                        double aa = sumAA * poly2[iTerm][jTerm];
+                        double bb = sumBB * poly2[iTerm][jTerm];
+                        double ab = sumAB * poly2[iTerm][jTerm];
+                        matrix->data.F64[iIndex][jIndex] = aa;
+                        matrix->data.F64[jIndex][iIndex] = aa;
+                        matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
+                        matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
+                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
+                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
+                    }
+                }
+            }
 …
                 for (int x = - footprint; x <= footprint; x++) {
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
+#ifdef USE_WEIGHT
+                    ab *= weight->kernel[y][x];
+#endif
+                    if (weight) {
+                        ab *= weight->kernel[y][x];
+                    }
+                    if (window) {
+                        ab *= window->kernel[y][x];
+                    }
                     sumAB += ab;
+                }
+            }
+            // Spatial variation
+            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                    double ab = sumAB * poly2[iTerm][jTerm];
+                    matrixX->data.F64[iIndex][jIndex] = ab;
+                }
+            }
+        }
+        double sumAI2 = 0.0;            // Sum of A.I_2 products (for vector 1)
+        double sumBI2 = 0.0;            // Sum of B.I_2 products (for vector 2)
+        double sumAI1 = 0.0;            // Sum of A.I_1 products (for matrix 1, normalisation)
+        double sumA = 0.0;              // Sum of A (for matrix 1, background)
+        double sumBI1 = 0.0;            // Sum of B.I_1 products (for matrix X, normalisation)
+        double sumB = 0.0;              // Sum of B products (for matrix X, background)
+        double sumI2 = 0.0;             // Sum of I_2 (for vector 1, background)
+        double sumI1I2 = 0.0;           // Sum of I_1.I_2 (for vector 1, normalisation)
+            // Spatial variation of kernel coeffs
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+                        double ab = sumAB * poly2[iTerm][jTerm];
+                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
+                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
+                    }
+                }
+            }
+        }
+        double sumAI2 = 0.0;            // Sum of A.I_2 products (for vector)
+        double sumBI2 = 0.0;            // Sum of B.I_2 products (for vector)
+        double sumAI1 = 0.0;            // Sum of A.I_1 products (for matrix, normalisation)
+        double sumA = 0.0;              // Sum of A (for matrix, background)
+        double sumBI1 = 0.0;            // Sum of B.I_1 products (for matrix, normalisation)
+        double sumB = 0.0;              // Sum of B products (for matrix, background)
+        double sumI2 = 0.0;             // Sum of I_2 (for vector, background)
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
                 float a = iConv1->kernel[y][x];
                 float b = iConv2->kernel[y][x];
+                double a = iConv1->kernel[y][x];
+                double b = iConv2->kernel[y][x];
                 float i1 = image1->kernel[y][x];
                 float i2 = image2->kernel[y][x];
 …
                 double ai1 = a * i1;
                 double bi1 = b * i1;
+                double i1i2 = i1 * i2;
+#ifdef USE_WEIGHT
+                float wtVal = weight->kernel[y][x];
+                ai2 *= wtVal;
+                bi2 *= wtVal;
+                ai1 *= wtVal;
+                bi1 *= wtVal;
+                i1i2 *= wtVal;
+                a *= wtVal;
+                b *= wtVal;
+                i2 *= wtVal;
+#endif
+                if (weight) {
+                    float wtVal = weight->kernel[y][x];
+                    ai2 *= wtVal;
+                    bi2 *= wtVal;
+                    ai1 *= wtVal;
+                    bi1 *= wtVal;
+                    a *= wtVal;
+                    b *= wtVal;
+                    i2 *= wtVal;
+                }
+                if (window) {
+                    float wtVal = window->kernel[y][x];
+                    ai2 *= wtVal;
+                    bi2 *= wtVal;
+                    ai1 *= wtVal;
+                    bi1 *= wtVal;
+                    a *= wtVal;
+                    b *= wtVal;
+                    i2 *= wtVal;
+                }
                 sumAI2 += ai2;
                 sumBI2 += bi2;
 …
                 sumB += b;
                 sumI2 += i2;
-                sumI1I2 += i1i2;
+            }
+        }
 …
             double bi2 = sumBI2 * poly[iTerm];
             double ai1 = sumAI1 * poly[iTerm];
             double a = sumA * poly[iTerm];
+            double a   = sumA * poly[iTerm];
             double bi1 = sumBI1 * poly[iTerm];
+            double b = sumB * poly[iTerm];
+            matrix1->data.F64[iIndex][normIndex] = ai1;
+            matrix1->data.F64[normIndex][iIndex] = ai1;
+            matrix1->data.F64[iIndex][bgIndex] = a;
+            matrix1->data.F64[bgIndex][iIndex] = a;
+            vector1->data.F64[iIndex] = ai2;
+            vector2->data.F64[iIndex] = bi2;
+            matrixX->data.F64[iIndex][normIndex] = bi1;
+            matrixX->data.F64[iIndex][bgIndex] = b;
+        }
+    }
+    double sumI1 = 0.0;                 // Sum of I_1 (for matrix 1, background-normalisation)
+    double sumI1I1 = 0.0;               // Sum of I_1^2 (for matrix 1, normalisation-normalisation)
+    double sum1 = 0.0;                  // Sum of 1 (for matrix 1, background-background)
+    double sumI2 = 0.0;                 // Sum of I_2 (for vector 1, background)
+    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector 1, normalisation)
+            double b   = sumB * poly[iTerm];
+            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
+                matrix->data.F64[iIndex][normIndex] = ai1;
+                matrix->data.F64[normIndex][iIndex] = ai1;
+                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
+                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
+            }
+            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
+                matrix->data.F64[iIndex][bgIndex] = a;
+                matrix->data.F64[bgIndex][iIndex] = a;
+                matrix->data.F64[iIndex + numParams][bgIndex] = b;
+                matrix->data.F64[bgIndex][iIndex + numParams] = b;
+            }
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+                vector->data.F64[iIndex] = ai2;
+                vector->data.F64[iIndex + numParams] = bi2;
+                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
+                    // subtract norm * sumRC * poly[iTerm]
+                    psAssert (kernels->solution1, "programming error: define solution first!");
+                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
+                    vector->data.F64[iIndex] -= norm * ai1;
+                    vector->data.F64[iIndex + numParams] -= norm * bi1;
+                }
+            }
+        }
+    }
+    double sumI1 = 0.0;                 // Sum of I_1 (for matrix, background-normalisation)
+    double sumI1I1 = 0.0;               // Sum of I_1^2 (for matrix, normalisation-normalisation)
+    double sum1 = 0.0;                  // Sum of 1 (for matrix, background-background)
+    double sumI2 = 0.0;                 // Sum of I_2 (for vector, background)
+    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
+    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
     for (int y = - footprint; y <= footprint; y++) {
         for (int x = - footprint; x <= footprint; x++) {
             float i1 = image1->kernel[y][x];
             float i2 = image2->kernel[y][x];
+            double i1 = image1->kernel[y][x];
+            double i2 = image2->kernel[y][x];
             double i1i1 = i1 * i1;
 …
             double i1i2 = i1 * i2;
+#ifdef USE_WEIGHT
+            float wtVal = weight->kernel[y][x];
+            i1 *= wtVal;
+            i1i1 *= wtVal;
+            one *= wtVal;
+            i2 *= wtVal;
+            i1i2 *= wtVal;
+#endif
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
+                normI1 += i1;
+                normI2 += i2;
+            }
+            if (weight) {
+                float wtVal = weight->kernel[y][x];
+                i1 *= wtVal;
+                i1i1 *= wtVal;
+                one *= wtVal;
+                i2 *= wtVal;
+                i1i2 *= wtVal;
+            }
+            if (window) {
+                float wtVal = window->kernel[y][x];
+                i1 *= wtVal;
+                i1i1 *= wtVal;
+                one *= wtVal;
+                i2 *= wtVal;
+                i1i2 *= wtVal;
+            }
             sumI1 += i1;
             sumI1I1 += i1i1;
 …
+        }
+    }
+    matrix1->data.F64[bgIndex][normIndex] = sumI1;
+    matrix1->data.F64[normIndex][bgIndex] = sumI1;
+    matrix1->data.F64[normIndex][normIndex] = sumI1I1;
+    matrix1->data.F64[bgIndex][bgIndex] = sum1;
+    vector1->data.F64[bgIndex] = sumI2;
+    vector1->data.F64[normIndex] = sumI1I2;
+    *norm = normI2 / normI1;
+    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
+        matrix->data.F64[normIndex][normIndex] = sumI1I1;
+        vector->data.F64[normIndex] = sumI1I2;
+    }
+    if (mode & PM_SUBTRACTION_EQUATION_BG) {
+        matrix->data.F64[bgIndex][bgIndex] = sum1;
+        vector->data.F64[bgIndex] = sumI2;
+    }
+    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
+        matrix->data.F64[bgIndex][normIndex] = sumI1;
+        matrix->data.F64[normIndex][bgIndex] = sumI1;
+    }
+    // check for any NAN values in the result, skip if found:
+    for (int iy = 0; iy < matrix->numRows; iy++) {
+        for (int ix = 0; ix < matrix->numCols; ix++) {
+            if (!isfinite(matrix->data.F64[iy][ix])) {
+                fprintf (stderr, "WARNING: NAN in matrix\n");
+                return false;
+            }
+        }
+    }
+    for (int ix = 0; ix < vector->n; ix++) {
+        if (!isfinite(vector->data.F64[ix])) {
+            fprintf (stderr, "WARNING: NAN in vector\n");
+            return false;
+        }
+    }
     return true;
+}
+// Merge dual matrices and vectors into single matrix equation
+// Have: Aa = Ct.b + d
+// Have: Ca = Bb + e
+// Set: F = ( A -Ct ;  C -B )
+// Set: g = ( a ; b )
+// Set: h = ( d ; e )
+// So that we combine the above two equations: Fg = h
+static bool calculateEquationDual(psImage **outMatrix,
+                                  psVector **outVector,
+                                  const psImage *sumMatrix1,
+                                  const psImage *sumMatrix2,
+                                  const psImage *sumMatrixX,
+                                  const psVector *sumVector1,
+                                  const psVector *sumVector2
+                                  )
+{
+    psAssert(sumMatrix1 && sumMatrix2 && sumMatrixX, "Require input matrices");
+    psAssert(sumVector1 && sumVector2, "Require input vectors");
+    int num1 = sumVector1->n;   // Number of parameters in first set
+    int num2 = sumVector2->n;   // Number of parameters in second set
+    int num = num1 + num2;      // Number of parameters in new set
+    psAssert(sumMatrix1->type.type == PS_TYPE_F64 &&
+             sumMatrix2->type.type == PS_TYPE_F64 &&
+             sumMatrixX->type.type == PS_TYPE_F64 &&
+             sumVector1->type.type == PS_TYPE_F64 &&
+             sumVector2->type.type == PS_TYPE_F64,
+             "Require input type is F64");
+    psAssert(outMatrix, "Require output matrix");
+    psAssert(outVector, "Require output vector");
+    if (!*outMatrix) {
+        *outMatrix = psImageAlloc(num, num, PS_TYPE_F64);
+    }
+    if (!*outVector) {
+        *outVector = psVectorAlloc(num, PS_TYPE_F64);
+    }
+    psImage *matrix = *outMatrix;
+    psVector *vector = *outVector;
+    psAssert(sumMatrix1->numCols == num1 && sumMatrix1->numRows == num1, "Require size NxN");
+    psAssert(sumMatrix2->numCols == num2 && sumMatrix2->numRows == num2, "Require size MxM");
+    psAssert(sumMatrixX->numCols == num1 && sumMatrixX->numRows == num2, "Require size MxN");
+    memcpy(vector->data.F64, sumVector1->data.F64, num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+    memcpy(&vector->data.F64[num1], sumVector2->data.F64, num2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+    for (int i = 0; i < num1; i++) {
+        memcpy(matrix->data.F64[i], sumMatrix1->data.F64[i], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        for (int j = 0, k = num1; j < num2; j++, k++) {
+            matrix->data.F64[i][k] = - sumMatrixX->data.F64[j][i];
+        }
+    }
+    for (int i1 = 0, i2 = num1; i1 < num2; i1++, i2++) {
+        memcpy(matrix->data.F64[i2], sumMatrixX->data.F64[i1], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        for (int j = 0, k = num1; j < num2; j++, k++) {
+            matrix->data.F64[i2][k] = - sumMatrix2->data.F64[i1][j];
+        }
+    }
+    return true;
+}
+#if 1
 // Add in penalty term to least-squares vector
 static bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
+bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
                              psVector *vector,                    // Vector to which to add in penalty term
                              const pmSubtractionKernels *kernels, // Kernel parameters
 …
     int spatialOrder = kernels->spatialOrder; // Order of spatial variations
     int numKernels = kernels->num; // Number of kernel components
+    int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations
+    int numParams = numKernels * numSpatial;                 // Number of kernel parameters
+    // order is :
+    // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0]
+    // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0]
+    // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1]
+    // [norm]
+    // [bg]
+    // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0]
+    // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0]
+    // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1]
     for (int i = 0; i < numKernels; i++) {
         for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
                 // Contribution to chi^2: a_i^2 P_i
+                matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
+                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
+                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
+                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
+                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
+                    // penalties scale with second moments
+                    //
+                }
+            }
+        }
 …
     return true;
+}
+# endif
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 …
 // Calculate the value of a polynomial, specified by coefficients and polynomial values
 double p_pmSubtractionCalculatePolynomial(const psVector *coeff, // Coefficients
                                                  const psImage *polyValues, // Polynomial values
                                                  int order, // Order of polynomials
                                                  int index, // Index at which to begin
                                                  int step // Step between subsequent indices
+                                                 )
+                                          const psImage *polyValues, // Polynomial values
+                                          int order, // Order of polynomials
+                                          int index, // Index at which to begin
+                                          int step // Step between subsequent indices
+                                          )
+{
     double sum = 0.0;                   // Value of the polynomial sum
 …
 double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, const psImage *polyValues,
                                            int index, bool wantDual)
+                                    int index, bool wantDual)
+{
 #if 0
 …
     const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
+    return pmSubtractionCalculateEquationStamp(stamps, kernels, index);
+    pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
+    return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode);
+}
 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
                                          int index)
+                                         int index, const pmSubtractionEquationCalculationMode mode)
+{
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
 …
     int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
+    // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial).
+    // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3
     // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the
     // number of coefficients for the spatial polynomial, normalisation and a constant background offset.
     int numParams = numKernels * numSpatial + 1 + numBackground;
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        // An additional image is convolved
+        numParams += numKernels * numSpatial;
+    }
     pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest
     psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state.");
     // Generate convolutions
+    // Generate convolutions: these are generated once and saved
     if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", index);
 …
 #endif
+    // XXX visualize the set of convolved stamps
     psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
                                                     stamp->xNorm, stamp->yNorm); // Polynomial terms
     bool new = stamp->vector1 ? false : true; // Is this a new run?
+    bool new = stamp->vector ? false : true; // Is this a new run?
     if (new) {
         stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
         stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64);
+        stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+        stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64);
+    }
 #ifdef TESTING
     psImageInit(stamp->matrix1, NAN);
     psVectorInit(stamp->vector1, NAN);
+    psImageInit(stamp->matrix, NAN);
+    psVectorInit(stamp->vector, NAN);
 #endif
     bool status;                    // Status of least-squares matrix/vector calculation
+    psKernel *weight = NULL;
+    psKernel *window = NULL;
+#ifdef USE_WEIGHT
+    weight = stamp->weight;
+#endif
+#ifdef USE_WINDOW
+    window = stamps->window;
+#endif
     switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_1:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
                                        stamp->weight, stamp->convolutions1, kernels, polyValues,
                                        footprint);
+        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
+                                       weight, window, stamp->convolutions1, kernels,
+                                       polyValues, footprint, stamps->normWindow, mode);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
                                        stamp->weight, stamp->convolutions2, kernels, polyValues,
                                        footprint);
+        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
+                                       weight, window, stamp->convolutions2, kernels,
+                                       polyValues, footprint, stamps->normWindow, mode);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
+        if (new) {
+            stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
+            stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64);
+            stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64);
+        }
+#ifdef TESTING
+        psImageInit(stamp->matrix2, NAN);
+        psImageInit(stamp->matrixX, NAN);
+        psVectorInit(stamp->vector2, NAN);
+#endif
+        status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2,
+                                           stamp->matrixX, stamp->image1, stamp->image2, stamp->weight,
+                                           stamp->convolutions1, stamp->convolutions2, kernels, polyValues,
+                                           footprint);
+        status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm,
+                                           stamp->image1, stamp->image2,
+                                           weight, window, stamp->convolutions1, stamp->convolutions2,
+                                           kernels, polyValues, footprint, stamps->normWindow, mode);
         break;
       default:
 …
         stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
         psWarning("Rejecting stamp %d (%d,%d) because of bad equation",
                   index, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+                  index, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
     } else {
         stamp->status = PM_SUBTRACTION_STAMP_USED;
 …
+    {
         psString matrixName = NULL;
         psStringAppend(&matrixName, "matrix1_%d.fits", index);
+        psStringAppend(&matrixName, "matrix_%d.fits", index);
         psFits *matrixFile = psFitsOpen(matrixName, "w");
         psFree(matrixName);
         psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL);
+        psFitsWriteImage(matrixFile, NULL, stamp->matrix, 0, NULL);
         psFitsClose(matrixFile);
         matrixName = NULL;
         psStringAppend(&matrixName, "vector1_%d.fits", index);
         psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64);
         memcpy(dummy->data.F64[0], stamp->vector1->data.F64,
                PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n);
+        psStringAppend(&matrixName, "vector_%d.fits", index);
+        psImage *dummy = psImageAlloc(stamp->vector->n, 1, PS_TYPE_F64);
+        memcpy(dummy->data.F64[0], stamp->vector->data.F64,
+               PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector->n);
         matrixFile = psFitsOpen(matrixName, "w");
         psFree(matrixName);
 …
         psFree(dummy);
         psFitsClose(matrixFile);
-        if (stamp->vector2) {
-            matrixName = NULL;
-            psStringAppend(&matrixName, "vector2_%d.fits", index);
-            dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);
-            memcpy(dummy->data.F64[0], stamp->vector2->data.F64,
-                   PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);
-            matrixFile = psFitsOpen(matrixName, "w");
-            psFree(matrixName);
-            psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
-            psFree(dummy);
-            psFitsClose(matrixFile);
+        }
-        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            matrixName = NULL;
-            psStringAppend(&matrixName, "matrix2_%d.fits", index);
-            matrixFile = psFitsOpen(matrixName, "w");
-            psFree(matrixName);
-            psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL);
-            psFitsClose(matrixFile);
-            matrixName = NULL;
-            psStringAppend(&matrixName, "matrixX_%d.fits", index);
-            matrixFile = psFitsOpen(matrixName, "w");
-            psFree(matrixName);
-            psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL);
-            psFitsClose(matrixFile);
+        }
+    }
 #endif
 …
+}
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels)
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
+                                    const pmSubtractionEquationCalculationMode mode)
+{
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
 …
+        }
+        if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
+            psAbort ("bad stamp");
+        }
+        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
+            psAbort ("bad stamp");
+        }
         if (pmSubtractionThreaded()) {
             psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION");
 …
             psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32);
             if (!psThreadJobAddPending(job)) {
                 psFree(job);
 …
             psFree(job);
         } else {
             pmSubtractionCalculateEquationStamp(stamps, kernels, i);
+            pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode);
+        }
+    }
 …
     pmSubtractionVisualPlotLeastSquares(stamps);
+    pmSubtractionVisualShowKernels((pmSubtractionKernels  *)kernels);
+    pmSubtractionVisualShowBasis(stamps);
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
 …
+}
+bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps)
+// private functions used on pmSubtractionSolveEquation
+bool psVectorWriteFile (char *filename, const psVector *vector);
+bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header);
+psImage *p_pmSubSolve_wUt (psVector *w, psImage *U);
+psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt);
+bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask);
+bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B);
+bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB);
+bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB);
+bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x);
+bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y);
+bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
+psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w);
+double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
+bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
+                                const pmSubtractionStampList *stamps,
+                                const pmSubtractionEquationCalculationMode mode)
+{
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
 …
     // Check inputs
+    int numParams = -1;                // Number of parameters
+    int numParams2 = 0;                // Number of parameters for part solution (DUAL mode)
+    int numKernels = kernels->num;      // Number of kernel basis functions
+    int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
+    int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
+    int numParams = numKernels * numSpatial + 1 + numBackground;    // Number of parameters being solved for
+    int numSolution1 = numParams, numSolution2 = 0;                 // Number of parameters for each solution
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        // An additional image is convolved
+        numSolution2 = numKernels * numSpatial;
+        numParams += numSolution2;
+    }
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
 …
+        }
+        PS_ASSERT_VECTOR_NON_NULL(stamp->vector1, false);
+        if (numParams == -1) {
+            numParams = stamp->vector1->n;
+        }
+        PS_ASSERT_VECTOR_SIZE(stamp->vector1, (long)numParams, false);
+        PS_ASSERT_VECTOR_TYPE(stamp->vector1, PS_TYPE_F64, false);
+        PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false);
+        PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false);
+        PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false);
+        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+            PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false);
+            PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false);
+            if (numParams2 == 0) {
+                numParams2 = stamp->matrix2->numCols;
+            }
+            PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false);
+            PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, 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->vector2, false);
+            PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false);
+            PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false);
+        }
+    }
+    if (numParams == -1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found.");
+        return NULL;
+        PS_ASSERT_VECTOR_NON_NULL(stamp->vector, 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->matrix, false);
+        PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, false);
+        PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, false);
+    }
 …
         psVectorInit(sumVector, 0.0);
         psImageInit(sumMatrix, 0.0);
+        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
         int numStamps = 0;              // Number of good stamps
         for (int i = 0; i < stamps->num; i++) {
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
+#ifdef TESTING
+              // XXX double-check for NAN in data:
+                for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
+                    for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
+                        if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
+                            fprintf (stderr, "WARNING: NAN in matrix1\n");
+                        }
+                    }
+                }
+                for (int ix = 0; ix < stamp->vector1->n; ix++) {
+                    if (!isfinite(stamp->vector1->data.F64[ix])) {
+                        fprintf (stderr, "WARNING: NAN in vector1\n");
+                    }
+                }
+#endif
+                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
+                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
+                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
+                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
+                psVectorAppend(norms, stamp->norm);
                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
                 numStamps++;
 …
+        }
-#ifdef TESTING
-        for (int ix = 0; ix < sumVector->n; ix++) {
-            if (!isfinite(sumVector->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
-#endif
 #if 0
         int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
 …
 #endif
+#ifdef TESTING
+        for (int ix = 0; ix < sumVector->n; ix++) {
+            if (!isfinite(sumVector->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
+        psVector *solution = NULL;                       // Solution to equation!
+        solution = psVectorAlloc(numParams, PS_TYPE_F64);
+        psVectorInit(solution, 0);
+#if 0
+        // Regular, straight-forward solution
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+#else
+        {
+            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
+            psFits *fits = psFitsOpen("matrixInv.fits", "w");
+            psFitsWriteImage(fits, NULL, inverse, 0, NULL);
+            psFitsClose(fits);
+            psFree(inverse);
+        }
+        {
+            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
+            psImage *Xt = psMatrixTranspose(NULL, X);
+            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
+            psFits *fits = psFitsOpen("matrixErr.fits", "w");
+            psFitsWriteImage(fits, NULL, XtX, 0, NULL);
+            psFitsClose(fits);
+            psFree(X);
+            psFree(Xt);
+            psFree(XtX);
+        }
+#endif
+        psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
+        psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
+            // Solve normalisation and background separately
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
+            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
+            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
+                psFree(stats);
+                psFree(sumMatrix);
+                psFree(sumVector);
+                psFree(norms);
+                return false;
+            }
+            double normValue = stats->robustMedian;
+            // double bgValue = 0.0;
+            psFree(stats);
+            fprintf(stderr, "Norm: %lf\n", normValue);
+            // Solve kernel components
+            for (int i = 0; i < numSolution1; i++) {
+                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
+                sumMatrix->data.F64[i][normIndex] = 0.0;
+                sumMatrix->data.F64[normIndex][i] = 0.0;
+            }
+            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
+            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
+            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
+            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
+            sumVector->data.F64[normIndex] = 0.0;
+            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+            solution->data.F64[normIndex] = normValue;
+        }
+# endif
+        if (!kernels->solution1) {
+            kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+            psVectorInit(kernels->solution1, 0.0);
+        }
+        // only update the solutions that we chose to calculate:
+        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_BG) {
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+            int numKernels = kernels->num;
+            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
+            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+            for (int i = 0; i < numKernels * numPoly; i++) {
+                kernels->solution1->data.F64[i] = solution->data.F64[i];
+            }
+        }
+        psFree(solution);
+        psFree(sumVector);
         psFree(sumMatrix);
-        if (!luMatrix) {
-            psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
-            psFree(sumVector);
-            psFree(luMatrix);
-            psFree(permutation);
-            return NULL;
+        }
-        kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation);
 #ifdef TESTING
 …
         for (int ix = 0; ix < kernels->solution1->n; ix++) {
             if (!isfinite(kernels->solution1->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
+#endif
+        psFree(sumVector);
+        psFree(luMatrix);
+        psFree(permutation);
+        if (!kernels->solution1) {
+            psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
+            return NULL;
+        }
+                fprintf (stderr, "WARNING: NAN in vector\n");
+            }
+        }
+#endif
     } else {
         // Dual convolution solution
         // Accumulation of stamp matrices/vectors
+        psImage *sumMatrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+        psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64);
+        psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64);
+        psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64);
+        psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64);
+        psImageInit(sumMatrix1, 0.0);
+        psImageInit(sumMatrix2, 0.0);
+        psImageInit(sumMatrixX, 0.0);
+        psVectorInit(sumVector1, 0.0);
+        psVectorInit(sumVector2, 0.0);
+        psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
+        psImageInit(sumMatrix, 0.0);
+        psVectorInit(sumVector, 0.0);
+        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
         int numStamps = 0;              // Number of good stamps
 …
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
                 (void)psBinaryOp(sumMatrix1, sumMatrix1, "+", stamp->matrix1);
                 (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2);
+                (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX);
                 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);
+                (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2);
+                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
+                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
+                psVectorAppend(norms, stamp->norm);
                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
                 numStamps++;
 …
+        }
-        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
-        calculatePenalty(sumMatrix1, sumVector1, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]);
-        calculatePenalty(sumMatrix2, sumVector2, kernels, -sumMatrix1->data.F64[bgIndex][bgIndex]);
-        psImage *sumMatrix = NULL;      // Combined matrix
-        psVector *sumVector = NULL;     // Combined vector
-        calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
-                              sumMatrixX, sumVector1, sumVector2);
 #ifdef TESTING
+        psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
+        psVectorWriteFile("sumVector.dat", sumVector);
+#endif
+#if 1
+        // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
+        // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0);
+#endif
+        psVector *solution = NULL;                       // Solution to equation!
+        solution = psVectorAlloc(numParams, PS_TYPE_F64);
+        psVectorInit(solution, 0);
+#if 0
+        // Regular, straight-forward solution
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+#else
+        {
+            psFits *fits = psFitsOpen("sumMatrix.fits", "w");
+            psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
+            psFitsClose(fits);
+        }
+        {
+            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
+            psFits *fits = psFitsOpen("sumVector.fits", "w");
+            for (int i = 0; i < numParams + numParams2; i++) {
+                image->data.F64[0][i] = sumVector->data.F64[i];
+            }
+            psFitsWriteImage(fits, NULL, image, 0, NULL);
+            psFree(image);
+            psFitsClose(fits);
+        }
+#endif
+        psVector *solution = NULL;                       // Solution to equation!
+        {
+            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
+            if (!solution) {
+                psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
+            // Solve normalisation and background separately
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
+#if 0
+            psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64);
+            psVector *normVector = psVectorAlloc(2, PS_TYPE_F64);
+            normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex];
+            normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex];
+            normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex];
+            normVector->data.F64[0] = sumVector->data.F64[normIndex];
+            normVector->data.F64[1] = sumVector->data.F64[bgIndex];
+            psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN);
+            double normValue = normSolution->data.F64[0];
+            double bgValue = normSolution->data.F64[1];
+            psFree(normMatrix);
+            psFree(normVector);
+            psFree(normSolution);
+#endif
+            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
+            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
+                psFree(stats);
                 psFree(sumMatrix);
                 psFree(sumVector);
+                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 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 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;                  \
+                }                                                       \
+            }
+            #define TOL 1.0e-5
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            double norm = fabs(solution->data.F64[normIndex]);  // Normalisation
+            double thresh = norm * TOL;                         // Threshold for low parameters
+            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", i);
+                    MASK_BASIS_1(i);
+                    mask1 = true;
+                }
+                if (fabs(param2) < thresh) {
+                    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", i);
+                        MASK_BASIS_1(i);
+                    } else {
+                        psTrace("psModules.imcombine", 7, "Parameter %d: 2 < 1\n", i);
+                        MASK_BASIS_2(i);
+                    }
+                }
+            }
+        }
+        calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
+                              sumMatrixX, sumVector1, sumVector2);
+                psFree(norms);
+                return false;
+            }
+            double normValue = stats->robustMedian;
+            psFree(stats);
+            fprintf(stderr, "Norm: %lf\n", normValue);
+            // Solve kernel components
+            for (int i = 0; i < numSolution2; i++) {
+                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
+                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1];
+                sumMatrix->data.F64[i][normIndex] = 0.0;
+                sumMatrix->data.F64[normIndex][i] = 0.0;
+                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
+                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
+            }
+            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
+            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
+            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
+            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
+            sumVector->data.F64[normIndex] = 0.0;
+            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+            solution->data.F64[normIndex] = normValue;
+        }
+#endif
 #ifdef TESTING
+        {
+            psFits *fits = psFitsOpen("sumMatrixFix.fits", "w");
+            psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
+            psFitsClose(fits);
+        }
+        {
+            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
+            psFits *fits = psFitsOpen("sumVectorFix.fits", "w");
+            for (int i = 0; i < numParams + numParams2; i++) {
+                image->data.F64[0][i] = sumVector->data.F64[i];
+            }
+            psFitsWriteImage(fits, NULL, image, 0, NULL);
+            psFree(image);
+            psFitsClose(fits);
+        }
+#endif
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
+        if (!solution) {
+            psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
+            psFree(sumMatrix);
+            psFree(sumVector);
+            return NULL;
+        }
+        psFree(sumMatrix1);
+        psFree(sumMatrix2);
+        psFree(sumMatrixX);
+        psFree(sumVector1);
+        psFree(sumVector2);
+        for (int i = 0; i < solution->n; i++) {
+            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
+        }
+#endif
         psFree(sumMatrix);
         psFree(sumVector);
+#ifdef TESTING
+        {
+            psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
+            psFits *fits = psFitsOpen("solnVector.fits", "w");
+            for (int i = 0; i < numParams + numParams2; i++) {
+                image->data.F64[0][i] = solution->data.F64[i];
+            }
+            psFitsWriteImage(fits, NULL, image, 0, NULL);
+            psFree(image);
+            psFitsClose(fits);
+        }
+#endif
+        psFree(norms);
         if (!kernels->solution1) {
+            kernels->solution1 = psVectorAlloc(numParams, PS_TYPE_F64);
+            kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64);
+            psVectorInit (kernels->solution1, 0.0);
+        }
         if (!kernels->solution2) {
+            kernels->solution2 = psVectorAlloc(numParams2, PS_TYPE_F64);
+        }
+        memcpy(kernels->solution1->data.F64, solution->data.F64, numParams * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        memcpy(kernels->solution2->data.F64, &solution->data.F64[numParams],
+               numParams2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+            kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
+            psVectorInit (kernels->solution2, 0.0);
+        }
+        // only update the solutions that we chose to calculate:
+        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_BG) {
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+            int numKernels = kernels->num;
+            for (int i = 0; i < numKernels * numSpatial; i++) {
+                // XXX fprintf (stderr, "keep\n");
+                kernels->solution1->data.F64[i] = solution->data.F64[i];
+                kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
+            }
+        }
+        memcpy(kernels->solution1->data.F64, solution->data.F64,
+               numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1],
+               numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
         psFree(solution);
 …
+     }
     pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
+    // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
     return true;
+}
+bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
+    // XXX measure some useful stats on the residuals
+    float sum = 0.0;
+    float peak = 0.0;
+    for (int y = - footprint; y <= footprint; y++) {
+        for (int x = - footprint; x <= footprint; x++) {
+            sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+            peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
+        }
+    }
+    // only count pixels with more than X% of the source flux
+    // calculate stdev(dflux)
+    float dflux1 = 0.0;
+    float dflux2 = 0.0;
+    int npix = 0;
+    float dmax = 0.0;
+    float dmin = 0.0;
+    for (int y = - footprint; y <= footprint; y++) {
+        for (int x = - footprint; x <= footprint; x++) {
+            float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+            if (dflux < 0.02*sum) continue;
+            dflux1 += residual->kernel[y][x];
+            dflux2 += PS_SQR(residual->kernel[y][x]);
+            dmax = PS_MAX(residual->kernel[y][x], dmax);
+            dmin = PS_MIN(residual->kernel[y][x], dmin);
+            npix ++;
+        }
+    }
+    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
+    if (!isfinite(sum))  return false;
+    if (!isfinite(dmax)) return false;
+    if (!isfinite(dmin)) return false;
+    if (!isfinite(peak)) return false;
+    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
+    psVectorAppend(fSigRes, sigma/sum);
+    psVectorAppend(fMaxRes, dmax/peak);
+    psVectorAppend(fMinRes, dmin/peak);
+    return true;
+}
 psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
                                            const pmSubtractionKernels *kernels)
+                                           pmSubtractionKernels *kernels)
+{
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
 …
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    // set up holding images for the visualization
+    pmSubtractionVisualShowFitInit (stamps);
+    psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    // we want to save the residual images for the 9 brightest stamps.
+    // identify the 9 brightest stamps
+    psVector *keepStamps  = psVectorAlloc(stamps->num, PS_TYPE_S32);
+    psVectorInit (keepStamps, 0);
+    {
+        psVector *flux  = psVectorAlloc(stamps->num, PS_TYPE_F32);
+        psVectorInit (flux, 0.0);
+        for (int i = 0; i < stamps->num; i++) {
+            pmSubtractionStamp *stamp = stamps->stamps->data[i];
+            if (!isfinite(stamp->flux)) continue;
+            flux->data.F32[i] = stamp->flux;
+        }
+        psVector *index = psVectorSortIndex(NULL, flux);
+        for (int i = 0; (i < stamps->num) && (i < 9); i++) {
+            int n = stamps->num - i - 1;
+            keepStamps->data.S32[index->data.S32[n]] = 1;
+            fprintf (stderr, "keeping %d (%d of %d)\n", index->data.S32[n], n, 9);
+        }
+        psFree (flux);
+        psFree (index);
+        // this function is called multiple times in the iteration, but
+        // we only know after the interation is done if we will try again.
+        // therefore we must save the sample each time, and blow away the old one
+        // if it exists.
+        psFree (kernels->sampleStamps);
+        kernels->sampleStamps = psArrayAllocEmpty(9);
+    }
+    psString log = psStringCopy("Deviations:\n");               // Log message with deviations
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
 …
                 for (int y = - footprint; y <= footprint; y++) {
                     for (int x = - footprint; x <= footprint; x++) {
                         residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
+                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
+                    }
+                }
+            }
+            // XXX visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
+                    residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm;
+                }
+            }
+                    residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x];
+                }
+            }
+            if (keepStamps->data.S32[i]) {
+                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+                psArrayAdd (kernels->sampleStamps, 9, sample);
+                psFree (sample);
+            }
+            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
         } else {
             // Dual convolution
 …
                 for (int y = - footprint; y <= footprint; y++) {
                     for (int x = - footprint; x <= footprint; x++) {
                         residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1;
+                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
+                    }
+                }
+            }
+            // XXX visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
+                    residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm;
+                }
+            }
+                    residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
+                }
+            }
+            if (keepStamps->data.S32[i]) {
+                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+                psArrayAdd (kernels->sampleStamps, 9, sample);
+                psFree (sample);
+            }
+            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
+        }
 …
         deviations->data.F32[i] = devNorm * deviation;
         psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
+                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
+                i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
+        psStringAppend(&log, "Stamp %d (%d,%d): %f\n",
+                       i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
         if (!isfinite(deviations->data.F32[i])) {
             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
             psTrace("psModules.imcombine", 5,
                     "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
                     i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+                    i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
             continue;
+        }
 …
+    }
+    psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log);
+    psFree(log);
+    // calculate and report the normalization and background for the image center
+    {
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
+        pmSubtractionVisualShowFit(norm);
+        pmSubtractionVisualPlotFit(kernels);
+        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+        psVectorStats (stats, fSigRes, NULL, NULL, 0);
+        kernels->fSigResMean = stats->robustMedian;
+        kernels->fSigResStdev = stats->robustStdev;
+        psStatsInit (stats);
+        psVectorStats (stats, fMaxRes, NULL, NULL, 0);
+        kernels->fMaxResMean = stats->robustMedian;
+        kernels->fMaxResStdev = stats->robustStdev;
+        psStatsInit (stats);
+        psVectorStats (stats, fMinRes, NULL, NULL, 0);
+        kernels->fMinResMean = stats->robustMedian;
+        kernels->fMinResStdev = stats->robustStdev;
+        // XXX save these values somewhere
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",
+                 kernels->fSigResMean, kernels->fSigResStdev,
+                 kernels->fMaxResMean, kernels->fMaxResStdev,
+                 kernels->fMinResMean, kernels->fMinResStdev);
+        psFree (fSigRes);
+        psFree (fMaxRes);
+        psFree (fMinRes);
+        psFree (stats);
+    }
     psFree(residual);
     psFree(polyValues);
     return deviations;
+}
+// we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w)
+psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) {
+    psAssert (w->n == U->numCols, "w and U dimensions do not match");
+    // wUt has dimensions transposed relative to Ut.
+    psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64);
+    psImageInit (wUt, 0.0);
+    for (int i = 0; i < wUt->numCols; i++) {
+        for (int j = 0; j < wUt->numRows; j++) {
+            if (!isfinite(w->data.F64[j])) continue;
+            if (w->data.F64[j] == 0.0) continue;
+            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
+        }
+    }
+    return wUt;
+}
+// XXX this is just standard matrix multiplication: use psMatrixMultiply?
+psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) {
+    psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match");
+    psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64);
+    for (int i = 0; i < Ainv->numCols; i++) {
+        for (int j = 0; j < Ainv->numRows; j++) {
+            double sum = 0.0;
+            for (int k = 0; k < V->numCols; k++) {
+                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
+            }
+            Ainv->data.F64[j][i] = sum;
+        }
+    }
+    return Ainv;
+}
+// we are supplied U, not Ut
+bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) {
+    psAssert (U->numRows == B->n, "U and B dimensions do not match");
+    UtB[0] = psVectorRecycle (UtB[0], U->numCols, PS_TYPE_F64);
+    for (int i = 0; i < U->numCols; i++) {
+        double sum = 0.0;
+        for (int j = 0; j < U->numRows; j++) {
+            sum += B->data.F64[j] * U->data.F64[j][i];
+        }
+        UtB[0]->data.F64[i] = sum;
+    }
+    return true;
+}
+// w is diagonal
+bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) {
+    psAssert (w->n == UtB->n, "w and UtB dimensions do not match");
+    // wUt has dimensions transposed relative to Ut.
+    wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64);
+    psVectorInit (wUtB[0], 0.0);
+    for (int i = 0; i < w->n; i++) {
+        if (!isfinite(w->data.F64[i])) continue;
+        if (w->data.F64[i] == 0.0) continue;
+        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
+    }
+    return true;
+}
+// this is basically matrix * vector
+bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) {
+    psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match");
+    VwUtB[0] = psVectorRecycle (*VwUtB, V->numRows, PS_TYPE_F64);
+    for (int j = 0; j < V->numRows; j++) {
+        double sum = 0.0;
+        for (int i = 0; i < V->numCols; i++) {
+            sum += V->data.F64[j][i] * wUtB->data.F64[i];
+        }
+        VwUtB[0]->data.F64[j] = sum;
+    }
+    return true;
+}
+// this is basically matrix * vector
+bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) {
+    psAssert (A->numCols == x->n, "A and x dimensions do not match");
+    B[0] = psVectorRecycle (*B, A->numRows, PS_TYPE_F64);
+    for (int j = 0; j < A->numRows; j++) {
+        double sum = 0.0;
+        for (int i = 0; i < A->numCols; i++) {
+            sum += A->data.F64[j][i] * x->data.F64[i];
+        }
+        B[0]->data.F64[j] = sum;
+    }
+    return true;
+}
+// this is basically Vector * vector
+bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y) {
+    psAssert (x->n == y->n, "x and y dimensions do not match");
+    double sum = 0.0;
+    for (int i = 0; i < x->n; i++) {
+        sum += x->data.F64[i] * y->data.F64[i];
+    }
+    *value = sum;
+    return true;
+}
+bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
+    int footprint = stamps->footprint; // Half-size of stamps
+    double sum = 0.0;
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i];
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+        psKernel *weight = NULL;
+        psKernel *window = NULL;
+        psKernel *input = NULL;
+#ifdef USE_WEIGHT
+        weight = stamp->weight;
+#endif
+#ifdef USE_WINDOW
+        window = stamps->window;
+#endif
+        switch (kernels->mode) {
+            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
+          case PM_SUBTRACTION_MODE_1:
+            input = stamp->image2;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            input = stamp->image1;
+            break;
+          default:
+            psAbort ("programming error");
+        }
+        for (int y = - footprint; y <= footprint; y++) {
+            for (int x = - footprint; x <= footprint; x++) {
+                double in = input->kernel[y][x];
+                double value = in*in;
+                if (weight) {
+                    float wtVal = weight->kernel[y][x];
+                    value *= wtVal;
+                }
+                if (window) {
+                    float  winVal = window->kernel[y][x];
+                    value *= winVal;
+                }
+                sum += value;
+            }
+        }
+    }
+    *y2 = sum;
+    return true;
+}
+double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
+    int footprint = stamps->footprint; // Half-size of stamps
+    int numKernels = kernels->num;      // Number of kernels
+    double sum = 0.0;
+    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    psImageInit(residual->image, 0.0);
+    psImage *polyValues = NULL;         // Polynomial values
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i];
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+        psKernel *weight = NULL;
+        psKernel *window = NULL;
+        psKernel *target = NULL;
+        psKernel *source = NULL;
+        psArray *convolutions = NULL;
+#ifdef USE_WEIGHT
+        weight = stamp->weight;
+#endif
+#ifdef USE_WINDOW
+        window = stamps->window;
+#endif
+        switch (kernels->mode) {
+            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
+          case PM_SUBTRACTION_MODE_1:
+            target = stamp->image2;
+            source = stamp->image1;
+            convolutions = stamp->convolutions1;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            target = stamp->image1;
+            source = stamp->image2;
+            convolutions = stamp->convolutions2;
+            break;
+          default:
+            psAbort ("programming error");
+        }
+        // Calculate coefficients of the kernel basis functions
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+        psImageInit(residual->image, 0.0);
+        for (int j = 0; j < numKernels; j++) {
+            psKernel *convolution = convolutions->data[j]; // Convolution
+            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
+                }
+            }
+        }
+        for (int y = - footprint; y <= footprint; y++) {
+            for (int x = - footprint; x <= footprint; x++) {
+                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
+                double value = PS_SQR(resid);
+                if (weight) {
+                    float wtVal = weight->kernel[y][x];
+                    value *= wtVal;
+                }
+                if (window) {
+                    float  winVal = window->kernel[y][x];
+                    value *= winVal;
+                }
+                sum += value;
+            }
+        }
+    }
+    psFree (polyValues);
+    psFree (residual);
+    return sum;
+}
+bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) {
+    for (int i = 0; i < w->n; i++) {
+        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
+    }
+    return true;
+}
+// we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w)
+psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) {
+    psAssert (w->n == V->numCols, "w and U dimensions do not match");
+    psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
+    psImageInit (Vn, 0.0);
+    // generate Vn = V * w^{-1}
+    for (int j = 0; j < Vn->numRows; j++) {
+        for (int i = 0; i < Vn->numCols; i++) {
+            if (!isfinite(w->data.F64[i])) continue;
+            if (w->data.F64[i] == 0.0) continue;
+            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
+        }
+    }
+    psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
+    psImageInit (Xvar, 0.0);
+    // generate Xvar = Vn * Vn^T
+    for (int j = 0; j < Vn->numRows; j++) {
+        for (int i = 0; i < Vn->numCols; i++) {
+            double sum = 0.0;
+            for (int k = 0; k < Vn->numCols; k++) {
+                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
+            }
+            Xvar->data.F64[j][i] = sum;
+        }
+    }
+    return Xvar;
+}
+// I get confused by the index values between the image vs matrix usage:  In terms
+// of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix
+// multiplication is: A_k,j * B_i,k = C_i,j
+bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) {
+    psFits *fits = psFitsOpen(filename, "w");
+    psFitsWriteImage(fits, header, image, 0, NULL);
+    psFitsClose(fits);
+    return true;
+}
+bool psVectorWriteFile (char *filename, const psVector *vector) {
+    FILE *f = fopen (filename, "w");
+    int fd = fileno(f);
+    p_psVectorPrint (fd, vector, "unnamed");
+    fclose (f);
+    return true;
+}
+# if 0
+#ifdef TESTING
+        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
+        psVectorWriteFile ("B.dat", sumVector);
+#endif
+# define SVD_ANALYSIS 0
+# define COEFF_SIG 0.0
+# define SVD_TOL 0.0
+        // Use SVD to determine the kernel coeffs (and validate)
+        if (SVD_ANALYSIS) {
+            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
+            // sumMatrix * x = sumVector.
+            // we can use any standard matrix inversion to solve this.  However, the basis
+            // functions in general have substantial correlation, so that the solution may be
+            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
+            // system of equations may be statistically ill-conditioned.  Noise in the image
+            // will drive insignificant, but correlated, terms in the solution.  To avoid these
+            // problems, we can use SVD to identify numerically unconstrained values and to
+            // avoid statistically badly determined value.
+            // A = sumMatrix, B = sumVector
+            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
+            // x = V (1/w) (U^T B)
+            // \sigma_x = sqrt(diag(A^{-1}))
+            // solve for x and A^{-1} to get x & dx
+            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
+            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
+            // If I use the SVD trick to re-condition the matrix, I need to break out the
+            // kernel and normalization terms from the background term.
+            // XXX is this true?  or was this due to an error in the analysis?
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+            // now pull out the kernel elements into their own square matrix
+            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
+            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
+            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
+                if (ix == bgIndex) continue;
+                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
+                    if (iy == bgIndex) continue;
+                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
+                    ky++;
+                }
+                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
+                kx++;
+            }
+            psImage *U = NULL;
+            psImage *V = NULL;
+            psVector *w = NULL;
+            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
+                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
+                return NULL;
+            }
+            // calculate A_inverse:
+            // Ainv = V * w * U^T
+            psImage *wUt  = p_pmSubSolve_wUt (w, U);
+            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
+            psImage *Xvar = NULL;
+            psFree (wUt);
+# ifdef TESTING
+            // kernel terms:
+            for (int i = 0; i < w->n; i++) {
+                fprintf (stderr, "w: %f\n", w->data.F64[i]);
+            }
+# endif
+            // loop over w adding in more and more of the values until chisquare is no longer
+            // dropping significantly.
+            // XXX this does not seem to work very well: we seem to need all terms even for
+            // simple cases...
+            psVector *Xsvd = NULL;
+            {
+                psVector *Ax = NULL;
+                psVector *UtB = NULL;
+                psVector *wUtB = NULL;
+                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
+                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
+                psVectorInit (wMask, 1); // start by masking everything
+                double chiSquareLast = NAN;
+                int maxWeight = 0;
+                double Axx, Bx, y2;
+                // XXX this is an attempt to exclude insignificant modes.
+                // it was not successful with the ISIS kernel set: removing even
+                // the least significant mode leaves additional ringing / noise
+                // because the terms are so coupled.
+                for (int k = 0; false && (k < w->n); k++) {
+                    // unmask the k-th weight
+                    wMask->data.U8[k] = 0;
+                    p_pmSubSolve_SetWeights(wApply, w, wMask);
+                    // solve for x:
+                    // x = V * w * (U^T * B)
+                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
+                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
+                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
+                    // chi-square for this system of equations:
+                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
+                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
+                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
+                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
+                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
+                    p_pmSubSolve_y2 (&y2, kernels, stamps);
+                    // apparently, this works (compare with the brute force value below
+                    double chiSquare = Axx - 2.0*Bx + y2;
+                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
+                    chiSquareLast = chiSquare;
+                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
+                    if (k && !maxWeight && (deltaChi < 1.0)) {
+                        maxWeight = k;
+                    }
+                }
+                // keep all terms or we get extra ringing
+                maxWeight = w->n;
+                psVectorInit (wMask, 1);
+                for (int k = 0; k < maxWeight; k++) {
+                    wMask->data.U8[k] = 0;
+                }
+                p_pmSubSolve_SetWeights(wApply, w, wMask);
+                // solve for x:
+                // x = V * w * (U^T * B)
+                p_pmSubSolve_UtB (&UtB, U, kernelVector);
+                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
+                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
+                // chi-square for this system of equations:
+                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
+                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
+                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
+                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
+                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
+                p_pmSubSolve_y2 (&y2, kernels, stamps);
+                // apparently, this works (compare with the brute force value below
+                double chiSquare = Axx - 2.0*Bx + y2;
+                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
+                // re-calculate A^{-1} to get new variances:
+                // Ainv = V * w * U^T
+                // XXX since we keep all terms, this is identical to Ainv
+                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
+                Xvar = p_pmSubSolve_VwUt (V, wUt);
+                psFree (wUt);
+                psFree (Ax);
+                psFree (UtB);
+                psFree (wUtB);
+                psFree (wApply);
+                psFree (wMask);
+            }
+            // copy the kernel solutions to the full solution vector:
+            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
+                if (ix == bgIndex) {
+                    solution->data.F64[ix] = 0;
+                    solutionErr->data.F64[ix] = 0.001;
+                    continue;
+                }
+                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
+                solution->data.F64[ix] = Xsvd->data.F64[kx];
+                kx++;
+            }
+            psFree (kernelMatrix);
+            psFree (kernelVector);
+            psFree (U);
+            psFree (V);
+            psFree (w);
+            psFree (Ainv);
+            psFree (Xsvd);
+        } else {
+            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
+            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
+            if (!luMatrix) {
+                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
+                psFree(solution);
+                psFree(sumVector);
+                psFree(sumMatrix);
+                psFree(luMatrix);
+                psFree(permutation);
+                return NULL;
+            }
+            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
+            psFree(luMatrix);
+            psFree(permutation);
+            if (!solution) {
+                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
+                psFree(solution);
+                psFree(sumVector);
+                psFree(sumMatrix);
+                return NULL;
+            }
+            // XXX LUD does not provide A^{-1}?  fake the error for now
+            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+            for (int ix = 0; ix < sumVector->n; ix++) {
+                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
+            }
+        }
+        if (!kernels->solution1) {
+            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
+            psVectorInit (kernels->solution1, 0.0);
+        }
+        // only update the solutions that we chose to calculate:
+        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_BG) {
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+            int numKernels = kernels->num;
+            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
+            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+            for (int i = 0; i < numKernels * numPoly; i++) {
+                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
+                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
+                    // XXX fprintf (stderr, "drop\n");
+                    kernels->solution1->data.F64[i] = 0.0;
+                } else {
+                    // XXX fprintf (stderr, "keep\n");
+                    kernels->solution1->data.F64[i] = solution->data.F64[i];
+                }
+            }
+        }
+        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
+        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
+        psFree(solution);
+        psFree(sumVector);
+        psFree(sumMatrix);
+# endif
+#ifdef TESTING
+              // XXX double-check for NAN in data:
+                for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
+                    for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
+                        if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
+                            fprintf (stderr, "WARNING: NAN in matrix\n");
+                        }
+                    }
+                }
+                for (int ix = 0; ix < stamp->vector->n; ix++) {
+                    if (!isfinite(stamp->vector->data.F64[ix])) {
+                        fprintf (stderr, "WARNING: NAN in vector\n");
+                    }
+                }
+#endif
+#ifdef TESTING
+        for (int ix = 0; ix < sumVector->n; ix++) {
+            if (!isfinite(sumVector->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector\n");
+            }
+        }
+#endif
+#ifdef TESTING
+        for (int ix = 0; ix < sumVector->n; ix++) {
+            if (!isfinite(sumVector->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector\n");
+            }
+        }
+        {
+            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
+            psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
+            psFree(inverse);
+        }
+        {
+            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
+            psImage *Xt = psMatrixTranspose(NULL, X);
+            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
+            psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
+            psFree(X);
+            psFree(Xt);
+            psFree(XtX);
+        }
+#endif

trunk/psModules/src/imcombine/pmSubtractionEquation.h

-              r19299
+              r26893
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionKernels.h"
+typedef enum {
+    PM_SUBTRACTION_EQUATION_NONE    = 0x00,
+    PM_SUBTRACTION_EQUATION_NORM    = 0x01,
+    PM_SUBTRACTION_EQUATION_BG      = 0x02,
+    PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
+    PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
+} pmSubtractionEquationCalculationMode;
 /// Execute a thread job to calculate the least-squares equation for a stamp
 …
 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps
                                          const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                         int index ///< Index of stamp
+                                    );
+                                         int index, ///< Index of stamp
+                                         const pmSubtractionEquationCalculationMode mode
+    );
 /// Calculate the least-squares equation to match the image quality
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
+                                    const pmSubtractionKernels *kernels ///< Kernel parameters
+                                    );
+                                    const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                    const pmSubtractionEquationCalculationMode mode
+    );
 /// Solve the least-squares equation to match the image quality
 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
+                                const pmSubtractionStampList *stamps ///< Stamps
+                                const pmSubtractionStampList *stamps, ///< Stamps
+                                const pmSubtractionEquationCalculationMode mode
     );
 /// Calculate deviations
 psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps
                                            const pmSubtractionKernels *kernels ///< Kernel parameters
+                                           pmSubtractionKernels *kernels ///< Kernel parameters
     );

trunk/psModules/src/imcombine/pmSubtractionIO.c

-              r23378
+              r26893
         while ((item = psMetadataGetAndIncrement(iter))) {
             assert(item->type == PS_DATA_UNKNOWN);
-            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
-            // images by scans.
             pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
-            kernel->numCols = ro->image->numCols;
-            kernel->numRows = ro->image->numRows;
             kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
+        }
 …
                          kernel->bgOrder);
         psMetadataAddS32(row, PS_LIST_TAIL, NAME_MODE,  0, "Matching mode (enum)", kernel->mode);
-        psMetadataAddS32(row, PS_LIST_TAIL, NAME_COLS,  0, "Number of columns", kernel->numCols);
-        psMetadataAddS32(row, PS_LIST_TAIL, NAME_ROWS,  0, "Number of rows", kernel->numRows);
         if (kernel->mode == PM_SUBTRACTION_MODE_1 || kernel->mode == PM_SUBTRACTION_MODE_2) {
             psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL1, 0, "Solution vector 1", kernel->solution1);
 …
         TABLE_LOOKUP(int, S32, bg,      NAME_BG);
-        TABLE_LOOKUP(int, S32, numCols, NAME_COLS);
-        TABLE_LOOKUP(int, S32, numRows, NAME_ROWS);
         TABLE_LOOKUP(float, F32, mean,      NAME_MEAN);
 …
         TABLE_LOOKUP(int,   S32, numStamps, NAME_NUMSTAMPS);
+        pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode);
+        kernels->numCols = numCols;
+        kernels->numRows = numRows;
+        pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, *region, mode);
         kernels->mean = mean;
         kernels->rms = rms;

trunk/psModules/src/imcombine/pmSubtractionKernels.c

-              r26156
+              r26893
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
+#include "pmSubtractionHermitian.h"
+#include "pmSubtractionDeconvolve.h"
+#include "pmSubtractionVisual.h"
 #define RINGS_BUFFER 10                 // Buffer size for RINGS data
 // Free function for pmSubtractionKernels
 …
     psFree(kernels->solution1);
     psFree(kernels->solution2);
+    psFree(kernels->sampleStamps);
+}
+// Free function for pmSubtractionPreCalcKernel
+static void pmSubtractionKernelPreCalcFree(pmSubtractionKernelPreCalc *kernel)
+{
+    psFree(kernel->xKernel);
+    psFree(kernel->yKernel);
+    psFree(kernel->kernel);
+    psFree(kernel->uCoords);
+    psFree(kernel->vCoords);
+    psFree(kernel->poly);
+}
 …
 // Generate 1D convolution kernel for ISIS
 static psVector *subtractionKernelISIS(float sigma, // Gaussian width
+psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width
                                        int order, // Polynomial order
                                        int size // Kernel half-size
 …
     for (int i = 0, x = -size; x <= size; i++, x++) {
         kernel->data.F32[i] = norm * power(x, order) * expf(expNorm * PS_SQR(x));
+    }
+    return kernel;
+}
+// Generate 1D convolution kernel for HERM (normalized for 2D)
+psVector *pmSubtractionKernelHERM(float sigma, // Gaussian width
+                                       int order, // Polynomial order
+                                       int size // Kernel half-size
+    )
+{
+    int fullSize = 2 * size + 1;        // Full size of kernel
+    psVector *kernel = psVectorAlloc(fullSize, PS_TYPE_F32); // Kernel to return
+    // for now, we are only allowing equal orders and sigmas in X and Y
+    float nf = exp(lgamma(order + 1));
+    float norm = 1.0 / sqrt(nf*sigma*sqrt(M_2_PI));
+    for (int i = 0, x = -size; x <= size; i++, x++) {
+        float xf = x / sigma;
+        float z = -0.25*xf*xf;
+        kernel->data.F32[i] = norm * p_pmSubtractionHermitianPolynomial(xf, order) * exp(z);
+    }
+    return kernel;
+}
+// Generate 1D convolution kernel for HERM (normalized for 2D)
+psKernel *pmSubtractionKernelHERM_RADIAL(float sigma, // Gaussian width
+                                         int order, // Polynomial order
+                                         int size // Kernel half-size
+    )
+{
+    psKernel *kernel = psKernelAlloc(-size, size, -size, size); // 2D Kernel
+    // for now, we are only allowing equal orders and sigmas in X and Y
+    float nf = exp(lgamma(order + 1));
+    float norm = 1.0 / sqrt(nf*sigma*sqrt(M_2_PI));
+    // generate 2D radial hermitian
+    for (int v = -size; v <= size; v++) {
+        for (int u = -size; u <= size; u++) {
+            float r = hypot(u, v) / sigma;
+            float z = -0.25*r*r;
+            kernel->kernel[v][u] = norm * p_pmSubtractionHermitianPolynomial(r, order) * exp(z);
+        }
+    }
 …
             kernels->preCalc->data[index] = NULL;
             kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
+            psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
+        }
 …
+}
+bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
+                                         int index, int size, int uOrder, int vOrder, float fwhm,
+                                         bool AlardLuptonStyle, bool forceZeroNull)
+{
+    // we have 4 cases here:
+    // 1) for odd functions, normalize the kernel by the maximum swing / Npix
+    // 2) for even functions, normalize the kernel to unity
+    // 3) for alard-lupton style normalization, subtract 1 from the 0,0 pixel for all even functions
+    // 4) for deconvolved hermitians, subtract 1 from the 0,0 pixel for the 0,0 function(s)
+    // Calculate moments
+    double penalty = 0.0;                   // Moment, for penalty
+    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
+    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+    for (int v = -size; v <= size; v++) {
+        for (int u = -size; u <= size; u++) {
+            double value = preCalc->kernel->kernel[v][u];
+            double value2 = PS_SQR(value);
+            sum += value;
+            sum2 += value2;
+            penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+        }
+    }
+#if 0
+    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
+#endif
+    bool zeroNull = false;              // Zero out using the null position?
+    float scale2D = NAN;                // Scaling for 2-D kernels
+    if (AlardLuptonStyle) {
+        if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+            // Even functions: normalise to unit sum and subtract null pixel so that sum is zero
+            scale2D = 1.0 / fabs(sum);
+            zeroNull = true;
+        } else {
+            // Odd functions: choose normalisation so that parameters have about the same strength as for even
+            // functions, no subtraction of null pixel because the sum is already (near) zero
+            scale2D = 1.0 / sqrt(sum2);
+            zeroNull = false;
+        }
+    }
+    if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) {
+        zeroNull = true;
+    }
+    if (forceZeroNull) {
+        // Force rescaling and subtraction of null pixel even though the order doesn't indicate it's even
+        scale2D = 1.0 / fabs(sum);
+        zeroNull = true;
+    }
+    if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
+        // Odd function
+        scale2D = 1.0 / sqrt(sum2);
+    }
+    float scale1D = sqrtf(scale2D);     // Scaling for 1-D kernels
+    if (preCalc->xKernel) {
+        psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+    }
+    if (preCalc->yKernel) {
+        psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+    }
+    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
+    penalty *= 1.0 / sum2;
+    if (zeroNull) {
+        preCalc->kernel->kernel[0][0] -= 1.0;
+    }
+#if 0
+    {
+        double sum = 0.0;   // Sum of kernel component
+        float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+        for (int v = -size; v <= size; v++) {
+            for (int u = -size; u <= size; u++) {
+                sum += preCalc->kernel->kernel[v][u];
+                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
+                max = PS_MAX(preCalc->kernel->kernel[v][u], max);
+            }
+        }
+        fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
+    }
+#endif
+    kernels->widths->data.F32[index] = fwhm;
+    kernels->u->data.S32[index] = uOrder;
+    kernels->v->data.S32[index] = vOrder;
+    if (kernels->preCalc->data[index]) {
+        psFree(kernels->preCalc->data[index]);
+    }
+    kernels->preCalc->data[index] = preCalc;
+    kernels->penalties->data.F32[index] = kernels->penalty * penalty;
+    psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
+    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
+    return true;
+}
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
                                                     const psVector *fwhms, const psVector *orders,
                                                     float penalty, pmSubtractionMode mode)
+{
     PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
     PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL);
     PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
     PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL);
+                                                    const psVector *fwhmsIN, const psVector *ordersIN,
+                                                    float penalty, psRegion bounds, pmSubtractionMode mode)
+{
+    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL);
     PS_ASSERT_INT_POSITIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+    // check the requested fwhm values: any values <= 0.0 should be dropped
+    psVector *fwhms  = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32);
+    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
+    for (int i = 0; i < fwhmsIN->n; i++) {
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
+    }
     int numGaussians = fwhms->n;       // Number of Gaussians
 …
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
                                                               spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
 …
     // Set the kernel parameters
-    int fullSize = 2 * size + 1;        // Full size of kernels
     for (int i = 0, index = 0; i < numGaussians; i++) {
         float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
 …
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
+                psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
+                psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel
+                psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel
+                psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
+                // Calculate moments
+                double moment = 0.0;    // Moment, for penalty
+                for (int v = -size, y = 0; v <= size; v++, y++) {
+                    for (int u = -size, x = 0; u <= size; u++, x++) {
+                        double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
+                        kernel->kernel[v][u] = value;
+                        moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
+                    }
+                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
+                // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false);
+            }
+        }
+    }
+    return kernels;
+}
+pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder,
+                                                      const psVector *fwhmsIN, const psVector *ordersIN,
+                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
+{
+    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL);
+    PS_ASSERT_INT_POSITIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+    // check the requested fwhm values: any values <= 0.0 should be dropped
+    psVector *fwhms  = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32);
+    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
+    for (int i = 0; i < fwhmsIN->n; i++) {
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
+    }
+    int numGaussians = fwhms->n;       // Number of Gaussians
+    int num = 0;                        // Number of basis functions
+    psString params = NULL;             // List of parameters
+    for (int i = 0; i < numGaussians; i++) {
+        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
+        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
+        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
+        num += (11 - gaussOrder - 1);   // include all higher order radial terms
+    }
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size,
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
+    psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS_RADIAL kernel: %s,%d --> %d elements", params, spatialOrder, num);
+    psFree(params);
+    // Set the kernel parameters
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
+        // Iterate over (u,v) order
+        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
+                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
+            }
+        }
+        for (int order = orders->data.S32[i] + 1; order < 11; order ++, index ++) {
+            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
+            pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values
+            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true);
+        }
+    }
+    return kernels;
+}
+pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
+                                               const psVector *fwhmsIN, const psVector *ordersIN,
+                                               float penalty, psRegion bounds, pmSubtractionMode mode)
+{
+    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL);
+    PS_ASSERT_INT_POSITIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+    // check the requested fwhm values: any values <= 0.0 should be dropped
+    psVector *fwhms  = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32);
+    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
+    for (int i = 0; i < fwhmsIN->n; i++) {
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
+    }
+    int numGaussians = fwhms->n;       // Number of Gaussians
+    int num = 0;                        // Number of basis functions
+    psString params = NULL;             // List of parameters
+    for (int i = 0; i < numGaussians; i++) {
+        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
+        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
+        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
+    }
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size,
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
+    psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements",
+             params, spatialOrder, num);
+    psFree(params);
+    // Set the kernel parameters
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
+        // Iterate over (u,v) order
+        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
+                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
+            }
+        }
+    }
+    return kernels;
+}
+pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
+                                                      const psVector *fwhmsIN, const psVector *ordersIN,
+                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
+{
+    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL);
+    PS_ASSERT_INT_POSITIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+    // check the requested fwhm values: any values <= 0.0 should be dropped
+    psVector *fwhms  = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32);
+    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
+    for (int i = 0; i < fwhmsIN->n; i++) {
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
+    }
+    int numGaussians = fwhms->n;       // Number of Gaussians
+    int num = 0;                        // Number of basis functions
+    psString params = NULL;             // List of parameters
+    for (int i = 0; i < numGaussians; i++) {
+        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
+        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
+        num += PS_SQR(gaussOrder + 1);
+    }
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size,
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
+    psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
+    psFree(params);
+    // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest)
+    // generate the Gaussian deconvolution kernel
+    # define DECONV_SIGMA 1.6
+    psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECONV_SIGMA);
+# if 1
+    psArray *deconKernels = psArrayAllocEmpty(100);
+# endif
+    // Set the kernel parameters
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
+        // Iterate over (u,v) order
+        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i]; vOrder++, index++) {
+                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
+                // save the generated 2D kernel as the target, deconvolve it by Gaussian, replacing the generated 2D kernel
+                psKernel *kernelTarget = preCalc->kernel;
+                preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
+                // XXX do we use Alard-Lupton normalization (last param true) or not?
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
+                // XXXX test demo that deconvolved kernel is valid
+# if 1
+                psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
+                psArrayAdd (deconKernels, 100, kernelConv);
+                psFree (kernelConv);
+                if (!uOrder && !vOrder){
+                    pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
+                }
+                // Normalise sum of kernel component to unity for even functions
+                if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+                    double sum = 0.0;   // Sum of kernel component
+                    for (int v = 0; v < fullSize; v++) {
+                        for (int u = 0; u < fullSize; u++) {
+                            sum += xKernel->data.F32[u] * yKernel->data.F32[v];
+                        }
+                    }
+                    sum = 1.0 / sqrt(sum);
+                    psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
+                    kernel->kernel[0][0] -= 1.0;
+                    moment *= PS_SQR(sum);
+# endif
+            }
+        }
+    }
+# if 1
+    psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32);
+    for (int i = 0; i < deconKernels->n; i++) {
+        for (int j = 0; j <= i; j++) {
+            psImage *t1 = deconKernels->data[i];
+            psImage *t2 = deconKernels->data[j];
+            double sum = 0.0;
+            for (int iy = 0; iy < t1->numRows; iy++) {
+                for (int ix = 0; ix < t1->numCols; ix++) {
+                    sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix];
+                }
-#if 0
-                double sum = 0.0;   // Sum of kernel component
-                for (int v = -size; v <= size; v++) {
-                    for (int u = -size; u <= size; u++) {
-                        sum += kernel->kernel[v][u];
+                    }
+                }
-                fprintf(stderr, "%d sum: %lf\n", index, sum);
-#endif
-                kernels->widths->data.F32[index] = fwhms->data.F32[i];
-                kernels->u->data.S32[index] = uOrder;
-                kernels->v->data.S32[index] = vOrder;
-                if (kernels->preCalc->data[index]) {
-                    psFree(kernels->preCalc->data[index]);
+                }
-                kernels->preCalc->data[index] = preCalc;
-                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
-                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
-                        fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
+            }
+        }
+    }
+            dot->data.F32[j][i] = sum;
+            dot->data.F32[i][j] = sum;
+        }
+    }
+    pmSubtractionVisualShowSubtraction (dot, NULL, NULL);
+    psFree (dot);
+    psFree (deconKernels);
+# endif
     return kernels;
 …
 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
                                                 int size, int spatialOrder, float penalty,
+                                                int size, int spatialOrder, float penalty, psRegion bounds,
                                                 pmSubtractionMode mode)
+{
 …
     kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
     kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
+    kernels->uStop = NULL;
+    kernels->vStop = NULL;
+    kernels->xMin = bounds.x0;
+    kernels->xMax = bounds.x1;
+    kernels->yMin = bounds.y0;
+    kernels->yMax = bounds.y1;
     kernels->preCalc = psArrayAlloc(numBasisFunctions);
     kernels->penalty = penalty;
     kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
-    kernels->uStop = NULL;
-    kernels->vStop = NULL;
     kernels->size = size;
     kernels->inner = 0;
 …
     kernels->bgOrder = 0;
     kernels->mode = mode;
-    kernels->numCols = 0;
-    kernels->numRows = 0;
     kernels->solution1 = NULL;
     kernels->solution2 = NULL;
+    kernels->mean = NAN;
+    kernels->rms = NAN;
+    kernels->numStamps = 0;
+    kernels->sampleStamps = NULL;
+    kernels->fSigResMean  = NAN;
+    kernels->fSigResStdev = NAN;
+    kernels->fMaxResMean  = NAN;
+    kernels->fMaxResStdev = NAN;
+    kernels->fMinResMean  = NAN;
+    kernels->fMinResStdev = NAN;
     return kernels;
+}
+pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
+pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(pmSubtractionKernelsType type, int uOrder, int vOrder, int size, float sigma) {
+    pmSubtractionKernelPreCalc *preCalc = psAlloc(sizeof(pmSubtractionKernelPreCalc)); // Kernels, to return
+    psMemSetDeallocator(preCalc, (psFreeFunc)pmSubtractionKernelPreCalcFree);
+    // 1D kernel realizations:
+    switch (type) {
+      case PM_SUBTRACTION_KERNEL_ISIS:
+        preCalc->xKernel = pmSubtractionKernelISIS(sigma, uOrder, size);
+        preCalc->yKernel = pmSubtractionKernelISIS(sigma, vOrder, size);
+        preCalc->uCoords = NULL;
+        preCalc->vCoords = NULL;
+        preCalc->poly    = NULL;
+        break;
+      case PM_SUBTRACTION_KERNEL_HERM:
+        preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
+        preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
+        preCalc->uCoords = NULL;
+        preCalc->vCoords = NULL;
+        preCalc->poly    = NULL;
+        break;
+      case PM_SUBTRACTION_KERNEL_RINGS:
+        // the RINGS kernel uses the uCoords, vCoords, and poly elements of the structure
+        // we allocate these vectors here, but leave the kernel generation to the main function
+        preCalc->xKernel = NULL;
+        preCalc->yKernel = NULL;
+        preCalc->kernel  = NULL;
+        preCalc->uCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // u coords
+        preCalc->vCoords = psVectorAllocEmpty(size, PS_TYPE_S32); // v coords
+        preCalc->poly    = psVectorAllocEmpty(size, PS_TYPE_F32); // Polynomial
+        return preCalc;
+      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
+        preCalc->kernel  = pmSubtractionKernelHERM_RADIAL(sigma, uOrder, size);
+        preCalc->xKernel = NULL;
+        preCalc->yKernel = NULL;
+        preCalc->uCoords = NULL;
+        preCalc->vCoords = NULL;
+        preCalc->poly    = NULL;
+        return preCalc;
+      default:
+        psAbort("programming error: invalid type for PreCalc kernel");
+    }
+    preCalc->kernel = psKernelAlloc(-size, size, -size, size); // 2D Kernel
+    // generate 2D kernel from 1D realizations
+    for (int v = -size, y = 0; v <= size; v++, y++) {
+        for (int u = -size, x = 0; u <= size; u++, x++) {
+            preCalc->kernel->kernel[v][u] = preCalc->xKernel->data.F32[x] * preCalc->yKernel->data.F32[y]; // Value of kernel
+        }
+    }
+    return preCalc;
+}
+pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, psRegion bounds,
                                                pmSubtractionMode mode)
+{
 …
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
                                                               spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
 …
 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
                                                const psVector *fwhms, const psVector *orders,
                                                float penalty, pmSubtractionMode mode)
+                                               float penalty, psRegion bounds, pmSubtractionMode mode)
+{
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
                                                                   penalty, mode); // Kernels
+                                                                  penalty, bounds, mode); // Kernels
     if (!kernels) {
         return NULL;
 …
 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
                                                float penalty, pmSubtractionMode mode)
+                                               float penalty, psRegion bounds, pmSubtractionMode mode)
+{
     PS_ASSERT_INT_POSITIVE(size, NULL);
 …
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
                                                               spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
     psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
 …
 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
                                                 pmSubtractionMode mode)
+                                                psRegion bounds, pmSubtractionMode mode)
+{
     PS_ASSERT_INT_POSITIVE(size, NULL);
 …
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
                                                               spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
     psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
 …
+}
 // Grid United with Normal Kernel
+// Grid United with Normal Kernel [description: GUNK=ISIS(...)+POIS(...)]
 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
                                                const psVector *orders, int inner, float penalty,
                                                pmSubtractionMode mode)
+                                               psRegion bounds, pmSubtractionMode mode)
+{
     PS_ASSERT_INT_POSITIVE(size, NULL);
 …
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
                                                                   penalty, mode); // Kernels
+                                                                  penalty, bounds, mode); // Kernels
     kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
     psStringPrepend(&kernels->description, "GUNK=");
 …
 // RINGS --- just what it says
 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
                                                 float penalty, pmSubtractionMode mode)
+                                                float penalty, psRegion bounds, pmSubtractionMode mode)
+{
     PS_ASSERT_INT_POSITIVE(size, NULL);
 …
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
                                                               spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
     psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
 …
             for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) {
+                psArray *data = psArrayAlloc(3); // Container for data
+                psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords
+                psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords
+                psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial
+                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
                 double moment = 0.0;    // Moment, for penalty
                 if (i == 0) {
                     // Central pixel is easy
+                    uCoords->data.S32[0] = vCoords->data.S32[0] = 0;
+                    poly->data.F32[0] = 1.0;
+                    uCoords->n = vCoords->n = poly->n = 1;
+                    preCalc->uCoords->data.S32[0] = 0;
+                    preCalc->vCoords->data.S32[0] = 0;
+                    preCalc->poly->data.F32[0] = 1.0;
+                    preCalc->uCoords->n = 1;
+                    preCalc->vCoords->n = 1;
+                    preCalc->poly->n = 1;
                     radiusLast = 0;
                     moment = 0.0;
 …
                                 float polyVal = uPoly * vPoly; // Value of polynomial
                                 if (polyVal != 0) { // No point adding it otherwise
                                     uCoords->data.S32[j] = u;
                                     vCoords->data.S32[j] = v;
                                     poly->data.F32[j] = polyVal;
+                                    preCalc->uCoords->data.S32[j] = u;
+                                    preCalc->vCoords->data.S32[j] = v;
+                                    preCalc->poly->data.F32[j] = polyVal;
                                     norm += polyVal;
                                     moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v));
                                     psVectorExtend(uCoords, RINGS_BUFFER, 1);
                                     psVectorExtend(vCoords, RINGS_BUFFER, 1);
                                     psVectorExtend(poly, RINGS_BUFFER, 1);
+                                    moment += PS_SQR(polyVal) * PS_SQR(PS_SQR(u) + PS_SQR(v));
+                                    psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1);
+                                    psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1);
+                                    psVectorExtend(preCalc->poly, RINGS_BUFFER, 1);
                                     psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n",
                                             u, v, poly->data.F32[j]);
+                                            u, v, preCalc->poly->data.F32[j]);
                                     j++;
+                                }
 …
                     // Normalise kernel component to unit sum
                     if (uOrder % 2 == 0 && vOrder % 2 == 0) {
                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
+                        psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
                         // Add subtraction of 0,0 component to preserve photometric scaling
                         uCoords->data.S32[j] = 0;
                         vCoords->data.S32[j] = 0;
                         poly->data.F32[j] = -1.0;
                         psVectorExtend(uCoords, RINGS_BUFFER, 1);
                         psVectorExtend(vCoords, RINGS_BUFFER, 1);
                         psVectorExtend(poly, RINGS_BUFFER, 1);
+                        preCalc->uCoords->data.S32[j] = 0;
+                        preCalc->vCoords->data.S32[j] = 0;
+                        preCalc->poly->data.F32[j] = -1.0;
+                        psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1);
+                        psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1);
+                        psVectorExtend(preCalc->poly, RINGS_BUFFER, 1);
                     } else {
                         norm = powf(size, uOrder) * powf(size, vOrder);
                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
+                        psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
+                    }
                     moment /= norm;
+                    moment /= PS_SQR(norm);
+                }
                 psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", uCoords->n);
                 kernels->preCalc->data[index] = data;
+                psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", preCalc->uCoords->n);
+                kernels->preCalc->data[index] = preCalc;
                 kernels->u->data.S32[index] = uOrder;
                 kernels->v->data.S32[index] = vOrder;
                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
+                if (!isfinite(kernels->penalties->data.F32[index])) {
+                    psAbort ("invalid penalty");
+                }
                 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,
 …
 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
                                                    const psVector *fwhms, const psVector *orders, int inner,
                                                    int binning, int ringsOrder, float penalty,
+                                                   int binning, int ringsOrder, float penalty, psRegion bounds,
                                                    pmSubtractionMode mode)
+{
     switch (type) {
       case PM_SUBTRACTION_KERNEL_POIS:
         return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode);
+        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_ISIS:
+        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
+        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
+      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
+        return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
+      case PM_SUBTRACTION_KERNEL_HERM:
+        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
+      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
+        return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_SPAM:
         return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
+        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_FRIES:
         return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode);
+        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_GUNK:
         return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode);
+        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_RINGS:
         return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode);
+        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode);
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
 …
 pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder,
                                                           pmSubtractionMode mode)
+                                                          psRegion bounds, pmSubtractionMode mode)
+{
     PS_ASSERT_STRING_NON_EMPTY(description, NULL);
 …
     float penalty = 0.0;                // Penalty for wideness
+    if (strncmp(description, "ISIS", 4) == 0) {
+        // XXX Support for GUNK
+        if (strstr(description, "+POIS")) {
+            type = PM_SUBTRACTION_KERNEL_GUNK;
+            psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description);
+        } else {
+            type = PM_SUBTRACTION_KERNEL_ISIS;
+            char *ptr = (char*)description + 5;    // Eat "ISIS("
+            PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
+            // Count the number of Gaussians
+            int numGauss = 0;
+            for (char *string = ptr; string; string = strchr(string + 1, '(')) {
+                numGauss++;
+            }
+            fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
+            orders = psVectorAlloc(numGauss, PS_TYPE_S32);
+            for (int i = 0; i < numGauss; i++) {
+                ptr++;                  // Eat the '('
+                PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
+                PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)"
+            }
+            ptr++;                      // Eat ','
+            PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
+            penalty = parseStringFloat(ptr);
+        }
+    } else if (strncmp(description, "RINGS", 5) == 0) {
+        type = PM_SUBTRACTION_KERNEL_RINGS;
+        char *ptr = (char*)description + 6;
+    // currently known descriptions:
+    // ISIS(...), ISIS_RADIAL(...), HERM(...), DECONV_HERM(...), POIS(...), SPAM(...),
+    // FRIES(...), GUNK=ISIS(...)+POIS(...), RINGS(...),
+    // the descriptive name is the set of characters before the (
+    type = pmSubtractionKernelsTypeFromString (description);
+    char *ptr = strchr(description, '(') + 1;
+    psAssert (ptr, "description is missing kernel parameters");
+    switch (type) {
+      case PM_SUBTRACTION_KERNEL_ISIS:
+      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
+      case PM_SUBTRACTION_KERNEL_HERM:
+      case PM_SUBTRACTION_KERNEL_DECONV_HERM:
+        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
+        // Count the number of Gaussians
+        int numGauss = 0;
+        for (char *string = ptr; string; string = strchr(string + 1, '(')) {
+            numGauss++;
+        }
+        fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
+        orders = psVectorAlloc(numGauss, PS_TYPE_S32);
+        for (int i = 0; i < numGauss; i++) {
+            ptr++;                                                               // Eat the '('
+            PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
+            PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt);  // Eat "3)"
+        }
+        ptr++;                      // Eat ','
+        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
+        penalty = parseStringFloat(ptr);
+        break;
+      case PM_SUBTRACTION_KERNEL_RINGS:
         PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
         PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt);
 …
         PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
         PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
+    } else {
+        psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported.");
+    }
+    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders,
+                                        inner, binning, ringsOrder, penalty, mode);
+}
+        break;
+      default:
+        psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
+    }
+    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning,
+                                        ringsOrder, penalty, bounds, mode);
+}
+// the input string can either be just the name or the description string.  Currently known
+// descriptions: ISIS(...), ISIS_RADIAL(...), HERM(...), DECONV_HERM(...), POIS(...),
+// SPAM(...), FRIES(...), GUNK=ISIS(...)+POIS(...), RINGS(...),
 pmSubtractionKernelsType pmSubtractionKernelsTypeFromString(const char *type)
+{
+    if (strcasecmp(type, "POIS") == 0) {
+    // for a bare name (ISIS, HERM), use the full string length.
+    // otherwise, use the length up to the first '('
+    int nameLength = strlen(type);
+    char *ptr = strchr(type, '(');
+    if (ptr) {
+        nameLength = ptr - type;
+    }
+    if (strncasecmp(type, "POIS", nameLength) == 0) {
         return PM_SUBTRACTION_KERNEL_POIS;
+    }
     if (strcasecmp(type, "ISIS") == 0) {
+    if (strncasecmp(type, "ISIS", nameLength) == 0) {
         return PM_SUBTRACTION_KERNEL_ISIS;
+    }
+    if (strcasecmp(type, "SPAM") == 0) {
+    if (strncasecmp(type, "ISIS_RADIAL", nameLength) == 0) {
+        return PM_SUBTRACTION_KERNEL_ISIS_RADIAL;
+    }
+    if (strncasecmp(type, "HERM", nameLength) == 0) {
+        return PM_SUBTRACTION_KERNEL_HERM;
+    }
+    if (strncasecmp(type, "DECONV_HERM", nameLength) == 0) {
+        return PM_SUBTRACTION_KERNEL_DECONV_HERM;
+    }
+    if (strncasecmp(type, "SPAM", nameLength) == 0) {
         return PM_SUBTRACTION_KERNEL_SPAM;
+    }
     if (strcasecmp(type, "FRIES") == 0) {
+    if (strncasecmp(type, "FRIES", nameLength) == 0) {
         return PM_SUBTRACTION_KERNEL_FRIES;
+    }
     if (strcasecmp(type, "GUNK") == 0) {
+    if (strncasecmp(type, "GUNK", nameLength) == 0) {
         return PM_SUBTRACTION_KERNEL_GUNK;
+    }
+    if (strcasecmp(type, "RINGS") == 0) {
+    // note that GUNK has a somewhat different description
+    if (strncasecmp(type, "GUNK=ISIS", nameLength) == 0) {
+        return PM_SUBTRACTION_KERNEL_GUNK;
+    }
+    if (strncasecmp(type, "RINGS", nameLength) == 0) {
         return PM_SUBTRACTION_KERNEL_RINGS;
+    }
 …
     out->bgOrder = in->bgOrder;
     out->mode = in->mode;
+    out->numCols = in->numCols;
+    out->numRows = in->numRows;
+    out->xMin = in->xMin;
+    out->xMax = in->xMax;
+    out->yMin = in->yMin;
+    out->yMax = in->yMax;
     out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
     out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
+    out->sampleStamps = psMemIncrRefCounter(in->sampleStamps);
     return out;

trunk/psModules/src/imcombine/pmSubtractionKernels.h

-              r26157
+              r26893
     PM_SUBTRACTION_KERNEL_POIS,         ///< Pan-STARRS Optimal Image Subtraction --- delta functions
     PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
+    PM_SUBTRACTION_KERNEL_ISIS_RADIAL,  ///< ISIS + higher-order radial Hermitians
+    PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
+    PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
     PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
     PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
 …
     pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
     psString description;               ///< Description of the kernel parameters
+    int xMin, xMax, yMin, yMax;         ///< Bounds of image (for normalisation)
     long num;                           ///< Number of kernel components (not including the spatial ones)
     psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
     psVector *widths;                   ///< Gaussian FWHMs (ISIS)
+    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM)
+    psVector *widths;                   ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM)
     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
+    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM)
     float penalty;                      ///< Penalty for wideness
     psVector *penalties;                ///< Penalty for each kernel component
 …
     int bgOrder;                        ///< The order for the background fitting
     pmSubtractionMode mode;             ///< Mode for subtraction
-    int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
     psVector *solution1, *solution2;    ///< Solution for the PSF matching
     // Quality information
     float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
     int numStamps;                      ///< Number of good stamps
+    float fSigResMean;                  ///< mean fractional stdev of residuals
+    float fSigResStdev;                 ///< stdev of fractional stdev of residuals
+    float fMaxResMean;                  ///< mean fractional positive swing in residuals
+    float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
+    float fMinResMean;                  ///< mean fractional negative swing in residuals
+    float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
+    psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
 } pmSubtractionKernels;
+// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures
+typedef struct {
+    psVector *uCoords;                  // used by RINGS
+    psVector *vCoords;                  // used by RINGS
+    psVector *poly;                     // used by RINGS
+    psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
+    psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
+    psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
+} pmSubtractionKernelPreCalc;
 // Assertion to check pmSubtractionKernels
 …
         PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
     } \
+    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_ISIS_RADIAL) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
+    } \
+    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_HERM) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
+    } \
+    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_DECONV_HERM) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
+    } \
     if ((KERNELS)->uStop || (KERNELS)->vStop) { \
         PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \
 …
+}
+// Generate 1D convolution kernel for ISIS
+psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width
+                                       int order, // Polynomial order
+                                       int size // Kernel half-size
+    );
+// Generate 1D convolution kernel for HERM (normalized for 2D)
+psVector *pmSubtractionKernelHERM(float sigma, // Gaussian width
+                                       int order, // Polynomial order
+                                       int size // Kernel half-size
+    );
 /// Generate a delta-function grid for subtraction kernels (like the POIS kernel)
 bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, ///< The subtraction kernels to append to
 …
                                                 int spatialOrder, ///< Order of spatial variations
                                                 float penalty, ///< Penalty for wideness
+                                                psRegion bounds,       ///< Bounds for validity
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
+/// Allocator for pre-calculated kernel data structure
+pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(
+    pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated)
+    int uOrder,                    ///< order in x-direction
+    int vOrder,                    ///< order in x-direction
+    int size,                      ///< Half-size of the kernel
+    float sigma                    ///< sigma of gaussian kernel
+    );
 /// Generate POIS kernels
 …
                                                int spatialOrder, ///< Order of spatial variations
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 …
                                                     const psVector *orders, ///< Polynomial order of gaussians
                                                     float penalty, ///< Penalty for wideness
+                                                    psRegion bounds,       ///< Bounds for validity
                                                     pmSubtractionMode mode ///< Mode for subtraction
     );
 …
                                                const psVector *orders, ///< Polynomial order of gaussians
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
                                                );
+/// Generate ISIS + RADIAL_HERM kernels
+pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, ///< Half-size of the kernel
+                                                      int spatialOrder, ///< Order of spatial variations
+                                                      const psVector *fwhms, ///< Gaussian FWHMs
+                                                      const psVector *orders, ///< Polynomial order of gaussians
+                                                      float penalty, ///< Penalty for wideness
+                                                      psRegion bounds,       ///< Bounds for validity
+                                                      pmSubtractionMode mode ///< Mode for subtraction
+                                               );
+/// Generate HERM kernels
+pmSubtractionKernels *pmSubtractionKernelsHERM(int size, ///< Half-size of the kernel
+                                               int spatialOrder, ///< Order of spatial variations
+                                               const psVector *fwhms, ///< Gaussian FWHMs
+                                               const psVector *orders, ///< order of hermitian polynomials
+                                               float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
+                                               pmSubtractionMode mode ///< Mode for subtraction
+                                               );
+/// Generate DECONV_HERM kernels
+pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, ///< Half-size of the kernel
+                                                      int spatialOrder, ///< Order of spatial variations
+                                                      const psVector *fwhms, ///< Gaussian FWHMs
+                                                      const psVector *orders, ///< order of hermitian polynomials
+                                                      float penalty, ///< Penalty for wideness
+                                                      psRegion bounds,       ///< Bounds for validity
+                                                      pmSubtractionMode mode ///< Mode for subtraction
+    );
 /// Generate SPAM kernels
 …
                                                int binning, ///< Kernel binning factor
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 …
                                                 int inner, ///< Inner radius to preserve unbinned
                                                 float penalty, ///< Penalty for wideness
+                                                psRegion bounds,       ///< Bounds for validity
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
 …
                                                int inner, ///< Inner radius containing grid of delta functions
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 …
                                                 int ringsOrder, ///< Polynomial order
                                                 float penalty, ///< Penalty for wideness
+                                                psRegion bounds,       ///< Bounds for validity
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
 …
                                                    int ringsOrder, ///< Polynomial order for RINGS
                                                    float penalty, ///< Penalty for wideness
+                                                   psRegion bounds,       ///< Bounds for validity
                                                    pmSubtractionMode mode ///< Mode for subtraction
     );
 …
     const char *description,            ///< Description of kernel
     int bgOrder,                        ///< Polynomial order for background fitting
+    psRegion bounds,                    ///< Bounds for validity
     pmSubtractionMode mode              ///< Mode for subtraction
     );

trunk/psModules/src/imcombine/pmSubtractionMask.c

-              r24257
+              r26893
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
+psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psImageMaskType maskVal,
+                           int size, int footprint, float badFrac, pmSubtractionMode mode)
+{
+    PS_ASSERT_IMAGE_NON_NULL(mask1, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_IMAGE_MASK, NULL);
+    if (mask2) {
+        PS_ASSERT_IMAGE_NON_NULL(mask2, NULL);
+        PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_IMAGE_MASK, NULL);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL);
+    }
+psImage *pmSubtractionMask(psRegion *bounds, const pmReadout *ro1, const pmReadout *ro2,
+                           psImageMaskType maskVal, int size, int footprint, float badFrac,
+                           pmSubtractionMode mode)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro1, NULL);
+    PM_ASSERT_READOUT_IMAGE(ro1, NULL);
+    PM_ASSERT_READOUT_MASK(ro1, NULL);
+    PM_ASSERT_READOUT_NON_NULL(ro2, NULL);
+    PM_ASSERT_READOUT_IMAGE(ro2, NULL);
+    PM_ASSERT_READOUT_MASK(ro2, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL);
     PS_ASSERT_INT_NONNEGATIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
 …
+    }
     int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images
+    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of the images
     // Dereference inputs for convenience
+    psImageMaskType **data1 = mask1->data.PS_TYPE_IMAGE_MASK_DATA;
+    psImageMaskType **data2 = NULL;
+    if (mask2) {
+        data2 = mask2->data.PS_TYPE_IMAGE_MASK_DATA;
+    }
+    psF32 **imageData1 = ro1->image->data.F32, **imageData2 = ro2->image->data.F32;
+    psImageMaskType **maskData1 = ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA,
+        **maskData2 = ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA;
     // First, a pass through to determine the fraction of bad pixels
+    if (isfinite(badFrac) && badFrac != 1.0) {
+    if (bounds || (isfinite(badFrac) && badFrac != 1.0)) {
+        int xMin = numCols, xMax = 0, yMin = numRows, yMax = 0; // Bounds of good pixels
         int numBad = 0;                 // Number of bad pixels
         for (int y = 0; y < numRows; y++) {
             for (int x = 0; x < numCols; x++) {
                 if (data1[y][x] & maskVal) {
+                if ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x])) {
                     numBad++;
                     continue;
+                }
                 if (data2 && data2[y][x] & maskVal) {
+                if ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x])) {
                     numBad++;
+                    continue;
+                }
+            }
+        }
+        if (numBad > badFrac * numCols * numRows) {
+                xMin = PS_MIN(xMin, x);
+                xMax = PS_MAX(xMax, x);
+                yMin = PS_MIN(yMin, y);
+                yMax = PS_MAX(yMax, y);
+            }
+        }
+        if (bounds) {
+            bounds->x0 = xMin;
+            bounds->x1 = xMax;
+            bounds->y0 = yMin;
+            bounds->y1 = yMax;
+        }
+        if (isfinite(badFrac) && badFrac != 1.0 && numBad > badFrac * numCols * numRows) {
             psError(PM_ERR_SMALL_AREA, true,
                     "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n",
 …
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
             if (data1[y][x] & maskVal) {
+            if (maskData1[y][x] & maskVal) {
                 maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1;
+            }
             if (data2 && data2[y][x] & maskVal) {
+            if (maskData2[y][x] & maskVal) {
                 maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2;
+            }

trunk/psModules/src/imcombine/pmSubtractionMask.h

-              r23851
+              r26893
 /// Generate a mask for use in the subtraction process
+psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved)
+                           const psImage *inMask, ///< Mask for the input image, or NULL
+                           psImageMaskType maskVal, ///< Value to mask out
+                           int size, ///< Half-size of the kernel (pmSubtractionKernels.size)
+                           int footprint, ///< Half-size of the kernel footprint
+                           float badFrac, ///< Maximum fraction of bad input pixels to accept
+                           pmSubtractionMode mode  ///< Subtraction mode
+psImage *pmSubtractionMask(
+    psRegion *bounds,                   ///< Bounds of valid pixels (or NULL), returned
+    const pmReadout *ro1,               ///< Readout 1
+    const pmReadout *ro2,               ///< Readout 2
+    psImageMaskType maskVal,            ///< Value to mask out
+    int size,                           ///< Half-size of the kernel (pmSubtractionKernels.size)
+    int footprint,                      ///< Half-size of the kernel footprint
+    float badFrac,                      ///< Maximum fraction of bad input pixels to accept
+    pmSubtractionMode mode              ///< Subtraction mode
     );

trunk/psModules/src/imcombine/pmSubtractionMatch.c

-              r26035
+              r26893
 static bool useFFT = true;              // Do convolutions using FFT
+# define SEPARATE 0
+# if (SEPARATE)
+# define SUBMODE PM_SUBTRACTION_EQUATION_NORM
+# else
+# define SUBMODE PM_SUBTRACTION_EQUATION_ALL
+# endif
 //#define TESTING
 …
                                  const psImage *subMask, // Mask for subtraction, or NULL
                                  psImage *variance,  // Variance map
                                  const psRegion *region, // Region of interest, or NULL
+                                 const psRegion *region, // Region of interest
                                  float thresh1,  // Threshold for stamp finding on readout 1
                                  float thresh2,  // Threshold for stamp finding on readout 2
                                  float stampSpacing, // Spacing between stamps
+                                 float normFrac,     // Fraction of flux in window for normalisation window
                                  float sysError,     // Relative systematic error in images
+                                 float skyError,     // Relative systematic error in images
                                  int size,         // Kernel half-size
                                  int footprint,     // Convolution footprint for stamps
 …
+    )
+{
+    PS_ASSERT_PTR_NON_NULL(stamps, false);
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
+    PS_ASSERT_IMAGE_NON_NULL(subMask, false);
+    PS_ASSERT_IMAGE_NON_NULL(variance, false);
+    PS_ASSERT_PTR_NON_NULL(region, false);
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
 …
     *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
                                       size, footprint, stampSpacing, sysError, mode);
+                                      size, footprint, stampSpacing, normFrac, sysError, skyError, mode);
     if (!*stamps) {
         psError(psErrorCodeLast(), false, "Unable to find stamps.");
 …
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
     if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
 …
                                   const pmReadout *ro1, const pmReadout *ro2, // Input images
                                   int stride, // Size for convolution patches
+                                  float normFrac,           // Fraction of window for normalisation window
                                   float sysError,           // Systematic error in images
+                                  float skyError,           // Systematic error in images
                                   float kernelError, // Systematic error in kernel
+                                  float covarFrac,   // Fraction for kernel truncation before covariance
                                   psImageMaskType maskVal, // Value to mask for input
                                   psImageMaskType maskBad, // Mask for output bad pixels
 …
     PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false);
     PS_ASSERT_INT_NONNEGATIVE(stride, false);
+    if (isfinite(normFrac)) {
+        PS_ASSERT_FLOAT_LARGER_THAN(normFrac, 0.0, false);
+        PS_ASSERT_FLOAT_LESS_THAN(normFrac, 1.0, false);
+    }
     if (isfinite(sysError)) {
         PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(sysError, 0.0, false);
         PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
+    }
+    if (isfinite(sysError)) {
+        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.0, false);
+    }
     if (isfinite(kernelError)) {
         PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false);
         PS_ASSERT_FLOAT_LESS_THAN(kernelError, 1.0, false);
+    }
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(covarFrac, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN(covarFrac, 1.0, false);
     // Don't care about maskVal
     // Don't care about maskBad
 …
+}
+bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {
+    if (!readout) return true;
+    psImage *image = readout->image;
+    psImage *mask  = readout->mask;
+    psImage *variance = readout->variance;
+    for (int y = 0; y < image->numRows; y++) {
+        for (int x = 0; x < image->numCols; x++) {
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;
+            bool valid = false;
+            valid = isfinite(image->data.F32[y][x]);
+            if (variance) {
+                valid &= isfinite(variance->data.F32[y][x]);
+            }
+            if (valid) continue;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;
+        }
+    }
+    return true;
+}
 bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
                                psMetadata *analysis, int stride, float kernelError,
+                               psMetadata *analysis, int stride, float kernelError, float covarFrac,
                                psImageMaskType maskVal, psImageMaskType maskBad, psImageMaskType maskPoor,
                                float poorFrac, float badFrac)
 …
+    }
     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, kernelError,
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, NAN, kernelError, covarFrac,
                                maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
         psFree(kernels);
 …
+    }
+    psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
+    pmSubtractionMaskInvalid(ro1, maskVal);
+    pmSubtractionMaskInvalid(ro2, maskVal);
+    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
+    psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0,
                                          badFrac, mode); // Subtraction mask
     if (!subMask) {
 …
         if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
                                    kernelError, region, kernel, true, useFFT)) {
+                                   kernelError, covarFrac, region, kernel, true, useFFT)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
             psFree(outAnalysis);
 …
                         int inner, int ringsOrder, int binning, float penalty,
                         bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
                         int iter, float rej, float sysError, float kernelError, psImageMaskType maskVal,
                         psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac,
                         float badFrac, pmSubtractionMode subMode)
+                        int iter, float rej, float normFrac, float sysError, float skyError,
+                        float kernelError, float covarFrac, psImageMaskType maskVal, psImageMaskType maskBad,
+                        psImageMaskType maskPoor, float poorFrac, float badFrac, pmSubtractionMode subMode)
+{
     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, kernelError,
                                maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
+    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, normFrac, sysError, skyError, kernelError,
+                               covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
         return false;
+    }
 …
     // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
     psImage *subMask = NULL;            // Mask for subtraction
     psRegion *region = NULL;            // Iso-kernel region
+    psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region
     psString regionString = NULL;       // String for region
     pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
 …
     memCheck("start");
+    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
+                                badFrac, subMode);
+    pmSubtractionMaskInvalid(ro1, maskVal);
+    pmSubtractionMaskInvalid(ro2, maskVal);
+    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
+    subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode);
     if (!subMask) {
         psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
 …
     // Get region of interest
     int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
     float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
+    float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions
     if (isfinite(regionSize) && regionSize != 0.0) {
+        xRegions = numCols / regionSize + 1;
+        yRegions = numRows / regionSize + 1;
+        xRegionSize = (float)numCols / (float)xRegions;
+        yRegionSize = (float)numRows / (float)yRegions;
+        region = psRegionAlloc(NAN, NAN, NAN, NAN);
+    }
+        xRegions = (bounds.x1 - bounds.x0) / regionSize + 1;
+        yRegions = (bounds.y1 - bounds.y0) / regionSize + 1;
+        xRegionSize = (float)(bounds.x1 - bounds.x0) / (float)xRegions;
+        yRegionSize = (float)(bounds.y1 - bounds.y0) / (float)yRegions;
+    } else {
+        xRegionSize = bounds.x1 - bounds.x0;
+        yRegionSize = bounds.y1 - bounds.y0;
+    }
+    // General background subtraction and measurement of stamp threshold
     float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
+    {
         psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
+        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
         if (ro1) {
+            psStatsInit(bg);
             if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
 …
                 goto MATCH_ERROR;
+            }
+            stampThresh1 = bg->robustMedian + threshold * bg->robustStdev;
+            stampThresh1 = threshold * bg->robustStdev;
+            psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+        }
         if (ro2) {
+            psStatsInit(bg);
             if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics.");
 …
                 goto MATCH_ERROR;
+            }
+            stampThresh2 = bg->robustMedian + threshold * bg->robustStdev;
+        }
+            stampThresh2 = threshold * bg->robustStdev;
+            psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32));
+       }
         psFree(bg);
+    }
+    // Just in case the iso-kernel region doesn't cover the entire image...
+    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
+        subMode == PM_SUBTRACTION_MODE_DUAL) {
+        if (!conv1->image) {
+            conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        psImageInit(conv1->image, NAN);
+        if (ro1->variance) {
+            if (!conv1->variance) {
+                conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
+            psImageInit(conv1->variance, NAN);
+        }
+        if (subMask) {
+            if (!conv1->mask) {
+                conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            }
+            psImageInit(conv1->mask, maskBad);
+        }
+    }
+    if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
+        subMode == PM_SUBTRACTION_MODE_DUAL) {
+        if (!conv2->image) {
+            conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        psImageInit(conv2->image, NAN);
+        if (ro2->variance) {
+            if (!conv2->variance) {
+                conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
+            psImageInit(conv2->variance, NAN);
+        }
+        if (subMask) {
+            if (!conv2->mask) {
+                conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            }
+            psImageInit(conv2->mask, maskBad);
+        }
+    }
     // Iterate over iso-kernel regions
 …
             psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
                     j * xRegions + i + 1, xRegions * yRegions);
             if (region) {
                 *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
                                       (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
                 psFree(regionString);
                 regionString = psRegionToString(*region);
                 psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
                         regionString, numCols, numRows);
+            }
+            *region = psRegionSet(bounds.x0 + (int)(i * xRegionSize),
+                                  bounds.x0 + (int)((i + 1) * xRegionSize),
+                                  bounds.y0 + (int)(j * yRegionSize),
+                                  bounds.y0 + (int)((j + 1) * yRegionSize));
+            psFree(regionString);
+            regionString = psRegionToString(*region);
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
+                    regionString, numCols, numRows);
             if (stampsName && strlen(stampsName) > 0) {
                 stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
+                                                        footprint, stampSpacing, sysError, subMode);
+                                                        footprint, stampSpacing, normFrac,
+                                                        sysError, skyError, subMode);
             } else if (sources) {
                 stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
+                                                           footprint, stampSpacing, sysError, subMode);
+                                                           footprint, stampSpacing, normFrac,
+                                                           sysError, skyError, subMode);
+            }
             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
             // doesn't matter.
             if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
                                       stampSpacing, sysError, size, footprint, subMode)) {
+            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
+                                      stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
                 goto MATCH_ERROR;
+            }
+            // generate the window function from the set of stamps
+            if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to get stamp window.");
+                goto MATCH_ERROR;
+            }
             // Define kernel basis functions
             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
+                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
+                                                          stamps, footprint, optThreshold, penalty, subMode);
+                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
+                                                          optFWHMs, optOrder, stamps, footprint,
+                                                          optThreshold, penalty, bounds, subMode);
                 if (!kernels) {
                     psErrorClear();
 …
                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
+                                                       inner, binning, ringsOrder, penalty, subMode);
+                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
+                // pmSubtractionVisualShowKernels(kernels);
+            }
 …
                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
+                                          stampThresh1, stampThresh2, stampSpacing, sysError,
+                                          size, footprint, subMode)) {
+                    goto MATCH_ERROR;
+                }
+                psTrace("psModules.imcombine", 3, "Calculating equation...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+                                          stampThresh1, stampThresh2, stampSpacing, normFrac,
+                                          sysError, skyError, size, footprint, subMode)) {
+                    goto MATCH_ERROR;
+                }
+                // generate the window function from the set of stamps
+                if (!pmSubtractionStampsGetWindow(stamps, size)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to get stamps window.");
+                    goto MATCH_ERROR;
+                }
+                // XXX step 1: calculate normalization
+                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
+                }
+                psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+                memCheck("  solve equation");
+# if (SEPARATE)
+                // set USED -> CALCULATE
+                pmSubtractionStampsResetStatus (stamps);
+                // XXX step 2: calculate kernel parameters
+                psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
                 memCheck("  calculate equation");
+                psTrace("psModules.imcombine", 3, "Solving equation...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps)) {
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
+                }
                 memCheck("  solve equation");
+# endif
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
 …
+            }
+            // if we hit the max number of iterations and we have rejected stamps, re-solve
             if (numRejected > 0) {
+                psTrace("psModules.imcombine", 3, "Solving equation...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps)) {
+                // XXX step 1: calculate normalization
+                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
+                }
+                // solve normalization
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+# if (SEPARATE)
+                // set USED -> CALCULATE
+                pmSubtractionStampsResetStatus (stamps);
+                // XXX step 2: calculate kernel parameters
+                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+                // solve kernel parameters
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+                memCheck("  solve equation");
+# endif
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
 …
             psTrace("psModules.imcombine", 2, "Convolving...\n");
             if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac,
                                        kernelError, region, kernels, true, useFFT)) {
+                                       kernelError, covarFrac, region, kernels, true, useFFT)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
                 goto MATCH_ERROR;
 …
     assert(kernels);
     psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description);
     if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
         return false;
+    }
     psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description);
     if (!pmSubtractionSolveEquation(kernels, stamps)) {
+    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
+    if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
         return false;
+    }
+# if (SEPARATE)
+    // set USED -> CALCULATE
+    pmSubtractionStampsResetStatus (stamps);
+    psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+    psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description);
+    if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+# endif
     psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
 …
     if (numRejected > 0) {
         // Allow re-fit with reduced stamps set
         psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
         if (!pmSubtractionSolveEquation(kernels, stamps)) {
+        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
             return false;
+        }
+        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
+        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+            return false;
+        }
         psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
+# if (SEPARATE)
+        // set USED -> CALCULATE
+        pmSubtractionStampsResetStatus (stamps);
+        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
+        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
+# endif
         psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
         if (!deviations) {
 …
+}
+bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
+                              float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax)
+{
+    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
+    PS_ASSERT_PTR_NON_NULL(stampSize, false);
+    PS_ASSERT_VECTOR_NON_NULL(widths, false);
+    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
+//    float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
+    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
+    if (isfinite(scaleMin) && scale < scaleMin) {
+        scale = scaleMin;
+    }
+    if (isfinite(scaleMax) && scale > scaleMax) {
+        scale = scaleMax;
+    }
+    for (int i = 0; i < widths->n; i++) {
+        widths->data.F32[i] *= scale;
+    }
+    *kernelSize = *kernelSize * scale + 0.5;
+    *stampSize = *stampSize * scale + 0.5;
+    psLogMsg("psModules.imcombine", PS_LOG_INFO,
+             "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
+    return true;
+}

trunk/psModules/src/imcombine/pmSubtractionMatch.h

-              r26035
+              r26893
                         int iter,       ///< Rejection iterations
                         float rej,      ///< Rejection threshold
+                        float normFrac, ///< Fraction of flux in window for normalisation window
                         float sysError, ///< Relative systematic error in images
+                        float skyError, ///< Relative systematic error in images
                         float kernelError, ///< Relative systematic error in kernel
+                        float covarFrac,   ///< Fraction for kernel truncation before covariance calculation
                         psImageMaskType maskVal, ///< Value to mask for input
                         psImageMaskType maskBad, ///< Mask for output bad pixels
 …
                                int stride, ///< Size for convolution patches
                                float kernelError, ///< Relative systematic error in kernel
+                               float covarFrac,   ///< Fraction for kernel truncation before covariance calc.
                                psImageMaskType maskVal, ///< Value to mask for input
                                psImageMaskType maskBad, ///< Mask for output bad pixels
 …
     );
+/// Scale subtraction parameters according to the FWHMs of the inputs
+bool pmSubtractionParamsScale(
+    int *kernelSize,                    ///< Half-size of the kernel
+    int *stampSize,                     ///< Half-size of the stamp (footprint)
+    psVector *widths,                   ///< ISIS widths
+    float fwhm1, float fwhm2,           ///< FWHMs for inputs
+    float scaleRef,                     ///< Reference width for scaling
+    float scaleMin,                     ///< Minimum scaling ratio, or NAN
+    float scaleMax                      ///< Maximum scaling ratio, or NAN
+    );
 #endif

trunk/psModules/src/imcombine/pmSubtractionParams.c

-              r26035
+              r26893
                                                       int spatialOrder, const psVector *fwhms, int maxOrder,
                                                       const pmSubtractionStampList *stamps, int footprint,
+                                                      float tolerance, float penalty, pmSubtractionMode mode)
+                                                      float tolerance, float penalty, psRegion bounds,
+                                                      pmSubtractionMode mode)
+{
     if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
 …
     psVectorInit(orders, maxOrder);
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
                                                                   penalty, mode); // Kernels
+                                                                  penalty, bounds, mode); // Kernels
     psFree(orders);
     psFree(kernels->description);
 …
     // Maintain photometric scaling
     if (type == PM_SUBTRACTION_KERNEL_ISIS) {
+        psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
+        // XXX in r26035, this code was just wrong.  we had:
+        // psKernel *subtract = kernels->preCalc->data[0]
+        // but, kernels->preCalc was an array of psArray, not an array of kernels.  It is now
+        // an array of pmSubtractionKernelPreCalc.
+        pmSubtractionKernelPreCalc *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
         for (int i = 1; i < newSize; i++) {
             if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
                 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
                 psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+                pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest
+                psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "-", subtract->kernel->image);
+            }
+        }
 …
         for (int i = 0; i < newSize; i++) {
             if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
                 psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
                 kernel->kernel[0][0] -= 1.0;
+                pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest
+                preCalc->kernel->kernel[0][0] -= 1.0;
+            }
+        }

trunk/psModules/src/imcombine/pmSubtractionParams.h

-              r18287
+              r26893
                                                       float tolerance, ///< Maximum difference in chi^2
                                                       float penalty, ///< Penalty for wideness
+                                                      pmSubtractionMode mode // Mode for subtraction
+                                                      psRegion bounds,       ///< Bounds of validity
+                                                      pmSubtractionMode mode ///< Mode for subtraction
     );

trunk/psModules/src/imcombine/pmSubtractionStamps.c

-              r26046
+              r26893
     psFree(list->y);
     psFree(list->flux);
+    psFree(list->window);
+}
 …
     psFree(stamp->convolutions1);
     psFree(stamp->convolutions2);
+    psFree(stamp->matrix1);
+    psFree(stamp->matrix2);
+    psFree(stamp->matrixX);
+    psFree(stamp->vector1);
+    psFree(stamp->vector2);
+    psFree(stamp->matrix);
+    psFree(stamp->vector);
+}
 …
 // Search a region for a suitable stamp
 bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return
+bool stampSearch(float *xStamp, float *yStamp, // Coordinates of stamp, to return
                  float *fluxStamp, // Flux of stamp, to return
                  const psImage *image1, const psImage *image2, // Images to search
 …
     *fluxStamp = -INFINITY;             // Flux of best stamp
+    // fprintf (stderr, "xMin, xMax: %d, %d -> ", xMin, xMax);
     // Ensure we're not going to go outside the bounds of the image
     xMin = PS_MAX(border, xMin);
 …
     yMax = PS_MIN(numRows - border - 1, yMax);
+    if (xMax < xMin) return false;
+    if (yMax < yMin) return false;
+    psAssert (xMin <= xMax, "x mismatch?");
+    psAssert (yMin <= yMax, "y mismatch?");
     for (int y = yMin; y <= yMax; y++) {
         for (int x = xMin; x <= xMax; x++) {
 …
                 ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel
             if (flux > *fluxStamp) {
                 *xStamp = x;
                 *yStamp = y;
+                *xStamp = x + 0.5;
+                *yStamp = y + 0.5;
                 *fluxStamp = flux;
                 found = true;
 …
 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
+                                                    int footprint, float spacing, float sysErr)
+                                                    int footprint, float spacing, float normFrac,
+                                                    float sysErr, float skyErr)
+{
     pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
 …
     list->y = NULL;
     list->flux = NULL;
+    list->window = NULL;
+    list->normFrac = normFrac;
+    list->normWindow = 0;
     list->footprint = footprint;
     list->sysErr = sysErr;
+    list->skyErr = skyErr;
     return list;
 …
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->footprint = in->footprint;
+    out->normWindow = in->normWindow;
     for (int i = 0; i < num; i++) {
 …
+        }
+        outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL;
+        outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL;
+        outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL;
+        outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL;
+        outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL;
+        outStamp->matrix = inStamp->matrix ? psImageCopy(NULL, inStamp->matrix, PS_TYPE_F64) : NULL;
+        outStamp->vector = inStamp->vector ? psVectorCopy(NULL, inStamp->vector, PS_TYPE_F64) : NULL;
         out->stamps->data[i] = outStamp;
 …
     stamp->convolutions2 = NULL;
+    stamp->matrix1 = NULL;
+    stamp->matrix2 = NULL;
+    stamp->matrixX = NULL;
+    stamp->vector1 = NULL;
+    stamp->vector2 = NULL;
+    stamp->matrix = NULL;
+    stamp->vector = NULL;
+    stamp->norm = NAN;
     return stamp;
 …
                                                 const psImage *image2, const psImage *subMask,
                                                 const psRegion *region, float thresh1, float thresh2,
                                                 int size, int footprint, float spacing, float sysErr,
                                                 pmSubtractionMode mode)
+                                                int size, int footprint, float spacing, float normFrac,
+                                                float sysErr, float skyErr, pmSubtractionMode mode)
+{
     if (!image1 && !image2) {
 …
     if (!stamps) {
+        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr);
+        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
+                                             normFrac, sysErr, skyErr);
+    }
 …
             numSearch++;
             int xStamp = 0, yStamp = 0; // Coordinates of stamp
+            float xStamp = 0.0, yStamp = 0.0; // Coordinates of stamp
             float fluxStamp = -INFINITY; // Flux of stamp
             bool goodStamp = false;     // Found a good stamp?
 …
                 // Take stamps off the top of the (sorted) list
                 for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
-                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
                     // Chop off the top of the list
                     xList->n = j;
 …
                     fluxList->n = j;
+#if 0
                     // Fish around a bit to see if we can find a pixel that isn't masked
+                    // This is not a good idea if we're using the window feature
                     psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n",
                             i, xCentre, yCentre);
                     // Search bounds
+                    int xCentre = xList->data.F32[j] - 0.5, yCentre = yList->data.F32[j] - 0.5;// Stamp centre
                     int search = footprint - size; // Search radius
                     int xMin = PS_MAX(border, xCentre - search);
 …
                     goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
                                             subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
+                    // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f (\n", xCentre, yCentre, xStamp, yStamp);
+#else
+                    // Only search the exact centre pixel
+                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
+                                            subMask, xList->data.F32[j], xList->data.F32[j],
+                                            yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
+#endif
+                    if (0 && goodStamp) {
+                        fprintf (stderr, "find: %6.1f < %6.1f < %6.1f, %6.1f < %6.1f %6.1f\n",
+                                 region->x0 + size, xStamp, region->x1 - size,
+                                 region->y0 + size, yStamp, region->y1 - size);
+                    }
+                }
             } else {
 …
                                                const psImage *image, const psImage *subMask,
                                                const psRegion *region, int size, int footprint,
+                                               float spacing, float sysErr, pmSubtractionMode mode)
+                                               float spacing, float normFrac, float sysErr, float skyErr,
+                                               pmSubtractionMode mode)
+{
 …
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
                                                                  region, footprint, spacing,
                                                                  sysErr); // Stamp list
+                                                                 normFrac, sysErr, skyErr); // Stamp list
     int numStamps = stamps->num;        // Number of stamps
 …
     for (int i = 0; i < numStars; i++) {
         float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp
         int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp
+        int xPix = xStamp - 0.5, yPix = yStamp - 0.5; // Pixel coordinate of stamp
         if (!checkStampRegion(xPix, yPix, region)) {
             // It's not in the big region
 …
             continue;
+        }
+        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
         bool found = false;
 …
+bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
+    int size = stamps->footprint; // Size of postage stamps
+    psFree (stamps->window);
+    stamps->window = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window->image, 0.0);
+    // generate normalizations for each stamp
+    psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    psVector *norm2 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (!stamp) continue;
+        if (!stamp->image1) continue;
+        if (!stamp->image2) continue;
+        float sum1 = 0.0;
+        float sum2 = 0.0;
+        for (int y = -size; y <= size; y++) {
+            for (int x = -size; x <= size; x++) {
+                sum1 += stamp->image1->kernel[y][x];
+                sum2 += stamp->image2->kernel[y][x];
+            }
+        }
+        norm1->data.F32[i] = sum1;
+        norm2->data.F32[i] = sum2;
+    }
+    // storage vector for flux data
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
+    psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
+    // generate the window pixels
+    double sum = 0.0;                   // Sum inside the window
+    float maxValue = 0.0;               // Maximum value, for normalisation
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            flux1->n = 0;
+            flux2->n = 0;
+            for (int i = 0; i < stamps->num; i++) {
+                pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+                if (!stamp) continue;
+                if (!stamp->image1) continue;
+                if (!stamp->image2) continue;
+                psVectorAppend(flux1, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
+                psVectorAppend(flux2, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
+            }
+            psStatsInit (stats);
+            if (!psVectorStats (stats, flux1, NULL, NULL, 0)) {
+                psAbort ("failed to generate stats");
+            }
+            float f1 = stats->robustMedian;
+            psStatsInit (stats);
+            if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
+                psAbort ("failed to generate stats");
+            }
+            float f2 = stats->robustMedian;
+            stamps->window->kernel[y][x] = f1 + f2;
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
+                sum += stamps->window->kernel[y][x];
+            }
+            maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
+        }
+    }
+    psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
+            sum, (1.0 - stamps->normFrac) * sum);
+    bool done = false;
+    for (int radius = 1; radius <= size && !done; radius++) {
+        double within = 0.0;
+        for (int y = -radius; y <= radius; y++) {
+            for (int x = -radius; x <= radius; x++) {
+                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
+                    within += stamps->window->kernel[y][x];
+                }
+            }
+        }
+        psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
+        if (within > (1.0 - stamps->normFrac) * sum) {
+            stamps->normWindow = radius;
+            done = true;
+        }
+    }
+    psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
+    if (stamps->normWindow == 0 || stamps->normWindow >= size) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to determine normalisation window size.");
+        return false;
+    }
+    // re-normalize so chisquare values are sensible
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            stamps->window->kernel[y][x] /= maxValue;
+        }
+    }
+#if 0
+    {
+        psFits *fits = psFitsOpen ("window.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
+        psFitsClose (fits);
+    }
+#endif
+    psFree (stats);
+    psFree (flux1);
+    psFree(flux2);
+    psFree (norm1);
+    psFree (norm2);
+    return true;
+}
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
                                 psImage *variance, int kernelSize)
+                                psImage *variance, int kernelSize, psRegion bounds)
+{
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
 …
+    }
-    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
     int size = kernelSize + stamps->footprint; // Size of postage stamps
 …
+        }
+        if (isnan(stamp->xNorm)) {
+            stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
+        }
+        if (isnan(stamp->yNorm)) {
+            stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
+        }
+        int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
+        if (x < size || x > numCols - size || y < size || y > numRows - size) {
+            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
+        p_pmSubtractionPolynomialNormCoords(&stamp->xNorm, &stamp->yNorm, stamp->x, stamp->y,
+                                            bounds.x0, bounds.x1, bounds.y0, bounds.y1);
+        int x = stamp->x - 0.5, y = stamp->y - 0.5; // Stamp coordinates
+        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d (size: %d)\n", stamp->x, stamp->y, x, y, size);
+        if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) {
+            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
             return false;
+        }
 …
             psImage *varSub = psImageSubset(variance, region); // Subimage with stamp
             psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp
+            if (isfinite(stamps->sysErr) && stamps->sysErr > 0) {
+            if ((isfinite(stamps->skyErr) && (stamps->skyErr > 0)) ||
+                (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) {
                 float sysErr = 0.25 * PS_SQR(stamps->sysErr); // Systematic error
+                float skyErr = stamps->skyErr;
                 psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images
                 for (int y = -size; y <= size; y++) {
                     for (int x = -size; x <= size; x++) {
                         float additional = image1->kernel[y][x] + image2->kernel[y][x];
                         weight->kernel[y][x] = 1.0 / (var->kernel[y][x] + sysErr * PS_SQR(additional));
+                        weight->kernel[y][x] = 1.0 / (skyErr + var->kernel[y][x] + sysErr * PS_SQR(additional));
+                    }
+                }
 …
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
                                                           const psImage *subMask, const psRegion *region,
+                                                          int size, int footprint, float spacing, float sysErr,
+                                                          int size, int footprint, float spacing,
+                                                          float normFrac, float sysErr, float skyErr,
                                                           pmSubtractionMode mode)
+{
 …
             y->data.F32[numOut] = source->peak->yf;
+        }
+        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
         numOut++;
+    }
 …
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
+                                                            footprint, spacing, sysErr, mode); // Stamps
+                                                            footprint, spacing, normFrac,
+                                                            sysErr, skyErr, mode); // Stamps
     psFree(x);
     psFree(y);
 …
 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
                                                        const psImage *subMask, const psRegion *region,
                                                        int size, int footprint, float spacing, float sysErr,
                                                        pmSubtractionMode mode)
+                                                       int size, int footprint, float spacing, float normFrac,
+                                                       float sysErr, float skyErr, pmSubtractionMode mode)
+{
     PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
 …
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
                                                             spacing, sysErr, mode);
+                                                            spacing, normFrac, sysErr, skyErr, mode);
     psFree(data);
 …
+}
+bool pmSubtractionStampsResetStatus (pmSubtractionStampList *stamps) {
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (!stamp) continue;
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
+    }
+    return true;
+}

trunk/psModules/src/imcombine/pmSubtractionStamps.h

-              r26035
+              r26893
     psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
     int footprint;                      ///< Half-size of stamps
+    psKernel *window;                   ///< window function generated from ensemble of stamps
+    float normFrac;                     ///< Fraction of flux in window for normalisation window
+    int normWindow;                     ///< Size of window for measuring normalisation
     float sysErr;                       ///< Systematic error
+    float skyErr;                       ///< increase effective readnoise
 } pmSubtractionStampList;
 /// Allocate a list of stamps
+pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image
+                                                    int numRows, // Number of rows in image
+                                                    const psRegion *region, // Region for stamps, or NULL
+                                                    int footprint, // Half-size of stamps
+                                                    float spacing, // Rough average spacing between stamps
+                                                    float sysErr  // Relative systematic error or NAN
+pmSubtractionStampList *pmSubtractionStampListAlloc(
+    int numCols, // Number of columns in image
+    int numRows, // Number of rows in image
+    const psRegion *region, // Region for stamps, or NULL
+    int footprint, // Half-size of stamps
+    float spacing, // Rough average spacing between stamps
+    float normFrac, // Fraction of flux in window for normalisation window
+    float sysErr,  // Relative systematic error or NAN
+    float skyErr  // Relative systematic error or NAN
     );
 …
     psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
     psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
     psImage *matrix1, *matrix2;         ///< Least-squares matrices for each image, or NULL
     psImage *matrixX;                   ///< Cross-matrix (for mode DUAL), or NULL
     psVector *vector1, *vector2;        ///< Least-squares vectors for each image, or NULL
+    psImage *matrix;                    ///< Least-squares matrix, or NULL
+    psVector *vector;                   ///< Least-squares vector, or NULL
+    double norm;                        ///< Normalisation difference
     pmSubtractionStampStatus status;    ///< Status of stamp
 } pmSubtractionStamp;
 …
 /// Find stamps on an image
+pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL
+                                                const psImage *image1, ///< Image for which to find stamps
+                                                const psImage *image2, ///< Image for which to find stamps
+                                                const psImage *mask, ///< Mask, or NULL
+                                                const psRegion *region, ///< Region to search, or NULL
+                                                float thresh1, ///< Threshold for stamps in image 1
+                                                float thresh2, ///< Threshold for stamps in image 2
+                                                int size, ///< Kernel half-size
+                                                int footprint, ///< Half-size for stamps
+                                                float spacing, ///< Rough spacing for stamps
+                                                float sysErr,  ///< Relative systematic error in images
+                                                pmSubtractionMode mode ///< Mode for subtraction
+pmSubtractionStampList *pmSubtractionStampsFind(
+    pmSubtractionStampList *stamps, ///< Output stamps, or NULL
+    const psImage *image1, ///< Image for which to find stamps
+    const psImage *image2, ///< Image for which to find stamps
+    const psImage *mask, ///< Mask, or NULL
+    const psRegion *region, ///< Region to search, or NULL
+    float thresh1, ///< Threshold for stamps in image 1
+    float thresh2, ///< Threshold for stamps in image 2
+    int size, ///< Kernel half-size
+    int footprint, ///< Half-size for stamps
+    float spacing, ///< Rough spacing for stamps
+    float normFrac, // Fraction of flux in window for normalisation window
+    float sysErr,  ///< Relative systematic error in images
+    float skyErr,  ///< Relative systematic error in images
+    pmSubtractionMode mode ///< Mode for subtraction
     );
 /// Set stamps based on a list of x,y
+pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp
+                                               const psVector *y, ///< y coordinates for each stamp
+                                               const psImage *image, ///< Image for flux of stamp
+                                               const psImage *mask, ///< Mask, or NULL
+                                               const psRegion *region, ///< Region to search, or NULL
+                                               int size, ///< Kernel half-size
+                                               int footprint, ///< Half-size for stamps
+                                               float spacing, ///< Rough spacing for stamps
+                                               float sysErr,  ///< Systematic error in images
+                                               pmSubtractionMode mode ///< Mode for subtraction
+pmSubtractionStampList *pmSubtractionStampsSet(
+    const psVector *x, ///< x coordinates for each stamp
+    const psVector *y, ///< y coordinates for each stamp
+    const psImage *image, ///< Image for flux of stamp
+    const psImage *mask, ///< Mask, or NULL
+    const psRegion *region, ///< Region to search, or NULL
+    int size, ///< Kernel half-size
+    int footprint, ///< Half-size for stamps
+    float spacing, ///< Rough spacing for stamps
+    float normFrac, // Fraction of flux in window for normalisation window
+    float sysErr,  ///< Systematic error in images
+    float skyErr,  ///< Systematic error in images
+    pmSubtractionMode mode ///< Mode for subtraction
     );
 …
     int footprint,                      ///< Half-size for stamps
     float spacing,                      ///< Rough spacing for stamps
+    float normFrac, // Fraction of flux in window for normalisation window
     float sysErr,                       ///< Systematic error in images
+    float skyErr,                       ///< Systematic error in images
     pmSubtractionMode mode              ///< Mode for subtraction
     );
 …
     int footprint,                      ///< Half-size for stamps
     float spacing,                      ///< Rough spacing for stamps
+    float normFrac, // Fraction of flux in window for normalisation window
     float sysErr,                       ///< Systematic error in images
+    float skyErr,                       ///< Systematic error in images
     pmSubtractionMode mode              ///< Mode for subtraction
+    );
+/// Calculate the window and normalisation window from the stamps
+bool pmSubtractionStampsGetWindow(
+    pmSubtractionStampList *stamps,     ///< List of stamps
+    int kernelSize                      ///< Half-size of kernel
     );
 …
                                 psImage *image2, ///< Input image (or NULL)
                                 psImage *variance, ///< Variance map
+                                int kernelSize ///< Kernel half-size
+                                int kernelSize, ///< Kernel half-size
+                                psRegion bounds ///< Bounds of validity
     );
 …
     );
+bool pmSubtractionStampsResetStatus (pmSubtractionStampList *stamps);
 #endif

trunk/psModules/src/imcombine/pmSubtractionThreads.c

r21363	r26893
69	69
70	70	{
71		psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 3);
	71	psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 4);
72	72	task->function = &pmSubtractionCalculateEquationThread;
73	73	psThreadTaskAdd(task);

trunk/psModules/src/imcombine/pmSubtractionVisual.c

-              r23487
+              r26893
 #include "pmKapaPlots.h"
+#include "pmSubtraction.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionEquation.h"
 #include "pmVisual.h"
 …
 static bool plotLeastSquares     = true;
 static bool plotImage            = true;
 // variables to store plotting window indices
 static int kapa  = -1;
+static int kapa1 = -1;
 static int kapa2 = -1;
+static int kapa3 = -1;
 /** function prototypes*/
 …
 bool pmSubtractionVisualClose(void)
+{
+    if(kapa != -1)
+        KiiClose(kapa);
+    if(kapa2 != -1)
+        KiiClose(kapa2);
+    if(kapa1 != -1) KiiClose(kapa1);
+    if(kapa2 != -1) KiiClose(kapa2);
     return true;
+}
 …
 bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {
     if (!pmVisualIsVisual() || !plotConvKernels) return true;
     if (!pmVisualInitWindow(&kapa, "ppSub:Images")) {
         return false;
+    }
     pmVisualScaleImage(kapa, convKernels, "Convolution_Kernels", 0, true);
+    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) {
+        return false;
+    }
+    pmVisualScaleImage(kapa1, convKernels, "Convolution_Kernels", 0, true);
     pmVisualAskUser(&plotConvKernels);
     return true;
 …
 bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
     if (!pmVisualIsVisual() || !plotStamps) return true;
     if (!pmVisualInitWindow (&kapa, "ppSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "ppSub:Images")) {
         return false;
+    }
 …
         stampNum++;
+    }
     pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true);
+    pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);
     pmVisualAskUser(&plotStamps);
 …
     if (!pmVisualIsVisual() || !plotLeastSquares) return true;
     if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
         return false;
+    }
 …
         if (stamp == NULL) continue;
+        psImage *im = stamp->matrix1;
+        if (im == NULL) im = stamp->matrix2;
+        if (im == NULL) im = stamp->matrixX;
+        psImage *im = stamp->matrix;
         if (im == NULL) continue;
 …
         if (stamp == NULL) continue;
+        psImage *im = stamp->matrix1;
+        if (im == NULL) im = stamp->matrix2;
+        if (im == NULL) im = stamp->matrixX;
+        psImage *im = stamp->matrix;
         if (im == NULL) continue;
 …
     psImage *canvas32 = pmVisualImageToFloat(canvas);
     pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true);
+    pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
     pmVisualAskUser(&plotLeastSquares);
 …
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
     if (!pmVisualIsVisual() || !plotImage) return true;
+    if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
+        return false;
+    }
+    pmVisualScaleImage(kapa, image, "Image", 0, true);
+    pmVisualScaleImage(kapa, ref, "Reference", 1, true);
+    pmVisualScaleImage(kapa, sub, "Subtraction", 2, true);
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
+        return false;
+    }
+    pmVisualScaleImage(kapa1, image, "Image", 0, true);
+    pmVisualScaleImage(kapa1, ref, "Reference", 1, true);
+    pmVisualScaleImage(kapa1, sub, "Subtraction", 2, true);
+    pmVisualAskUser(&plotImage);
+    return true;
+}
+bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) {
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {
+        return false;
+    }
+    // get the kernel sizes
+    int footprint = kernels->size;
+    // output image is a grid of NXsub by NYsub sub-images
+    int NXsub = sqrt(kernels->num);
+    int NYsub = kernels->num / NXsub;
+    if (kernels->num % NXsub) NYsub++;
+    int NXpix = NXsub * (2*footprint + 1 + 3);
+    int NYpix = NYsub * (2*footprint + 1 + 3);
+    psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
+    psImageInit (output, 0.0);
+    for (int i = 0; i < kernels->num; i++) {
+        pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i];
+        psKernel *kernel = preCalc->kernel;
+        int xSub = i % NXsub;
+        int ySub = i / NXsub;
+        int xPix = xSub * (2*footprint + 1 + 3) + footprint;
+        int yPix = ySub * (2*footprint + 1 + 3) + footprint;
+        double sum = 0.0;
+        for (int y = -footprint; y <= footprint; y++) {
+            for (int x = -footprint; x <= footprint; x++) {
+                output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                sum += kernel->kernel[y][x];
+            }
+        }
+        fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+    }
+    pmVisualScaleImage(kapa1, output, "Image", 0, true);
+    pmVisualAskUser(&plotImage);
+    return true;
+}
+bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps) {
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) {
+        return false;
+    }
+    // get the kernel sizes
+    int footprint = stamps->footprint;
+    // choose the brightest stamp
+    pmSubtractionStamp *maxStamp = NULL;
+    float maxFlux = NAN;
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i];
+        if (!isfinite(stamp->flux)) continue;
+        if (!stamp->convolutions1 && !stamp->convolutions2) continue;
+        if (!maxStamp) {
+            maxFlux = stamp->flux;
+            maxStamp = stamp;
+            continue;
+        }
+        if (stamp->flux > maxFlux) {
+            maxFlux = stamp->flux;
+            maxStamp = stamp;
+        }
+    }
+    if (!isfinite(maxStamp->flux)) {
+        fprintf (stderr, "no valid stamps?\n");
+    }
+    int nKernels = 0;
+    if (maxStamp->convolutions1) {
+        // output image is a grid of NXsub by NYsub sub-images
+        nKernels = maxStamp->convolutions1->n;
+        int NXsub = sqrt(nKernels);
+        int NYsub = nKernels / NXsub;
+        if (nKernels % NXsub) NYsub++;
+        int NXpix = NXsub * (2*footprint + 1 + 3);
+        int NYpix = NYsub * (2*footprint + 1 + 3);
+        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
+        psImageInit (output, 0.0);
+        for (int i = 0; i < nKernels; i++) {
+            psKernel *kernel = maxStamp->convolutions1->data[i];
+            int xSub = i % NXsub;
+            int ySub = i / NXsub;
+            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
+            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
+            double sum = 0.0;
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                    sum += kernel->kernel[y][x];
+                }
+            }
+            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+        }
+        pmVisualScaleImage(kapa2, output, "Image", 0, true);
+    }
+    if (maxStamp->convolutions2) {
+        // output image is a grid of NXsub by NYsub sub-images
+        nKernels = maxStamp->convolutions2->n;
+        int NXsub = sqrt(nKernels);
+        int NYsub = nKernels / NXsub;
+        if (nKernels % NXsub) NYsub++;
+        int NXpix = NXsub * (2*footprint + 1 + 3);
+        int NYpix = NYsub * (2*footprint + 1 + 3);
+        psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32);
+        psImageInit (output, 0.0);
+        for (int i = 0; i < nKernels; i++) {
+            psKernel *kernel = maxStamp->convolutions2->data[i];
+            int xSub = i % NXsub;
+            int ySub = i / NXsub;
+            int xPix = xSub * (2*footprint + 1 + 3) + footprint;
+            int yPix = ySub * (2*footprint + 1 + 3) + footprint;
+            double sum = 0.0;
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                    sum += kernel->kernel[y][x];
+                }
+            }
+            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+        }
+        pmVisualScaleImage(kapa2, output, "Image", 1, true);
+    }
     pmVisualAskUser(&plotImage);
     return true;
 …
         overlay[Noverlay].type = KII_OVERLAY_BOX;
+        if ((stamp->x < 1.0) && (stamp->y < 1.0)) {
+            // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
+            continue;
+        }
+        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
+            // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
+            continue;
+        }
         overlay[Noverlay].x = stamp->x;
         overlay[Noverlay].y = stamp->y;
 …
+}
+static int footprint = 0;
+static int NX = 0;
+static int NY = 0;
+static psImage *sourceImage      = NULL;
+static psImage *targetImage      = NULL;
+static psImage *residualImage    = NULL;
+static psImage *fresidualImage   = NULL;
+static psImage *differenceImage  = NULL;
+static psImage *convolutionImage = NULL;
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
+    if (!pmVisualIsVisual()) return true;
+    // generate 4 storage images large enough to hold the N stamps:
+    footprint = stamps->footprint;
+    float NXf = sqrt(stamps->num);
+    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
+    float NYf = stamps->num / NX;
+    NY = (int) NYf == NY ? NYf : NYf + 1.0;
+    int NXpix = (2*footprint + 1) * NX;
+    NXpix += (NX > 1) ? 3 * NX : 0;
+    int NYpix = (2*footprint + 1) * NY;
+    NYpix += (NY > 1) ? 3 * NY : 0;
+    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    psImageInit (sourceImage,      0.0);
+    psImageInit (targetImage,      0.0);
+    psImageInit (residualImage,    0.0);
+    psImageInit (fresidualImage,   0.0);
+    psImageInit (differenceImage,  0.0);
+    psImageInit (convolutionImage, 0.0);
+    return true;
+}
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
+    if (!pmVisualIsVisual()) return true;
+    double sum;
+    int NXoff = index % NX;
+    int NYoff = index / NX;
+    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
+    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
+    // insert the (target) kernel into the (target) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
+            sum += targetImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    // insert the (source) kernel into the (source) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
+            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    // insert the (convolution) kernel into the (convolution) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x];
+            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    // insert the (difference) kernel into the (difference) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
+            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    // insert the (residual) kernel into the (residual) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
+            sum += residualImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    // insert the (fresidual) kernel into the (fresidual) image:
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
+        }
+    }
+    return true;
+}
+bool pmSubtractionVisualShowFit(double norm) {
+    // for testing, dump the residual image and exit
+    if (0) {
+        psMetadata *header = psMetadataAlloc();
+        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
+        psFits *fits = psFitsOpen("resid.fits", "w");
+        psFitsWriteImage(fits, header, residualImage, 0, NULL);
+        psFitsClose(fits);
+        // exit (0);
+    }
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
+    if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
+    KiiEraseOverlay (kapa1, "red");
+    KiiEraseOverlay (kapa2, "red");
+    pmVisualScaleImage(kapa1, targetImage, "Target Stamps", 0, true);
+    pmVisualScaleImage(kapa1, sourceImage, "Source Stamps", 1, true);
+    pmVisualScaleImage(kapa1, convolutionImage, "Convolution Stamps", 2, true);
+    KiiCenter (kapa1, 0.5*targetImage->numCols, 0.5*targetImage->numRows, 1);
+    pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true);
+    pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);
+    if (1) {
+        KiiImage image;
+        KapaImageData data;
+        Coords coords;
+        strcpy (coords.ctype, "RA---TAN");
+        image.data2d = residualImage->data.F32;
+        image.Nx = residualImage->numCols;
+        image.Ny = residualImage->numRows;
+        strcpy (data.name, "Residual Stamps");
+        strcpy (data.file, "Residual Stamps");
+        data.zero  = -300.0;
+        data.range = +600.0;
+        data.logflux = 0;
+        KiiSetChannel (kapa2, 1);
+        KiiNewPicture2D (kapa2, &image, &data, &coords);
+    } else {
+        pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true);
+    }
+    KiiCenter (kapa2, 0.5*residualImage->numCols, 0.5*residualImage->numRows, 1);
+    pmVisualAskUser(NULL);
+    psFree(targetImage);
+    psFree(sourceImage);
+    psFree(convolutionImage);
+    psFree(differenceImage);
+    psFree(residualImage);
+    psFree(fresidualImage);
+    targetImage = NULL;
+    sourceImage = NULL;
+    convolutionImage = NULL;
+    differenceImage = NULL;
+    residualImage = NULL;
+    fresidualImage = NULL;
+    return true;
+}
+bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels) {
+    Graphdata graphdata;
+    if (!pmVisualIsVisual()) return true;
+    if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false;
+    KapaClearSections (kapa3);
+    KapaInitGraph (&graphdata);
+    psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
+    psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
+    graphdata.xmin = -1.0;
+    graphdata.xmax = kernels->num + 1.0;
+    graphdata.ymin = +32.0;
+    graphdata.ymax = -32.0;
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0);
+    // construct the plot vectors
+    for (int i = 0; i < kernels->num; i++) {
+        x->data.F32[i] = i;
+        y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
+        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
+        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
+    }
+    x->n = y->n = kernels->num;
+    float range;
+    range = graphdata.xmax - graphdata.xmin;
+    graphdata.xmax += 0.05*range;
+    graphdata.xmin -= 0.05*range;
+    range = graphdata.ymax - graphdata.ymin;
+    graphdata.ymax += 0.05*range;
+    graphdata.ymin -= 0.05*range;
+    KapaSetLimits (kapa3, &graphdata);
+    KapaSetFont (kapa3, "helvetica", 14);
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "kernel number", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "coeff", KAPA_LABEL_YM);
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot   (kapa3, x->n, &graphdata);
+    KapaPlotVector (kapa3, x->n, x->data.F32, "x");
+    KapaPlotVector (kapa3, x->n, y->data.F32, "y");
+    psFree (x);
+    psFree (y);
+    pmVisualAskUser(NULL);
+    return true;
+}
 #else
 bool pmSubtractionVisualClose(void) {return true;}
 …
 bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps) {return true;}
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {return true;}
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {return true;}
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {return true;}
+bool pmSubtractionVisualShowFit() {return true;}
+bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
 #endif

trunk/psModules/src/imcombine/pmSubtractionVisual.h

-              r23487
+              r26893
 bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps);
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub);
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps);
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index);
+bool pmSubtractionVisualShowFit(double norm);
+bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
+bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels);
+bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps);
 #endif

Context Navigation

Legend:

trunk/psModules/src/imcombine/Makefile.am

trunk/psModules/src/imcombine/pmImageCombine.c

trunk/psModules/src/imcombine/pmPSFEnvelope.c

trunk/psModules/src/imcombine/pmPSFEnvelope.h

trunk/psModules/src/imcombine/pmStackReject.c

trunk/psModules/src/imcombine/pmSubtraction.c

trunk/psModules/src/imcombine/pmSubtraction.h

trunk/psModules/src/imcombine/pmSubtractionAnalysis.c

trunk/psModules/src/imcombine/pmSubtractionAnalysis.h

trunk/psModules/src/imcombine/pmSubtractionEquation.c

trunk/psModules/src/imcombine/pmSubtractionEquation.h

trunk/psModules/src/imcombine/pmSubtractionIO.c

trunk/psModules/src/imcombine/pmSubtractionKernels.c

trunk/psModules/src/imcombine/pmSubtractionKernels.h

trunk/psModules/src/imcombine/pmSubtractionMask.c

trunk/psModules/src/imcombine/pmSubtractionMask.h

trunk/psModules/src/imcombine/pmSubtractionMatch.c

trunk/psModules/src/imcombine/pmSubtractionMatch.h

trunk/psModules/src/imcombine/pmSubtractionParams.c

trunk/psModules/src/imcombine/pmSubtractionParams.h

trunk/psModules/src/imcombine/pmSubtractionStamps.c

trunk/psModules/src/imcombine/pmSubtractionStamps.h

trunk/psModules/src/imcombine/pmSubtractionThreads.c

trunk/psModules/src/imcombine/pmSubtractionVisual.c

trunk/psModules/src/imcombine/pmSubtractionVisual.h

Download in other formats: