Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26573) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26574) @@ -90,7 +90,31 @@ 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,6 +156,5 @@ 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,13 +170,13 @@ } -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,9 +192,9 @@ 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,5 +204,8 @@ // 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,20 +215,27 @@ 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,9 +245,9 @@ 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,11 +258,10 @@ 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,7 +284,7 @@ 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,4 +305,59 @@ 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(¶ms, "(%.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,15 +369,19 @@ 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,7 +397,7 @@ 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,5 +425,5 @@ 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,6 +434,6 @@ 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,7 +449,7 @@ 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,20 +488,20 @@ 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,17 +512,17 @@ 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,9 +551,9 @@ 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,4 +565,14 @@ 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,36 +587,37 @@ 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,9 +626,9 @@ // 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,5 +855,5 @@ } -// 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,5 +958,5 @@ 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,9 +964,9 @@ // 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,7 +1025,7 @@ 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,4 +1048,6 @@ 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,54 +1116,42 @@ 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,38 +1159,54 @@ 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; }