Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14708) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14709) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-30 03:50:28 $ + * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-30 21:45:26 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -249,17 +249,4 @@ } -// Generate the convolution given a precalculated kernel -static psKernel *convolvePrecalc(const psKernel *image, // Image to convolve (a kernel for convenience) - const psKernel *kernel, // Kernel by which to convolve - int footprint // Size of region of interest - ) -{ - psImage *conv = psImageConvolveFFT(image->image, kernel, 0.0); // Convolved image - int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image - psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version - psFree(conv); - return convolved; -} - // Generate the convolved reference image static psKernel *convolveRef(const pmSubtractionKernels *kernels, // Kernel basis functions @@ -308,5 +295,5 @@ if (index < kernels->inner) { // Photometric scaling is already built in to the precalculated kernel - return convolvePrecalc(image, kernels->preCalc->data[index], footprint); + return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]); } // Using delta function @@ -321,5 +308,5 @@ case PM_SUBTRACTION_KERNEL_ISIS: { // Photometric scaling is already built in to the precalculated kernel - return convolvePrecalc(image, kernels->preCalc->data[index], footprint); + return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]); } case PM_SUBTRACTION_KERNEL_RINGS: { @@ -355,4 +342,62 @@ } +// Convolve an image using FFT +static void convolveFFT(psImage *target,// Place the result in here + const psImage *image, // Image to convolve + const psKernel *kernel, // Kernel by which to convolve + psRegion region,// Region of interest + float background, // Background to add + int size // Size of (square) kernel + ) +{ + psRegion border = psRegionSet(region.x0 - size, region.x1 + size, + region.y0 - size, region.y1 + size); // Add a border + + // Casting away const so psImageSubset can add the child + psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve + psImage *convolved = psImageConvolveFFT(subImage, kernel, 0.0); // Convolution + psFree(subImage); + + // Now, we have to chop off the borders, and stick it in where it belongs + psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, + size, -size)); // Cut off the edges + psImage *subTarget = psImageSubset(target, region); // Target for this subregion + if (background != 0.0) { + psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32)); + } else { + int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy + for (int y = 0; y < subTarget->numRows; y++) { + memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes); + } + } + psFree(subConv); + psFree(convolved); + psFree(subTarget); + return; +} + +// Convolve an image directly +static void convolveDirect(psImage *target, // Put the result here + const psImage *image, // Image to convolve + const psKernel *kernel, // Kernel by which to convolve + psRegion region,// Region of interest + float background, // Background to add + int size // Size of (square) kernel + ) +{ + for (int y = region.y0; y < region.y1; y++) { + for (int x = region.x0; x < region.x1; x++) { + target->data.F32[y][x] = background; + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + target->data.F32[y][x] += kernel->kernel[v][u] * + image->data.F32[y - v][x - u]; + } + } + } + } + return; +} + // Mark a pixel as blank in the image, mask and weight static inline void markBlank(psImage *image, // Image to mark as blank @@ -372,4 +417,22 @@ return; } + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Semi-public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Generate the convolution given a precalculated kernel +psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel) +{ + PS_ASSERT_KERNEL_NON_NULL(image, NULL); + PS_ASSERT_KERNEL_NON_NULL(kernel, NULL); + + psImage *conv = psImageConvolveFFT(image->image, kernel, 0.0); // Convolved image + int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image + psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version + psFree(conv); + return convolved; +} + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -916,62 +979,4 @@ } -// Convolve an image using FFT -static void convolveFFT(psImage *target,// Place the result in here - const psImage *image, // Image to convolve - const psKernel *kernel, // Kernel by which to convolve - psRegion region,// Region of interest - float background, // Background to add - int size // Size of (square) kernel - ) -{ - psRegion border = psRegionSet(region.x0 - size, region.x1 + size, - region.y0 - size, region.y1 + size); // Add a border - - // Casting away const so psImageSubset can add the child - psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve - psImage *convolved = psImageConvolveFFT(subImage, kernel, 0.0); // Convolution - psFree(subImage); - - // Now, we have to chop off the borders, and stick it in where it belongs - psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, - size, -size)); // Cut off the edges - psImage *subTarget = psImageSubset(target, region); // Target for this subregion - if (background != 0.0) { - psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32)); - } else { - int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy - for (int y = 0; y < subTarget->numRows; y++) { - memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes); - } - } - psFree(subConv); - psFree(convolved); - psFree(subTarget); - return; -} - -// Convolve an image directly -static void convolveDirect(psImage *target, // Put the result here - const psImage *image, // Image to convolve - const psKernel *kernel, // Kernel by which to convolve - psRegion region,// Region of interest - float background, // Background to add - int size // Size of (square) kernel - ) -{ - for (int y = region.y0; y < region.y1; y++) { - for (int x = region.x0; x < region.x1; x++) { - target->data.F32[y][x] = background; - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - target->data.F32[y][x] += kernel->kernel[v][u] * - image->data.F32[y - v][x - u]; - } - } - } - } - return; -} - bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, const psImage *inImage, const psImage *inWeight, const psImage *subMask, Index: /trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.h (revision 14708) +++ /trunk/psModules/src/imcombine/pmSubtraction.h (revision 14709) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-30 03:50:28 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-30 21:45:26 $ * * Copyright 2004-207 Institute for Astronomy, University of Hawaii @@ -93,4 +93,11 @@ ); +/// Generate the convolution of an image, given a precalculated kernel +/// +/// The 'image' is a kernel for convenience --- intended to be a stamp +psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, ///< Image to convolve + const psKernel *kernel ///< Kernel by which to convolve + ); + /// @} #endif Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14708) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14709) @@ -216,5 +216,5 @@ // Subtract a reference kernel from all the others to maintain flux scaling -// if (spatialOrder > 0) { + if (spatialOrder > 0) { int subIndex = 0; // Index of kernel to subtract from others psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others @@ -230,5 +230,5 @@ } } -// } + } if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { @@ -454,5 +454,5 @@ // Subtract unity from the kernels to maintain photometric flux scaling int numGaussians = fwhms->n; // Number of Gaussians -// if (spatialOrder > 0) { + if (spatialOrder > 0) { for (int i = 0, index = 0; i < numGaussians; i++) { for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { @@ -465,5 +465,5 @@ } } -// } + } int numISIS = kernels->num; // Number of ISIS kernels Index: /trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 14708) +++ /trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 14709) @@ -9,7 +9,10 @@ #include "pmSubtractionStamps.h" +#include "pmSubtraction.h" +#include "pmSubtractionParams.h" ////////////////////////////////////////////////////////////////////////////////////////////////////////////// +#if 0 // Convolve the reference stamp by the kernel static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved @@ -45,5 +48,5 @@ return convolution; } - +#endif // Accumulate cross-term sums for a stamp @@ -53,20 +56,16 @@ const pmSubtractionStamp *stamp, // Stamp with weight const psKernel *input, // Input image, I(x) - const psKernel *convolution // Convolution, conv(x) - ) -{ - // Range of convolution - int xMin = convolution->xMin; - int xMax = convolution->xMax; - int yMin = convolution->yMin; - int yMax = convolution->yMax; - + int kernelIndex, // Index for kernel component + int footprint // Size of region of interest + ) +{ psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - - for (int y = yMin; y <= yMax; y++) { - psF32 *in = &input->kernel[y][xMin]; // Dereference input - psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight - psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution - for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) { + psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { double temp = *in / *wt; // Temporary product *sumI += temp; @@ -82,19 +81,15 @@ double *sumCC, // Sum of conv(x)^2/sigma(x)^2 const pmSubtractionStamp *stamp, // Stamp with input and weight - const psKernel *convolution // Convolution, conv(x) - ) -{ - // Range of convolution - int xMin = convolution->xMin; - int xMax = convolution->xMax; - int yMin = convolution->yMin; - int yMax = convolution->yMax; - + int kernelIndex, // Index for kernel component + int footprint // Size of region of interest + ) +{ psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - - for (int y = yMin; y <= yMax; y++) { - psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight - psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution - for (int x = xMin; x <= xMax; x++, wt++, conv++) { + psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + + for (int y = -footprint; y <= footprint; y++) { + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, wt++, conv++) { double convNoise = *conv / *wt; // Temporary product *sumC += convNoise; @@ -106,24 +101,20 @@ static double accumulateChi2(psKernel *input, // Input stamp - psKernel *convolution, // Convolution, R(x)*k(u) pmSubtractionStamp *stamp, // Stamp with weight + int kernelIndex, // Index for kernel component double coeff, // Coefficient of convolution - double bg // Background term - ) -{ - // Range of convolution - int xMin = convolution->xMin; - int xMax = convolution->xMax; - int yMin = convolution->yMin; - int yMax = convolution->yMax; - + double bg, // Background term + int footprint // Size of region of interest + ) +{ double chi2 = 0.0; - psKernel *weight = stamp->weight; - - for (int y = yMin; y <= yMax; y++) { - psF32 *in = &input->kernel[y][xMin]; // Dereference input - psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight - psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution - for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) { + psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt; } @@ -159,22 +150,16 @@ static void subtractConvolution(psKernel *input, // Input stamp const pmSubtractionStamp *stamp, // Stamp with weight - const psKernel *convolution, // Convolution, R(x)*k(u) + int kernelIndex, // Index for kernel component float coeff, // Coefficient of subtraction - float bg // Background term - ) -{ - // Range of convolution - int xMin = convolution->xMin; - int xMax = convolution->xMax; - int yMin = convolution->yMin; - int yMax = convolution->yMax; - - psKernel *weight = stamp->weight; // Weight map - - for (int y = yMin; y <= yMax; y++) { - psF32 *in = &input->kernel[y][xMin]; // Dereference input - psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution - psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight - for (int x = xMin; x <= xMax; x++, in++, conv++, wt++) { + float bg, // Background term + int footprint // Size of region of interest + ) +{ + psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, in++, conv++) { *in -= *conv * coeff + bg; } @@ -229,18 +214,18 @@ psKernel *input = stamp->input; // Input image of interest psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image - inputs->data[i] = psKernelAllocFromImage(copy, size + footprint + 1, size + footprint + 1); + inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); psFree(copy); // Drop reference } // Generate the convolutions - psArray *convolutions = psArrayAlloc(numKernels); // Convolutions per kernel component and stamp - for (int i = 0; i < numKernels; i++) { - psKernel *kernel = kernels->preCalc->data[i]; // Precalculated kernel - - psArray *stampConv = psArrayAlloc(numStamps); // Convolutions for each stamp - convolutions->data[i] = stampConv; - - for (int j = 0; j < numStamps; j++) { - stampConv->data[j] = convolveStamp(stamps->data[j], kernel, footprint); + for (int i = 0; i < numStamps; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + if (!stamp->convolutions) { + stamp->convolutions = psArrayAlloc(numKernels); + } + + for (int j = 0; j < numKernels; j++) { + psKernel *kernel = kernels->preCalc->data[j]; // Precalculated kernel + stamp->convolutions->data[j] = p_pmSubtractionConvolveStampPrecalc(stamp->reference, kernel); } } @@ -264,8 +249,6 @@ psVectorInit(sumCC, 0.0); for (int i = 0; i < numKernels; i++) { - psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp for (int j = 0; j < numStamps; j++) { - accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j], - stampsConv->data[j]); + accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j], i, footprint); } } @@ -298,8 +281,6 @@ double sumIC = 0.0; // sum of I(x)C(x)/sigma(x)^2 - psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp - for (int j = 0; j < numStamps; j++) { - accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], stampsConv->data[j]); + accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], i, footprint); } @@ -310,5 +291,5 @@ double chi2 = 0.0; // Chi^2 for (int j = 0; j < numStamps; j++) { - chi2 += accumulateChi2(inputs->data[j], stampsConv->data[j], stamps->data[j], coeff, bg); + chi2 += accumulateChi2(inputs->data[j], stamps->data[j], i, coeff, bg, footprint); } @@ -329,8 +310,6 @@ ranking->data.S32[bestIndex] = iter; // Remove its contribution, and don't include it in the future. - psArray *stampsConv = convolutions->data[bestIndex]; // Convolutions per stamp for (int j = 0; j < numStamps; j++) { - subtractConvolution(inputs->data[j], stamps->data[j], stampsConv->data[j], - bestCoeff, bestBG); + subtractConvolution(inputs->data[j], stamps->data[j], bestIndex, bestCoeff, bestBG, footprint); } @@ -349,5 +328,4 @@ } psFree(inputs); - psFree(convolutions); psFree(sumC); psFree(sumCC); @@ -392,6 +370,8 @@ psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest for (int i = 1; i < newSize; i++) { - psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest - psBinaryOp(kernel->image, kernel->image, "-", subtract->image); + 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); + } } } else if (type == PM_SUBTRACTION_KERNEL_GUNK) { @@ -400,6 +380,8 @@ for (int i = 0; i < newSize; i++) { - psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest - kernel->kernel[0][0] -= 1.0; + 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; + } }