Index: /branches/eam_branches/20091201/psModules/src/imcombine/Makefile.am =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/Makefile.am (revision 26478) +++ /branches/eam_branches/20091201/psModules/src/imcombine/Makefile.am (revision 26479) @@ -14,4 +14,5 @@ pmSubtractionKernels.c \ pmSubtractionHermitian.c \ + pmSubtractionDeconvolve.c \ pmSubtractionMask.c \ pmSubtractionMatch.c \ @@ -34,4 +35,5 @@ pmSubtractionKernels.h \ pmSubtractionHermitian.h \ + pmSubtractionDeconvolve.h \ pmSubtractionMask.h \ pmSubtractionMatch.h \ Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 26479) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 26479) @@ -0,0 +1,190 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmSubtractionDeconvolve.h" + +psKernel *pmSubtractionDeconvolveGauss (int size, float sigma) { + + psKernel *kernel = psKernelAlloc (-size, size, -size, size); + + // build the gaussian from 2 1-D Gaussians + psVector *vector = subtractionKernelISIS(sigma, 0, size); + + // generate 2D kernel, calculate moments + for (int v = -size, y = 0; v <= size; v++, y++) { + for (int u = -size, x = 0; u <= size; u++, x++) { + double value = vector->data.F32[x] * vector->data.F32[y]; // Value of kernel + kernel->kernel[v][u] = value; + } + } + + psFree (vector); + return kernel; +} + +// deconvolve kernelTarget by kernelConv to get the kernel which, when convolved +// by kernelConv results in kernelTarget... +psKernel *pmSubtractionDeconvolveKernel (psKernel *kernelTarg, psKernel *kernelConv) { + + PS_ASSERT_KERNEL_NON_NULL(kernelTarg, NULL); + PS_ASSERT_KERNEL_NON_NULL(kernelConv, NULL); + + bool threaded = psThreadPoolSize(); // Are we running threaded? + + // Size of image + int numCols = kernelConv->image->numCols; + int numRows = kernelConv->image->numRows; + + // kernel sizes + int xMin = kernelConv->xMin; + int xMax = kernelConv->xMax; + int yMin = kernelConv->yMin; + int yMax = kernelConv->yMax; + if (xMin != kernelTarg->xMin) goto escape; + if (xMax != kernelTarg->xMax) goto escape; + if (yMin != kernelTarg->yMin) goto escape; + if (yMax != kernelTarg->yMax) goto escape; + + int numPixels = numCols * numRows; // Number of pixels in padded image + + // operation is: Kt = FFT(kernelTarg), Kc = FFT(kernelConv) + // Kd = (Kt * Kc) / (Kc * Kc^*) + + // Create data array containing the image and kernel + FFTW_LOCK; + psF32 *dataTarg = fftwf_malloc(numPixels * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW + psF32 *dataConv = fftwf_malloc(numPixels * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW + FFTW_UNLOCK; + + // 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); + } + + // Do the forward FFTs + // 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); + FFTW_UNLOCK; + + fftwf_execute(forwardTarg); + fftwf_execute(forwardConv); + + FFTW_LOCK; + fftwf_destroy_plan(forwardTarg); + fftwf_destroy_plan(forwardConv); + FFTW_UNLOCK; + + // Combine the two transforms + // Targ = Tr + iTi, Conv = Cr + iCi + // Deco = Dr + iDi + // (Dr + i Di) = (Tr + iTi) / (Cr + iCi) + // (Dr + i Di) = (Tr + iTi) * (Cr - iCi) / (Cr^2 - Ci^2) + + // but anywhere Cr^2 - Ci^2 < 1e-7 of the max, mask it + + // generate Det = Cr^2 - Ci^2 + float maxValue = 0.0; + psImage *det = psImageAlloc(numCols, 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]; + det->data.F32[iy][ix] = convReal*convReal - convImag*convImag; + maxValue = PS_MAX(fabs(det->data.F32[iy][ix])); + } + } + 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]; + if (fabs(det->data.F32[iy][ix]) < limit) { + fftTarg[ix + iy*numCols][0] = 0.0; + fftTarg[ix + iy*numCols][1] = 0.0; + } else { + fftTarg[ix + iy*numCols][0] = targReal*convReal + targImag*convImag; + fftTarg[ix + iy*numCols][1] = targImag*convReal - targReal*convImag; + } + } + } + + // Do the backward FFT + FFTW_LOCK; + fftwf_plan backward = fftwf_plan_dft_c2r_2d(numRows, numCols, fftTarg, dataTarg, FFTW_PLAN_RIGOR); + FFTW_UNLOCK; + + fftwf_execute(backward); + + FFTW_LOCK; + fftwf_destroy_plan(backward); + fftwf_free(fftTarg); + fftwf_free(fftConv); + FFTW_UNLOCK; + + psKernel *output = psKernelAlloc (kernelTarg->xMin, kernelTarg->xMax, kernelTarg->yMin, kernelTarg->yMax); + for (int y = 0; y < numRows; y++, outData++, dataPtr += paddedCols) { + memcpy(output->image->data.F32[y], &dataTarg[y*numCols], numBytes); + } + + FFTW_LOCK; + fftwf_free(dataTarg); + fftwf_free(dataConv); + FFTW_UNLOCK; + + return output; +} + +bool pmSubtractionDeconvolutionTest () { + + int size = 15; + + // generate a Hermite polynomial + psVector *xKernel = subtractionKernelHERM(sigma, 2, size); // x Kernel + psVector *yKernel = subtractionKernelHERM(sigma, 2, size); // y Kernel + psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size); // Kernel + + // generate 2D kernel, calculate moments + for (int v = -size, y = 0; v <= size; v++, y++) { + for (int u = -size, x = 0; u <= size; u++, x++) { + double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel + kernelTarget->kernel[v][u] = value; + } + } + + // Gaussian convolution kernel + psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, 3.0); + + // deconvolve the target by the gaussian: + psKernel *kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel + + // re-convolve the kernel + psImage *kernelConv = psImageConvolveFFT(NULL, kernel, NULL, 0, kernelGauss); + pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv->image); + + return true; + +} Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.h (revision 26479) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.h (revision 26479) @@ -0,0 +1,15 @@ +#ifndef PM_SUBTRACTION_HERMITIAN_H +#define PM_SUBTRACTION_HERMITIAN_H + +/* these function support deconvolution operations used to generate deconvolved kernels. These + are kernels which, when convolved with the image sources, will yield a nearly orthonormal + basis set. The analysis starts with an orthonormal basis set (eg, Hermitian functions) and + deconvolves those basis functions with a Gaussian approximating the Gaussian of the image of + interest */ + +psKernel *pmSubtractionDeconvolveGauss (int size, float sigma); +psKernel *pmSubtractionDeconvolveKernel (psKernel *kernelTarg, psKernel *kernelConv); + +bool pmSubtractionDeconvolutionTest (); + +# endif Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26478) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26479) @@ -346,4 +346,125 @@ return kernels; } + +# if (0) +pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(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; + } + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECON_HERM, size, spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "DECON_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num); + psFree(params); + + // XXXXX hard-wired reference sigma for now of 1.7 pix + // generate the Gaussian deconvolution kernel + # define DECON_SIGMA 1.7 + psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECON_SIGMA); + + // Set the kernel parameters + int fullSize = 2 * size + 1; // Full size of kernels + 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++) { + 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 + psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size); // Kernel + + // generate 2D kernel, calculate moments + for (int v = -size, y = 0; v <= size; v++, y++) { + for (int u = -size, x = 0; u <= size; u++, x++) { + double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel + kernelTarget->kernel[v][u] = value; + } + } + + // deconvolve the target by the gaussian: + psKernel *kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel + preCalc->data[2] = kernel; + + psImage *kernelConv = psImageConvolveFFT(NULL, kernel, NULL, 0, kernelGauss); + pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv->image); + + // Normalise sum of kernel component to unity for even functions + if (uOrder % 2 == 0 && vOrder % 2 == 0) { + double sum = 0.0; // Sum of kernel component + for (int v = 0; v < fullSize; v++) { + for (int u = 0; u < fullSize; u++) { + sum += xKernel->data.F32[u] * yKernel->data.F32[v]; + } + } + sum = 1.0 / sqrt(sum); + psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); + psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); + psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); + +#if 1 + fprintf(stderr, "%d norm: %e, null: %e\n", index, sum, kernel->kernel[0][0]); +#endif + + kernel->kernel[0][0] -= 1.0; + moment *= PS_SQR(sum); + } + + +#if 1 + double sum = 0.0; // Sum of kernel component + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += kernel->kernel[v][u]; + } + } + fprintf(stderr, "%d sum: %e\n", index, sum); +#endif + + kernels->widths->data.F32[index] = fwhms->data.F32[i]; + 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; + kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); + + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, + fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); + } + } + } + + return kernels; +} +# endif ////////////////////////////////////////////////////////////////////////////////////////////////////////////// Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26478) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26479) @@ -11,4 +11,5 @@ PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels + PM_SUBTRACTION_KERNEL_DECON_HERM, ///< Deconvolved Hermitian polynomial kernels PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction @@ -31,8 +32,8 @@ psString description; ///< Description of the kernel parameters long num; ///< Number of kernel components (not including the spatial ones) - psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS or HERM) - psVector *widths; ///< Gaussian FWHMs (ISIS or HERM) + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECON_HERM) + psVector *widths; ///< Gaussian FWHMs (ISIS, HERM or DECON_HERM) psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS or HERM) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECON_HERM) float penalty; ///< Penalty for wideness psVector *penalties; ///< Penalty for each kernel component @@ -70,4 +71,9 @@ PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \ } \ + if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_DECON_HERM) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \ + } \ if ((KERNELS)->uStop || (KERNELS)->vStop) { \ PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \ @@ -157,4 +163,13 @@ ); +/// Generate DECON_HERM kernels +pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< order of hermitian polynomials + float penalty, ///< Penalty for wideness + pmSubtractionMode mode ///< Mode for subtraction + ); + /// Generate SPAM kernels pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel