Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26490) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (revision 26491) @@ -72,13 +72,11 @@ { int size = kernels->size; // Kernel half-size - psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values #if 0 - psVector *xKernel = preCalc->data[0]; // Kernel in x - psVector *yKernel = preCalc->data[1]; // Kernel in y // Iterating over the kernel for (int y = 0, v = -size; v <= size; y++, v++) { - float yValue = value * yKernel->data.F32[y]; + float yValue = value * preCalc->yKernel->data.F32[y]; for (int x = 0, u = -size; u <= size; x++, u++) { - kernel->kernel[v][u] += yValue * xKernel->data.F32[x]; + kernel->kernel[v][u] += yValue * preCalc->xKernel->data.F32[x]; } } @@ -88,8 +86,7 @@ } #else - psKernel *k = preCalc->data[2]; // Kernel image for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - kernel->kernel[v][u] += value * k->kernel[v][u]; + kernel->kernel[v][u] += value * preCalc->kernel->kernel[v][u]; } } @@ -168,14 +165,16 @@ break; } + case PM_SUBTRACTION_KERNEL_DECONV_HERM: { + solvedKernelISIS(kernel, kernels, value, i); + break; + } case PM_SUBTRACTION_KERNEL_RINGS: { - psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data - psVector *uCoords = preCalc->data[0]; // u coordinates - psVector *vCoords = preCalc->data[1]; // v coordinates - psVector *poly = preCalc->data[2]; // Polynomial values - int num = uCoords->n; // Number of pixels + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Precalculated kernels + int num = preCalc->uCoords->n; // Number of pixels for (int j = 0; j < num; j++) { - int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates - kernel->kernel[v][u] += poly->data.F32[j] * value; + int u = preCalc->uCoords->data.S32[j]; + int v = preCalc->vCoords->data.S32[j]; // Kernel coordinates + kernel->kernel[v][u] += preCalc->poly->data.F32[j] * value; } // Photometric scaling is built into the kernel --- no subtraction! @@ -466,9 +465,7 @@ ) { - psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data #if 0 // Convolving using separable convolution - psVector *xKernel = preCalc->data[0]; // Kernel in x - psVector *yKernel = preCalc->data[1]; // Kernel in y int size = kernels->size; // Size of kernel @@ -482,5 +479,5 @@ float value = 0.0; // Value of convolved pixel int uMin = x - size, uMax = x + size; // Range for u - psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values + psF32 *xKernelData = &preCalc->xKernel->data.F32[xKernel->n - 1]; // Kernel values psF32 *imageData = &image->kernel[y][uMin]; // Image values for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) { @@ -497,5 +494,5 @@ float value = 0.0; // Value of convolved pixel int vMin = y - size, vMax = y + size; // Range for v - psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values + psF32 *yKernelData = &preCalc->yKernel->data.F32[yKernel->n - 1]; // Kernel values psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way! for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) { @@ -514,5 +511,5 @@ #else // Convolving using precalculated kernel - return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]); + return p_pmSubtractionConvolveStampPrecalc(image, preCalc->kernel); #endif } @@ -687,14 +684,15 @@ return convolveStampISIS(image, kernels, index, footprint); } + case PM_SUBTRACTION_KERNEL_DECONV_HERM: { + return convolveStampISIS(image, kernels, index, footprint); + } case PM_SUBTRACTION_KERNEL_RINGS: { - psKernel *convolved = psKernelAlloc(-footprint, footprint, - -footprint, footprint); // Convolved image - psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data - psVector *uCoords = preCalc->data[0]; // u coordinates - psVector *vCoords = preCalc->data[1]; // v coordinates - psVector *poly = preCalc->data[2]; // Polynomial values - int num = uCoords->n; // Number of pixels - psS32 *uData = uCoords->data.S32, *vData = vCoords->data.S32; // Dereference u,v coordinates - psF32 *polyData = poly->data.F32; // Dereference polynomial values + psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data + + int num = preCalc->uCoords->n; // Number of pixels + psS32 *uData = preCalc->uCoords->data.S32; // Dereference v coordinate + psS32 *vData = preCalc->vCoords->data.S32; // Dereference u coordinate + psF32 *polyData = preCalc->poly->data.F32; // Dereference polynomial values psF32 **imageData = image->kernel; // Dereference image psF32 **convData = convolved->kernel; // Dereference convolved image Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 26490) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 26491) @@ -49,4 +49,5 @@ // deconvolve kernelTarget by kernelConv to get the kernel which, when convolved // by kernelConv results in kernelTarget... +// XXX using complex to complex, explicitly setting the imaginary part to zero psKernel *pmSubtractionDeconvolveKernel (psKernel *kernelTarg, psKernel *kernelConv) { @@ -77,17 +78,44 @@ // 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; - - size_t numBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row + // 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 + fftwf_complex *dataTarg = fftwf_malloc(numPixels * sizeof(fftwf_complex)); // Data for FFTW + fftwf_complex *dataConv = fftwf_malloc(numPixels * sizeof(fftwf_complex)); // Data for FFTW + 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 for (int y = 0; y < numRows; y++) { - memcpy(&dataTarg[y*numCols], kernelTarg->image->data.F32[y], numBytes); - } + for (int x = 0; x < numCols; x++) { + dataTarg[x + y*numCols][0] = kernelTarg->image->data.F32[y][x]; + dataTarg[x + y*numCols][1] = 0.0; + } + } + + // kernel must be copied to corners of image (0,0 pixel is center of kernel) // copy data from kernelConv image to dataConv array - for (int y = 0; y < numRows; y++) { - memcpy(&dataConv[y*numCols], kernelConv->image->data.F32[y], numBytes); + int oy = 0; + for (int iy = 0; iy <= yMax; iy++, oy++) { + int ox = 0; + for (int ix = 0; ix <= xMax; ix++, ox++) { + dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix]; + dataConv[ox + oy*numCols][1] = 0.0; + } + for (int ix = xMin; ix <= -1; ix++, ox++) { + dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix]; + dataConv[ox + oy*numCols][1] = 0.0; + } + } + for (int iy = yMin; iy <= -1; iy++, oy++) { + int ox = 0; + for (int ix = 0; ix <= xMax; ix++, ox++) { + dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix]; + dataConv[ox + oy*numCols][1] = 0.0; + } + for (int ix = xMin; ix <= -1; ix++, ox++) { + dataConv[ox + oy*numCols][0] = kernelConv->kernel[iy][ix]; + dataConv[ox + oy*numCols][1] = 0.0; + } } @@ -95,11 +123,11 @@ // Note that the FFT images have different size from the input FFTW_LOCK; - 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); + fftwf_complex *fftTarg = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // FFT + fftwf_complex *fftConv = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // FFT + FFTW_UNLOCK; + + FFTW_LOCK; + fftwf_plan forwardTarg = fftwf_plan_dft_2d(numRows, numCols, dataTarg, fftTarg, FFTW_FORWARD, FFTW_PLAN_RIGOR); + fftwf_plan forwardConv = fftwf_plan_dft_2d(numRows, numCols, dataConv, fftConv, FFTW_FORWARD, FFTW_PLAN_RIGOR); FFTW_UNLOCK; @@ -121,40 +149,75 @@ // the X dimension is halved by FFTW - int numColsOut = numCols / 2 + 1; + // int numColsOut = numCols / 2 + 1; // generate Det = Cr^2 - Ci^2 float maxValue = 0.0; - psImage *det = psImageAlloc(numColsOut, numRows, PS_TYPE_F32); - for (int iy = 0; iy < numRows; iy++) { - for (int ix = 0; ix < numColsOut; ix++) { - float convReal = fftConv[ix + iy*numColsOut][0]; - float convImag = fftConv[ix + iy*numColsOut][1]; + psImage *det = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImage *tR = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImage *tI = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImage *cR = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImage *cI = 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(maxValue, fabs(det->data.F32[iy][ix])); - } - } + + tR->data.F32[iy][ix] = fftTarg[ix + iy*numCols][0]; + tI->data.F32[iy][ix] = fftTarg[ix + iy*numCols][1]; + cR->data.F32[iy][ix] = fftConv[ix + iy*numCols][0]; + cI->data.F32[iy][ix] = fftConv[ix + iy*numCols][1]; + } + } + + // pmSubtractionVisualShowSubtraction (det, tR, tI); + // pmSubtractionVisualShowSubtraction (det, cR, cI); + +# if 1 # 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 < 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]; + 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*numColsOut][0] = 0.0; - fftTarg[ix + iy*numColsOut][1] = 0.0; + fftTarg[ix + iy*numCols][0] = 0.0; + fftTarg[ix + iy*numCols][1] = 0.0; } else { - fftTarg[ix + iy*numColsOut][0] = targReal*convReal + targImag*convImag; - fftTarg[ix + iy*numColsOut][1] = targImag*convReal - targReal*convImag; + fftTarg[ix + iy*numCols][0] = (targReal*convReal + targImag*convImag) / det->data.F32[iy][ix]; + fftTarg[ix + iy*numCols][1] = (targImag*convReal - targReal*convImag) / det->data.F32[iy][ix]; + // fftTarg[ix + iy*numCols][0] = (targReal*convReal + targImag*convImag); + // fftTarg[ix + iy*numCols][1] = (targImag*convReal - targReal*convImag); } } } +# else + 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]; + fftTarg[ix + iy*numCols][0] = targReal*convReal - targImag*convImag; + fftTarg[ix + iy*numCols][1] = targImag*convReal + targReal*convImag; + } + } +# endif + + for (int iy = 0; iy < numRows; iy++) { + for (int ix = 0; ix < numCols; ix++) { + tR->data.F32[iy][ix] = fftTarg[ix + iy*numCols][0]; + tI->data.F32[iy][ix] = fftTarg[ix + iy*numCols][1]; + } + } + // pmSubtractionVisualShowSubtraction (det, tR, tI); // Do the backward FFT FFTW_LOCK; - fftwf_plan backward = fftwf_plan_dft_c2r_2d(numRows, numCols, fftTarg, dataTarg, FFTW_PLAN_RIGOR); + fftwf_plan backward = fftwf_plan_dft_2d(numRows, numCols, fftTarg, dataTarg, FFTW_BACKWARD, FFTW_PLAN_RIGOR); FFTW_UNLOCK; @@ -169,5 +232,7 @@ psKernel *output = psKernelAlloc (kernelTarg->xMin, kernelTarg->xMax, kernelTarg->yMin, kernelTarg->yMax); for (int y = 0; y < numRows; y++) { - memcpy(output->image->data.F32[y], &dataTarg[y*numCols], numBytes); + for (int x = 0; x < numCols; x++) { + output->image->data.F32[y][x] = dataTarg[x + y*numCols][0]; + } } @@ -185,12 +250,12 @@ } -bool pmSubtractionDeconvolutionTest () { - - float sigma = 2.0; - int size = 15; +bool pmSubtractionDeconvolutionTest (int order) { + + float sigma = 1.0; + int size = 31; // generate a Hermite polynomial - psVector *xKernel = pmSubtractionKernelHERM(sigma, 2, size); // x Kernel - psVector *yKernel = pmSubtractionKernelHERM(sigma, 2, size); // y Kernel + psVector *xKernel = pmSubtractionKernelHERM(sigma, order, size); // x Kernel + psVector *yKernel = pmSubtractionKernelHERM(sigma, order, size); // y Kernel psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size); // Kernel Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.h (revision 26490) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionDeconvolve.h (revision 26491) @@ -1,4 +1,4 @@ -#ifndef PM_SUBTRACTION_HERMITIAN_H -#define PM_SUBTRACTION_HERMITIAN_H +#ifndef PM_SUBTRACTION_DECONVOLVE_H +#define PM_SUBTRACTION_DECONVOLVE_H /* these function support deconvolution operations used to generate deconvolved kernels. These @@ -11,5 +11,5 @@ psKernel *pmSubtractionDeconvolveKernel (psKernel *kernelTarg, psKernel *kernelConv); -bool pmSubtractionDeconvolutionTest (); +bool pmSubtractionDeconvolutionTest (int order); # endif Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26490) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (revision 26491) @@ -11,4 +11,6 @@ #include "pmSubtractionKernels.h" #include "pmSubtractionHermitian.h" +#include "pmSubtractionDeconvolve.h" +#include "pmSubtractionVisual.h" #define RINGS_BUFFER 10 // Buffer size for RINGS data @@ -27,4 +29,16 @@ psFree(kernels->solution1); psFree(kernels->solution2); +} + +// Free function for pmSubtractionPreCalcKernel +static void pmSubtractionKernelPreCalcFree(pmSubtractionKernelPreCalc *kernel) +{ + psFree(kernel->xKernel); + psFree(kernel->yKernel); + psFree(kernel->kernel); + + psFree(kernel->uCoords); + psFree(kernel->vCoords); + psFree(kernel->poly); } @@ -132,4 +146,61 @@ } +bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm) { + + // 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 += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); + } + } + + // Normalize 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 = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += preCalc->kernel->kernel[v][u]; + } + } + sum = 1.0 / sqrt(sum); + psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); + psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32)); + psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32)); + +#if 1 + fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum, preCalc->kernel->kernel[0][0]); +#endif + + preCalc->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 += preCalc->kernel->kernel[v][u]; + } + } + fprintf(stderr, "%d sum: %lf\n", index, sum); +#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; + 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)); + + return true; +} + pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, const psVector *fwhmsIN, const psVector *ordersIN, @@ -172,5 +243,4 @@ // 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 @@ -178,57 +248,7 @@ 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] = 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 - - // Calculate moments - double moment = 0.0; // Moment, for penalty - 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 - kernel->kernel[v][u] = value; - moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); - } - } - - // 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)); - kernel->kernel[0][0] -= 1.0; - moment *= PS_SQR(sum); - } - - -#if 0 - 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: %lf\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)); + + 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]); } } @@ -276,5 +296,4 @@ // 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 @@ -282,62 +301,6 @@ 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] = 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 - - // Calculate moments - double moment = 0.0; // Moment, for penalty - 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 - kernel->kernel[v][u] = value; - moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v))); - } - } - - // 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: %lf, null: %f\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: %lf\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)); + 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]); } } @@ -347,6 +310,5 @@ } -# if (0) -pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, int spatialOrder, +pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder, const psVector *fwhmsIN, const psVector *ordersIN, float penalty, pmSubtractionMode mode) @@ -379,17 +341,16 @@ } - 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); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, spatialOrder, penalty, mode); // The kernels + psStringAppend(&kernels->description, "DECONV_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 + // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest) // generate the Gaussian deconvolution kernel - # define DECON_SIGMA 1.7 - psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECON_SIGMA); + # define DECONV_SIGMA 1.7 + psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECONV_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 @@ -397,67 +358,16 @@ 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] = 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 - - // 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)); + + 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; + preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel + + // XXXX test demo that deconvolved kernel is valid + // psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss); + // pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv); + + pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]); } } @@ -466,5 +376,4 @@ return kernels; } -# endif ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -503,4 +412,58 @@ } +pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(pmSubtractionKernelsType type, int uOrder, int vOrder, int size, float sigma) { + + pmSubtractionKernelPreCalc *preCalc = psAlloc(sizeof(pmSubtractionKernelPreCalc)); // Kernels, to return + psMemSetDeallocator(preCalc, (psFreeFunc)pmSubtractionKernelPreCalcFree); + + // 1D kernel realizations: + 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; + 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; + 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; + default: + psAbort("programming error: invalid type for PreCalc kernel"); + } + + preCalc->kernel = psKernelAlloc(-size, size, -size, size); // 2D Kernel + + // 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 + } + } + + return preCalc; +} + pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, pmSubtractionMode mode) @@ -827,15 +790,15 @@ for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) { - psArray *data = psArrayAlloc(3); // Container for data - psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords - psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords - psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc (PM_SUBTRACTION_KERNEL_RINGS, 0, 0, RINGS_BUFFER, 0.0); double moment = 0.0; // Moment, for penalty if (i == 0) { // Central pixel is easy - uCoords->data.S32[0] = vCoords->data.S32[0] = 0; - poly->data.F32[0] = 1.0; - uCoords->n = vCoords->n = poly->n = 1; + preCalc->uCoords->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; radiusLast = 0; moment = 0.0; @@ -855,15 +818,15 @@ float polyVal = uPoly * vPoly; // Value of polynomial if (polyVal != 0) { // No point adding it otherwise - uCoords->data.S32[j] = u; - vCoords->data.S32[j] = v; - poly->data.F32[j] = polyVal; + preCalc->uCoords->data.S32[j] = u; + preCalc->vCoords->data.S32[j] = v; + preCalc->poly->data.F32[j] = polyVal; norm += polyVal; moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v)); - psVectorExtend(uCoords, RINGS_BUFFER, 1); - psVectorExtend(vCoords, RINGS_BUFFER, 1); - psVectorExtend(poly, RINGS_BUFFER, 1); + psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->poly, RINGS_BUFFER, 1); psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n", - u, v, poly->data.F32[j]); + u, v, preCalc->poly->data.F32[j]); j++; } @@ -873,22 +836,22 @@ // Normalise kernel component to unit sum if (uOrder % 2 == 0 && vOrder % 2 == 0) { - psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); + psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); // Add subtraction of 0,0 component to preserve photometric scaling - uCoords->data.S32[j] = 0; - vCoords->data.S32[j] = 0; - poly->data.F32[j] = -1.0; - psVectorExtend(uCoords, RINGS_BUFFER, 1); - psVectorExtend(vCoords, RINGS_BUFFER, 1); - psVectorExtend(poly, RINGS_BUFFER, 1); + preCalc->uCoords->data.S32[j] = 0; + preCalc->vCoords->data.S32[j] = 0; + preCalc->poly->data.F32[j] = -1.0; + psVectorExtend(preCalc->uCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->vCoords, RINGS_BUFFER, 1); + psVectorExtend(preCalc->poly, RINGS_BUFFER, 1); } else { norm = powf(size, uOrder) * powf(size, vOrder); - psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); + psBinaryOp(preCalc->poly, preCalc->poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); } moment /= norm; } - psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", uCoords->n); - - kernels->preCalc->data[index] = data; + psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", preCalc->uCoords->n); + + kernels->preCalc->data[index] = preCalc; kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; @@ -916,4 +879,6 @@ case PM_SUBTRACTION_KERNEL_HERM: return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode); + case PM_SUBTRACTION_KERNEL_DECONV_HERM: + return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode); case PM_SUBTRACTION_KERNEL_SPAM: return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); @@ -978,6 +943,6 @@ float penalty = 0.0; // Penalty for wideness - // ISIS and HERM have the same description layout - if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM")) { + // 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")) { @@ -987,7 +952,18 @@ type = pmSubtractionKernelsTypeFromString (description); - psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must be ISIS or HERM"); - - char *ptr = (char*)description + 5; // Eat "ISIS(" or "HERM(" + 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); @@ -1025,5 +1001,5 @@ } - psAbort("Deciphering kernels other than ISIS, HERM and RINGS is not currently supported."); + psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported."); return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, @@ -1042,4 +1018,7 @@ if (strcasecmp(type, "HERM") == 0) { return PM_SUBTRACTION_KERNEL_HERM; + } + if (strcasecmp(type, "DECONV_HERM") == 0) { + return PM_SUBTRACTION_KERNEL_DECONV_HERM; } if (strcasecmp(type, "SPAM") == 0) { Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26490) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h (revision 26491) @@ -11,5 +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_DECONV_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 @@ -32,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, HERM or DECON_HERM) - psVector *widths; ///< Gaussian FWHMs (ISIS, HERM or DECON_HERM) + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM) + psVector *widths; ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM) psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECON_HERM) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM) float penalty; ///< Penalty for wideness psVector *penalties; ///< Penalty for each kernel component @@ -49,4 +49,15 @@ int numStamps; ///< Number of good stamps } pmSubtractionKernels; + +// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures +typedef struct { + psVector *uCoords; // used by RINGS + psVector *vCoords; // used by RINGS + psVector *poly; // used by RINGS + + psVector *xKernel; // used by ISIS, HERM, DECONV_HERM + psVector *yKernel; // used by ISIS, HERM, DECONV_HERM + psKernel *kernel; // used by ISIS, HERM, DECONV_HERM +} pmSubtractionKernelPreCalc; // Assertion to check pmSubtractionKernels @@ -71,5 +82,5 @@ PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \ } \ - if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_DECON_HERM) { \ + if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_DECONV_HERM) { \ PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \ PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \ @@ -139,4 +150,14 @@ ); +/// Allocator for pre-calculated kernel data structure +pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc( + pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated) + int uOrder, ///< order in x-direction + int vOrder, ///< order in x-direction + int size, ///< Half-size of the kernel + float sigma ///< sigma of gaussian kernel + ); + + /// Generate POIS kernels pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims) @@ -173,6 +194,6 @@ ); -/// Generate DECON_HERM kernels -pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, ///< Half-size of the kernel +/// Generate DECONV_HERM kernels +pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, ///< Half-size of the kernel int spatialOrder, ///< Order of spatial variations const psVector *fwhms, ///< Gaussian FWHMs Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.c (revision 26490) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.c (revision 26491) @@ -482,9 +482,18 @@ // Maintain photometric scaling if (type == PM_SUBTRACTION_KERNEL_ISIS) { - psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest + + // XXX in r26035, this code was just wrong. we had: + + // psKernel *subtract = kernels->preCalc->data[0] + + // but, kernels->preCalc was an array of psArray, not an array of kernels. It is now + // an array of pmSubtractionKernelPreCalc. + + pmSubtractionKernelPreCalc *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest + for (int i = 1; i < newSize; i++) { if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { - psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest - psBinaryOp(kernel->image, kernel->image, "-", subtract->image); + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest + psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "-", subtract->kernel->image); } } @@ -495,6 +504,6 @@ for (int i = 0; i < newSize; i++) { if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { - psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest - kernel->kernel[0][0] -= 1.0; + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Kernel of interest + preCalc->kernel->kernel[0][0] -= 1.0; } } Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26490) +++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c (revision 26491) @@ -242,6 +242,6 @@ for (int i = 0; i < kernels->num; i++) { - psArray *preCalc = kernels->preCalc->data[i]; - psKernel *kernel = preCalc->data[2]; + pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; + psKernel *kernel = preCalc->kernel; int xSub = i % NXsub;