Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26586) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26587) @@ -90,6 +90,6 @@ 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); } @@ -100,6 +100,6 @@ // 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 + int order, // Polynomial order + int size // Kernel half-size ) { @@ -112,9 +112,9 @@ // 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); - } + 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); + } } @@ -156,5 +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); } @@ -170,15 +170,8 @@ } -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))); - } - } - +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 @@ -187,55 +180,64 @@ // 4) for deconvolved hermitians, subtract 1 from the 0,0 pixel for the 0,0 function(s) - double sum = 0.0; // Sum of kernel component - double min = FLT_MAX; - double max = FLT_MIN; - + // 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++) { - 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++) { + 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, moment); + 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 - // only even terms have non-zero sums - if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) { - moment /= PS_SQR(sum); - } else { - // XXX keep this? - moment = 0.0; - } - - bool zeroNull = false; - float scale1D = 1.0 / sqrt(sum); - float scale2D = 1.0 / sum; - - if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) { - zeroNull = true; - } + 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 / max; + zeroNull = false; + } + } + if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) { - zeroNull = true; + zeroNull = true; } if (forceZeroNull) { - zeroNull = true; + 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)); + // Odd function + scale2D = 1.0 / max; + } + + 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->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32)); + } + psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32)); + penalty *= PS_SQR(scale2D); + if (zeroNull) { - preCalc->kernel->kernel[0][0] -= 1.0; + preCalc->kernel->kernel[0][0] -= 1.0; } @@ -245,9 +247,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); @@ -258,10 +260,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); + 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, fabsf(moment)); + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty); return true; @@ -284,7 +286,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]); } @@ -315,6 +317,6 @@ 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); + 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); } } @@ -325,6 +327,6 @@ pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(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); @@ -340,7 +342,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]); } @@ -353,5 +355,5 @@ 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 + num += (11 - gaussOrder - 1); // include all higher order radial terms } @@ -369,12 +371,12 @@ 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], 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); - } + 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; @@ -382,6 +384,6 @@ 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); @@ -397,7 +399,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]); } @@ -425,5 +427,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, false); + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false); } } @@ -434,6 +436,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); @@ -449,7 +451,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]); } @@ -488,20 +490,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, false); - - // 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 } @@ -512,17 +514,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); @@ -587,37 +589,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; + 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; + 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"); } @@ -626,9 +628,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; } @@ -958,5 +960,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 @@ -964,9 +966,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; @@ -1025,7 +1027,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, @@ -1117,8 +1119,8 @@ // currently known descriptions: - // ISIS(...), ISIS_RADIAL(...), HERM(...), DECONV_HERM(...), POIS(...), SPAM(...), - // FRIES(...), GUNK=ISIS(...)+POIS(...), RINGS(...), + // 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, '('); @@ -1130,26 +1132,26 @@ 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); + 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: @@ -1159,8 +1161,8 @@ PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); - 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."); - } + psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported."); + } return NULL; } @@ -1177,5 +1179,5 @@ char *ptr = strchr(type, '('); if (ptr) { - nameLength = ptr - type; + nameLength = ptr - type; }