Changeset 26574

Timestamp:

Jan 12, 2010, 3:55:01 PM (16 years ago)

Author:

eugene

Message:

add ISIS_RADIAL (ISIS + radial Hermitian); modify pmSubtractionKernelPreCalcNormalize to force zero null on odd terms

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (modified) (28 diffs)

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

@@ -90 +90 @@
 
     for (int i = 0, x = -size; x <= size; i++, x++) {
-        float xf = x / sigma;
-        float z = -0.25*xf*xf;
+        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);
+        }
     }
 
@@ -132 +156 @@
             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");
-
+            psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
         }
@@ -147 +170 @@
 }
 
-bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm, bool AlardLuptonStyle) {
+bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm, bool AlardLuptonStyle, bool forceZeroNull) {
 
     // Calculate moments
     double moment = 0.0;    // Moment, for penalty
     for (int v = -size; v <= size; v++) {
-        for (int u = -size; u <= size; u++) {
-            double value = preCalc->kernel->kernel[v][u];
-            moment += PS_SQR(value) * PS_SQR((PS_SQR(u) + PS_SQR(v)));
-        }
+        for (int u = -size; u <= size; u++) {
+            double value = preCalc->kernel->kernel[v][u];
+            moment += PS_SQR(value) * PS_SQR((PS_SQR(u) + PS_SQR(v)));
+        }
     }
 
@@ -169 +192 @@
 
     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);
-        }
+        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);
+        }
     }
 #if 0
@@ -181 +204 @@
     // only even terms have non-zero sums
     if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) {
-        moment /= PS_SQR(sum);
+        moment /= PS_SQR(sum);
+    } else {
+        // XXX keep this?
+        moment = 0.0;
     }
 
@@ -189 +215 @@
 
     if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) {
-        zeroNull = true;
+        zeroNull = true;
     }
     if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) {
-        zeroNull = true;
-    }
-    if ((uOrder % 2) || (vOrder % 2)) {
-        // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max);
-        scale2D = 1.0 / max;
-        scale1D = sqrt(scale2D);
-    }
-
-    psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
-    psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+        zeroNull = true;
+    }
+    if (forceZeroNull) {
+        zeroNull = true;
+    }
+    if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
+        // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max);
+        scale2D = 1.0 / max;
+        scale1D = sqrt(scale2D);
+    }
+
+    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));
     if (zeroNull) {
-        preCalc->kernel->kernel[0][0] -= 1.0;
+        preCalc->kernel->kernel[0][0] -= 1.0;
     }
 
@@ -212 +245 @@
     max = FLT_MIN;
     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);
-        }
+        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);
@@ -225 +258 @@
     kernels->v->data.S32[index] = vOrder;
     if (kernels->preCalc->data[index]) {
-        psFree(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,
-            fwhm, uOrder, vOrder, fabsf(moment));
+    psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
+    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, fabsf(moment));
 
     return true;
@@ -252 +284 @@
     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]);
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
     }
 
@@ -273 +305 @@
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS 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);
+                // 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, 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, mode); // The 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);
 
@@ -282 +369 @@
             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);
+                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, pmSubtractionMode mode)
+                                               const psVector *fwhmsIN, const psVector *ordersIN,
+                                               float penalty, pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
@@ -306 +397 @@
     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]);
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
     }
 
@@ -334 +425 @@
             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);
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
             }
         }
@@ -343 +434 @@
 
 pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
-                                                     const psVector *fwhmsIN, const psVector *ordersIN,
-                                                     float penalty, pmSubtractionMode mode)
+                                                     const psVector *fwhmsIN, const psVector *ordersIN,
+                                                     float penalty, pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
@@ -358 +449 @@
     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]);
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
     }
 
@@ -397 +488 @@
                 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;
+                // 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);
-
-                // XXXX test demo that deconvolved kernel is valid
+                // 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);
-                }
+                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);
+                }
 # endif
             }
@@ -421 +512 @@
     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];
-                }
-            }
-            dot->data.F32[j][i] = sum;
-            dot->data.F32[i][j] = sum;
-        }
+        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];
+                }
+            }
+            dot->data.F32[j][i] = sum;
+            dot->data.F32[i][j] = sum;
+        }
     }
     pmSubtractionVisualShowSubtraction (dot, NULL, NULL);
