Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 29596) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 29597) @@ -183,10 +183,7 @@ } -# define CENTRAL_DELTA 0 - # if (CENTRAL_DELTA) - -// XXX *** this code used the central pixel to force zero net flux, -// Alard actually uses kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0]) +// This version of the code uses the central pixel to force zero net flux, Alard actually uses +// kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0]) static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int uOrder, int vOrder, float fwhm, @@ -301,5 +298,5 @@ static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int uOrder, int vOrder, float fwhm, - bool AlardLuptonStyle, bool forceZeroNull) + pmSubtractionKernelPreCalc *zeroKernel) { // 1) for odd functions: no renormalization @@ -339,9 +336,13 @@ } } - if (index == 0) { +# if (ZERO_KERNEL_ZERO_FLUX) + int firstZeroIndex = 0; +# else + int firstZeroIndex = 1; +# endif + if (index < firstZeroIndex) { zeroIsNormal = true; } else { psAssert(zeroIsNormal, "failed to normalize zero kernel first"); - pmSubtractionKernelPreCalc *zeroKernel = kernels->preCalc->data[0]; psAssert(zeroKernel, "failed to supply zero kernel"); for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) { @@ -379,11 +380,13 @@ #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; + if (kernels) { + 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; + } psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder); @@ -431,7 +434,23 @@ psFree(params); +# if (!CENTRAL_DELTA && !ZERO_KERNEL_ZERO_FLUX) + // in this case, subtract the 0th kernel from everyone else + float zeroFWHM = fwhms->data.F32[0]; + float zeroSigma = zeroFWHM / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma + pmSubtractionKernelPreCalc *zeroKernel = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, 0, 0, size, zeroSigma); // structure to hold precalculated values + pmSubtractionKernelPreCalcNormalize (NULL, zeroKernel, -1, 0, 0, zeroFWHM, NULL); +# 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 + +# if (!CENTRAL_DELTA && ZERO_KERNEL_ZERO_FLUX) + // in this case, subtract the a 1/2 width Gaussian from each series + float zeroFWHM = 0.50 * fwhms->data.F32[0]; + float zeroSigma = zeroFWHM / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma + pmSubtractionKernelPreCalc *zeroKernel = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, 0, 0, size, zeroSigma); // structure to hold precalculated values + pmSubtractionKernelPreCalcNormalize (NULL, zeroKernel, -1, 0, 0, zeroFWHM, NULL); +# endif // Iterate over (u,v) order for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { @@ -439,10 +458,19 @@ pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values - pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); - // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false); +# if (CENTRAL_DELTA) + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false); +# else + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], zeroKernel); +# endif } } - } - +# if (!CENTRAL_DELTA && ZERO_KERNEL_ZERO_FLUX) + psFree(zeroKernel); +# endif + } + +# if (!CENTRAL_DELTA && !ZERO_KERNEL_ZERO_FLUX) + psFree(zeroKernel); +# endif psFree(orders); psFree(fwhms); @@ -497,5 +525,9 @@ 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 +# if (CENTRAL_DELTA) pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); +# else + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], NULL); +# endif } } @@ -503,5 +535,9 @@ // 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 +# if (CENTRAL_DELTA) pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true); +# else + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], NULL); +# endif } } @@ -555,5 +591,9 @@ 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 +# if (CENTRAL_DELTA) pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); +# else + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], NULL); +# endif } } @@ -624,6 +664,9 @@ // XXX do we use Alard-Lupton normalization (last param true) or not? +# if (CENTRAL_DELTA) pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false); - +# else + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], NULL); +# endif // XXXX test demo that deconvolved kernel is valid # if 1 @@ -1396,2 +1439,61 @@ return out; } + +#define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image +psImage *pmSubtractionKernelsImageMosaic(pmSubtractionKernels *kernels) { + + psTrace("psModules.imcombine", 2, "Generating diagnostic image..\n"); + + // Generate image with convolution kernels + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1 + 1; // Full size of kernel + int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize; + psImage *convKernels = psImageAlloc((kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) * + imageSize - 1 + + (kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0), + imageSize - 1, PS_TYPE_F32); + psImageInit(convKernels, NAN); + for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) { + for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) { + psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, + (float)j / (float)KERNEL_MOSAIC, + false); // Image of the kernel + if (!kernel) { + psError(psErrorCodeLast(), false, "Unable to generate kernel image."); + psFree(convKernels); + return NULL; + } + + if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize, + (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { + psError(psErrorCodeLast(), false, "Unable to overlay kernel image."); + psFree(kernel); + psFree(convKernels); + return NULL; + } + psFree(kernel); + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, + (float)j / (float)KERNEL_MOSAIC, + true); // Image of the kernel + if (!kernel) { + psError(psErrorCodeLast(), false, "Unable to generate kernel image."); + psFree(convKernels); + return NULL; + } + + if (psImageOverlaySection(convKernels, kernel, + (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize + 4, + (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { + psError(psErrorCodeLast(), false, "Unable to overlay kernel image."); + psFree(kernel); + psFree(convKernels); + return NULL; + } + psFree(kernel); + } + } + } + return convKernels; +}