Changeset 27838 for branches/tap_branches/psModules/src/imcombine/pmSubtractionKernels.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionKernels.c (modified) (30 diffs)

branches/tap_branches/psModules/src/imcombine/pmSubtractionKernels.c

@@ -10 +10 @@
 #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
@@ -27 +29 @@
     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);
 }
 
@@ -45 +60 @@
 
 // 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
@@ -57 +72 @@
     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);
+        }
     }
 
@@ -81 +142 @@
     kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     kernels->inner = start;
+    kernels->num += numNew;
 
     // Generate a set of kernels for each (u,v)
@@ -94 +156 @@
             kernels->v->data.S32[index] = v;
             kernels->preCalc->data[index] = NULL;
-            kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v));
-
+            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);
         }
     }
 
+    kernels->widths->n = start + numNew;
+    kernels->u->n = start + numNew;
+    kernels->v->n = start + numNew;
+    kernels->preCalc->n = start + numNew;
+    kernels->penalties->n = start + numNew;
+
     return true;
 }
 
+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
@@ -126 +306 @@
 
     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);
 
@@ -134 +314 @@
 
     // 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
@@ -140 +319 @@
         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(2); // 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
-
-                // 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
-                        moment += value * (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);
+            }
+        }
+    }
+
+    psFree(orders);
+    psFree(fwhms);
+
+    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));
-                    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];
                 }
-
-                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;
@@ -190 +552 @@
 
 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
-                                                int size, int spatialOrder, float penalty,
+                                                int size, int spatialOrder, float penalty, psRegion bounds,
                                                 pmSubtractionMode mode)
 {
@@ -202 +564 @@
     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;
@@ -212 +578 @@
     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)
 {
@@ -228 +658 @@
     int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, size,
-                                                              spatialOrder, penalty, mode); // The kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
+                                                              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",
@@ -244 +674 @@
 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;
@@ -256 +686 @@
 
 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);
@@ -277 +707 @@
 
     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,
@@ -348 +778 @@
 
 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
-                                                pmSubtractionMode mode)
+                                                psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -375 +805 @@
 
     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);
@@ -441 +871 @@
 }
 
-// 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);
@@ -456 +886 @@
     PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
 
-    // XXX GUNK doesn't seem to work --- doesn't add the POIS components, or at least, they're not noticed
-
     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=");
     psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
@@ -474 +903 @@
 // 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);
@@ -505 +934 @@
 
     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,
@@ -545 +974 @@
             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;
@@ -573 +1002 @@
                                 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(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++;
                                 }
@@ -591 +1020 @@
                     // 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,
@@ -624 +1056 @@
 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);
@@ -675 +1113 @@
 
 pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder,
-                                                          pmSubtractionMode mode)
+                                                          psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_STRING_NON_EMPTY(description, NULL);
@@ -694 +1132 @@
     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);
@@ -731 +1173 @@
         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;
     }
@@ -789 +1255 @@
     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;