@@ -460 +551 @@
     kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32);
     kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
+    kernels->uStop = NULL;
+    kernels->vStop = NULL;
     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;
@@ -474 +565 @@
     kernels->solution1 = NULL;
     kernels->solution2 = NULL;
+    kernels->mean = NAN;
+    kernels->rms = NAN;
+    kernels->numStamps = 0;
+
+    kernels->fSigResMean  = NAN;
+    kernels->fSigResStdev = NAN;
+    kernels->fMaxResMean  = NAN;
+    kernels->fMaxResStdev = NAN;
+    kernels->fMinResMean  = NAN;
+    kernels->fMinResStdev = NAN;
 
     return kernels;
@@ -486 +587 @@
     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;
+        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_DECONV_HERM:
-        preCalc->xKernel = pmSubtractionKernelHERM(sigma, uOrder, size);
-        preCalc->yKernel = pmSubtractionKernelHERM(sigma, vOrder, size);
-        preCalc->uCoords = NULL;
-        preCalc->vCoords = NULL;
-        preCalc->poly    = NULL;
-        break;
+        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;
+        // 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");
+        psAbort("programming error: invalid type for PreCalc kernel");
     }
 
@@ -524 +626 @@
     // 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
-        }
-    }
-
+        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;
 }
@@ -753 +855 @@
 }
 
-// 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,
@@ -856 +958 @@
             for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) {
 
-                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
+                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0);
                 double moment = 0.0;    // Moment, for penalty
 
@@ -862 +964 @@
                     // Central pixel is easy
                     preCalc->uCoords->data.S32[0] = 0;
-                    preCalc->vCoords->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;
+                    preCalc->vCoords->n = 1;
+                    preCalc->poly->n = 1;
                     radiusLast = 0;
                     moment = 0.0;
@@ -923 +1025 @@
                 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");
-                }
+                if (!isfinite(kernels->penalties->data.F32[index])) {
+                    psAbort ("invalid penalty");
+                }
 
                 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,
@@ -946 +1048 @@
       case PM_SUBTRACTION_KERNEL_ISIS:
         return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
+      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
+        return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, mode);
       case PM_SUBTRACTION_KERNEL_HERM:
         return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode);
@@ -1012 +1116 @@
     float penalty = 0.0;                // Penalty for wideness
 
-    // ISIS, HERM, and DECONV_HERM have the same description layout
-    if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM") || !strcmp (description, "DECONV_HERM")) {
-        // XXX Support for GUNK (not yet supported)
-        if (strstr(description, "+POIS")) {
-            type = PM_SUBTRACTION_KERNEL_GUNK;
-            psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description);
-        }
-
-        type = pmSubtractionKernelsTypeFromString (description);
-        psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must  be ISIS, HERM or DECONV_HERM");
-
-        char *ptr = NULL;
-        switch (type) {
-          case PM_SUBTRACTION_KERNEL_ISIS:
-          case PM_SUBTRACTION_KERNEL_HERM:
-            ptr = (char*) description + 5;    // Eat "ISIS(" or "HERM("
-            break;
-          case PM_SUBTRACTION_KERNEL_DECONV_HERM:
-            ptr = (char*) description + 12;    // Eat "DECONV_HERM("
-            break;
-          default:
-            psAbort("programming error: invalid kernel type");
-        }
-        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);
-
-        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
-    }
-
-    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, '(');
+    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);
+
+        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
+
+      case PM_SUBTRACTION_KERNEL_RINGS:
         PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
         PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt);
@@ -1067 +1159 @@
         PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
         PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
-        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
-    }
-
-    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, mode);
-}
-
-
+        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
+      default:
+        psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
+    } 
+    return NULL;
+}
+
+
+// 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, "HERM") == 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 (strcasecmp(type, "DECONV_HERM") == 0) {
+    if (strncasecmp(type, "DECONV_HERM", nameLength) == 0) {
         return PM_SUBTRACTION_KERNEL_DECONV_HERM;
     }
-    if (strcasecmp(type, "SPAM") == 0) {
+    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;
     }