Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 26485) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 26486) @@ -6,7 +6,25 @@ #include #include +#include #include +#include "pmFPA.h" +#include "pmSubtractionKernels.h" #include "pmSubtractionDeconvolve.h" +#include "pmSubtractionStamps.h" +#include "pmSubtractionVisual.h" + +// Lock FFTW access +#define FFTW_LOCK \ +if (threaded) { \ + psFFTLock(); \ +} +// Unlock FFTW access +#define FFTW_UNLOCK \ +if (threaded) { \ + psFFTUnlock(); \ +} + +#define FFTW_PLAN_RIGOR FFTW_ESTIMATE // How rigorous the FFTW planning is psKernel *pmSubtractionDeconvolveGauss (int size, float sigma) { @@ -15,5 +33,5 @@ // build the gaussian from 2 1-D Gaussians - psVector *vector = subtractionKernelISIS(sigma, 0, size); + psVector *vector = pmSubtractionKernelISIS(sigma, 0, size); // generate 2D kernel, calculate moments @@ -63,12 +81,11 @@ FFTW_UNLOCK; + size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row + // copy data from kernelTarg image to dataTarg array - size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row for (int y = 0; y < numRows; y++) { memcpy(&dataTarg[y*numCols], kernelTarg->image->data.F32[y], numBytes); } - // copy data from kernelConv image to dataConv array - size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row for (int y = 0; y < numRows; y++) { memcpy(&dataConv[y*numCols], kernelConv->image->data.F32[y], numBytes); @@ -78,11 +95,11 @@ // Note that the FFT images have different size from the input FFTW_LOCK; - fftwf_complex *fftTarg = fftwf_malloc((numCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT - fftwf_complex *fftConv = fftwf_malloc((numCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT - FFTW_UNLOCK; - - FFTW_LOCK; - fftwf_plan forwardTarg = fftwf_plan_dft_r2c_2d(numRows, numCols, dataTarg, fftTarg, FFT_PLAN_RIGOR); - fftwf_plan forwardConv = fftwf_plan_dft_r2c_2d(numRows, numCols, dataConv, fftConv, FFT_PLAN_RIGOR); + fftwf_complex *fftTarg = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // FFT + fftwf_complex *fftConv = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // FFT + FFTW_UNLOCK; + + FFTW_LOCK; + fftwf_plan forwardTarg = fftwf_plan_dft_r2c_2d(numRows, numCols, dataTarg, fftTarg, FFTW_PLAN_RIGOR); + fftwf_plan forwardConv = fftwf_plan_dft_r2c_2d(numRows, numCols, dataConv, fftConv, FFTW_PLAN_RIGOR); FFTW_UNLOCK; @@ -103,30 +120,34 @@ // but anywhere Cr^2 - Ci^2 < 1e-7 of the max, mask it + // the X dimension is halved by FFTW + int numColsOut = numCols / 2 + 1; + // generate Det = Cr^2 - Ci^2 float maxValue = 0.0; - psImage *det = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImage *det = psImageAlloc(numColsOut, numRows, PS_TYPE_F32); for (int iy = 0; iy < numRows; iy++) { - for (int ix = 0; ix < numCols; ix++) { - float convReal = fftConv[ix + iy*numCols][0]; - float convImag = fftConv[ix + iy*numCols][1]; + for (int ix = 0; ix < numColsOut; ix++) { + float convReal = fftConv[ix + iy*numColsOut][0]; + float convImag = fftConv[ix + iy*numColsOut][1]; det->data.F32[iy][ix] = convReal*convReal - convImag*convImag; - maxValue = PS_MAX(fabs(det->data.F32[iy][ix])); - } - } + maxValue = PS_MAX(maxValue, fabs(det->data.F32[iy][ix])); + } + } +# define TOL 1e-7 float limit = TOL*maxValue; // generate Deco = targ * conv^* / (Cr^2 - Ci^2) for (int iy = 0; iy < numRows; iy++) { - for (int ix = 0; ix < numCols; ix++) { - float targReal = fftTarg[ix + iy*numCols][0]; - float targImag = fftTarg[ix + iy*numCols][1]; - float convReal = fftConv[ix + iy*numCols][0]; - float convImag = fftConv[ix + iy*numCols][1]; + for (int ix = 0; ix < numColsOut; ix++) { + float targReal = fftTarg[ix + iy*numColsOut][0]; + float targImag = fftTarg[ix + iy*numColsOut][1]; + float convReal = fftConv[ix + iy*numColsOut][0]; + float convImag = fftConv[ix + iy*numColsOut][1]; if (fabs(det->data.F32[iy][ix]) < limit) { - fftTarg[ix + iy*numCols][0] = 0.0; - fftTarg[ix + iy*numCols][1] = 0.0; + fftTarg[ix + iy*numColsOut][0] = 0.0; + fftTarg[ix + iy*numColsOut][1] = 0.0; } else { - fftTarg[ix + iy*numCols][0] = targReal*convReal + targImag*convImag; - fftTarg[ix + iy*numCols][1] = targImag*convReal - targReal*convImag; + fftTarg[ix + iy*numColsOut][0] = targReal*convReal + targImag*convImag; + fftTarg[ix + iy*numColsOut][1] = targImag*convReal - targReal*convImag; } } @@ -147,5 +168,5 @@ psKernel *output = psKernelAlloc (kernelTarg->xMin, kernelTarg->xMax, kernelTarg->yMin, kernelTarg->yMax); - for (int y = 0; y < numRows; y++, outData++, dataPtr += paddedCols) { + for (int y = 0; y < numRows; y++) { memcpy(output->image->data.F32[y], &dataTarg[y*numCols], numBytes); } @@ -157,13 +178,19 @@ return output; + + escape: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "mismatch between kernel and image"); + return NULL; + } bool pmSubtractionDeconvolutionTest () { + float sigma = 2.0; int size = 15; // generate a Hermite polynomial - psVector *xKernel = subtractionKernelHERM(sigma, 2, size); // x Kernel - psVector *yKernel = subtractionKernelHERM(sigma, 2, size); // y Kernel + psVector *xKernel = pmSubtractionKernelHERM(sigma, 2, size); // x Kernel + psVector *yKernel = pmSubtractionKernelHERM(sigma, 2, size); // y Kernel psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size); // Kernel @@ -183,6 +210,6 @@ // re-convolve the kernel - psImage *kernelConv = psImageConvolveFFT(NULL, kernel, NULL, 0, kernelGauss); - pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv->image); + psImage *kernelConv = psImageConvolveFFT(NULL, kernel->image, NULL, 0, kernelGauss); + pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv); return true; Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26485) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26486) @@ -45,5 +45,5 @@ // 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 @@ -63,5 +63,5 @@ // Generate 1D convolution kernel for HERM (normalized for 2D) -static psVector *subtractionKernelHERM(float sigma, // Gaussian width +psVector *pmSubtractionKernelHERM(float sigma, // Gaussian width int order, // Polynomial order int size // Kernel half-size @@ -179,6 +179,6 @@ for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { psArray *preCalc = psArrayAlloc(3); // 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 + psVector *xKernel = preCalc->data[0] = pmSubtractionKernelISIS(sigma, uOrder, size); // x Kernel + psVector *yKernel = preCalc->data[1] = pmSubtractionKernelISIS(sigma, vOrder, size); // y Kernel psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size); // Kernel @@ -283,6 +283,6 @@ for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values - psVector *xKernel = preCalc->data[0] = subtractionKernelHERM(sigma, uOrder, size); // x Kernel - psVector *yKernel = preCalc->data[1] = subtractionKernelHERM(sigma, vOrder, size); // y Kernel + psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel + psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size); // Kernel @@ -398,6 +398,6 @@ for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values - psVector *xKernel = preCalc->data[0] = subtractionKernelHERM(sigma, uOrder, size); // x Kernel - psVector *yKernel = preCalc->data[1] = subtractionKernelHERM(sigma, vOrder, size); // y Kernel + psVector *xKernel = preCalc->data[0] = pmSubtractionKernelHERM(sigma, uOrder, size); // x Kernel + psVector *yKernel = preCalc->data[1] = pmSubtractionKernelHERM(sigma, vOrder, size); // y Kernel psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size); // Kernel Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26485) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26486) @@ -111,4 +111,14 @@ } +psVector *pmSubtractionKernelISIS(float sigma, // Gaussian width + int order, // Polynomial order + int size // Kernel half-size + ); + +psVector *pmSubtractionKernelHERM(float sigma, // Gaussian width + int order, // Polynomial order + int size // Kernel half-size + ); + /// Generate a delta-function grid for subtraction kernels (like the POIS kernel) bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, ///< The subtraction kernels to append to Index: /branches/eam_branches/20091201/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/psmodules.h (revision 26485) +++ /branches/eam_branches/20091201/psModules/src/psmodules.h (revision 26486) @@ -96,4 +96,5 @@ #include #include +#include #include #include