Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 14648) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 14701) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-23 23:43:12 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-30 03:50:28 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -17,4 +17,5 @@ #include #include +#include #include @@ -24,4 +25,6 @@ #include "pmSubtraction.h" + +//#define TESTING #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time @@ -216,82 +219,116 @@ } -// Generate the convolved pixel value -static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions - int index, // Kernel basis function index - int x, int y, // Pixel around which to convolve - const psKernel *image // Image to convolve (a kernel for convenience) - ) +// Subtract the (0,0) element to preserve photometric scaling +static void convolveSub(psKernel *convolved, // Convolved image + const psKernel *image, // Image being convolved + int footprint // Size of region of interest + ) +{ + // Can't use psBinaryOp because the images are of different size + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + convolved->kernel[y][x] -= image->kernel[y][x]; + } + } + return; +} + +// Generate the convolution given some offset +static psKernel *convolveOffset(const psKernel *image, // Image to convolve (a kernel for convenience) + int u, int v, // Offset to apply + int footprint // Size of region of interest + ) +{ + psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image + int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy + for (int y = -footprint; y <= footprint; y++) { + // Convolution with a delta function is just the value specified by the offset + memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes); + } + return convolved; +} + +// 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 + int index, // Kernel basis function index + const psKernel *image, // Image to convolve (a kernel for convenience) + int footprint // Size of region of interest + ) { switch (kernels->type) { case PM_SUBTRACTION_KERNEL_POIS: { - // Convolution with a delta function is just the value specified by the offset int u = kernels->u->data.S32[index]; // Offset in x int v = kernels->v->data.S32[index]; // Offset in y - float value = image->kernel[y + v][x + u]; // Value of convolution + psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - // The (0,0) element is subtracted from most kernels to preserve photometric scaling - value -= image->kernel[y][x]; + convolveSub(convolved, image, footprint); } - return value; + return convolved; } // Method for SPAM and FRIES is the same case PM_SUBTRACTION_KERNEL_SPAM: case PM_SUBTRACTION_KERNEL_FRIES: { + psKernel *convolved = psKernelAlloc(-footprint, footprint, + -footprint, footprint); // Convolved image int uStart = kernels->u->data.S32[index]; int uStop = kernels->uStop->data.S32[index]; int vStart = kernels->v->data.S32[index]; int vStop = kernels->vStop->data.S32[index]; - double sum = 0.0; - for (int v = vStart; v <= vStop; v++) { - for (int u = uStart; u <= uStop; u++) { - sum += image->kernel[y + v][x + u]; + float norm = 1.0 / (uStop - uStart + 1) * (vStop - vStart + 1); // Normalisation + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + double sum = 0.0; + for (int v = vStart; v <= vStop; v++) { + for (int u = uStart; u <= uStop; u++) { + sum += image->kernel[y - v][x - u]; + } + } + convolved->kernel[y][x] = norm * sum; } } - sum /= (uStop - uStart + 1) * (vStop - vStart + 1); // Normalising sum of kernel component to unity if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - // The (0,0) element is subtracted from most kernels to preserve photometric scaling - sum -= image->kernel[y][x]; + convolveSub(convolved, image, footprint); } - return sum; + return convolved; } case PM_SUBTRACTION_KERNEL_GUNK: { if (index < kernels->inner) { - // Using pre-calculated function - psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel - int size = kernels->size; // Kernel half-size - double sum = 0.0; // Accumulated sum from convolution - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - sum += kernel->kernel[v][u] * image->kernel[y + v][x + u]; - } - } - return sum; + // Photometric scaling is already built in to the precalculated kernel + return convolvePrecalc(image, kernels->preCalc->data[index], footprint); } // Using delta function int u = kernels->u->data.S32[index]; // Offset in x int v = kernels->v->data.S32[index]; // Offset in y - float value = image->kernel[y + v][x + u]; // Value of convolution - // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling + psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - value -= image->kernel[y][x]; + convolveSub(convolved, image, footprint); } - return value; + return convolved; } case PM_SUBTRACTION_KERNEL_ISIS: { - psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel - int size = kernels->size; // Kernel half-size - double sum = 0.0; // Accumulated sum from convolution - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - sum += kernel->kernel[v][u] * image->kernel[y + v][x + u]; - // Photometric scaling is already built in to the precalculated kernel - } - } - return sum; + // Photometric scaling is already built in to the precalculated kernel + return convolvePrecalc(image, kernels->preCalc->data[index], footprint); } case PM_SUBTRACTION_KERNEL_RINGS: { + psKernel *convolved = psKernelAlloc(-footprint, footprint, + -footprint, footprint); // Convolved image if (index == kernels->subIndex) { - return image->kernel[y][x]; + // Simply copying over the image + return convolveOffset(image, 0, 0, footprint); } + psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data psVector *uCoords = preCalc->data[0]; // u coordinates @@ -299,16 +336,21 @@ psVector *poly = preCalc->data[2]; // Polynomial values int num = uCoords->n; // Number of pixels - double sum = 0.0; // Accumulated sum from convolution - for (int j = 0; j < num; j++) { - int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates - sum += image->kernel[y + v][x + u] * poly->data.F32[j]; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + double sum = 0.0; // Accumulated sum from convolution + for (int j = 0; j < num; j++) { + int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates + sum += image->kernel[y - v][x - u] * poly->data.F32[j]; + } + convolved->kernel[y][x] = sum; + // Photometric scaling is built into the kernel --- no subtraction! + } } - // Photometric scaling is built into the kernel --- no subtraction! - return sum; + return convolved; } default: psAbort("Should never get here."); } - return NAN; + return NULL; } @@ -460,5 +502,4 @@ int numParams = numKernels * numSpatial + numBackground; int bgIndex = numParams - numBackground; // Index in matrix for the background - psVector *convolutions = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); // Convolutions // We iterate over each stamp, allocate the matrix and vectors if @@ -466,106 +507,165 @@ for (int i = 0; i < stamps->n; i++) { pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_CALCULATE) { - psImage *stampMatrix = stamp->matrix; // Least-squares matrix for this stamp - psVector *stampVector = stamp->vector; // Least-squares vector for this stamp - - if (!stampMatrix) { - stamp->matrix = stampMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); - } - assert(stampMatrix->type.type == PS_TYPE_F64); - assert(stampMatrix->numCols == numParams && stampMatrix->numRows == numParams); - psImageInit(stampMatrix, 0.0); - - if (!stampVector) { - stamp->vector = stampVector = psVectorAlloc(numParams, PS_TYPE_F64); - } - assert(stampVector->type.type == PS_TYPE_F64); - assert(stampVector->n == numParams); - psVectorInit(stampVector, 0.0); - - // Spatial polynomial terms - psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); - - psKernel *reference = stamp->reference; // Reference postage stamp - psKernel *input = stamp->input; // Input postage stamp - psKernel *weight = stamp->weight; // Weight map postage stamp - - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - float invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse square noise - - // Generate the convolutions - for (int j = 0; j < numKernels; j++) { - double value = convolvePixel(kernels, j, x, y, reference); // Value from convolution - // Generate the pseudo-convolutions from the spatial polynomial terms - for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) { - for (int xOrder = 0; xOrder <= spatialOrder - yOrder; - xOrder++, index += numKernels) { - convolutions->data.F64[index] = value * polyValues->data.F64[yOrder][xOrder]; + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { + continue; + } + + psImage *stampMatrix = stamp->matrix; // Least-squares matrix for this stamp + psVector *stampVector = stamp->vector; // Least-squares vector for this stamp + + if (!stampMatrix) { + stamp->matrix = stampMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + } + assert(stampMatrix->type.type == PS_TYPE_F64); + assert(stampMatrix->numCols == numParams && stampMatrix->numRows == numParams); + psImageInit(stampMatrix, 0.0); + + if (!stampVector) { + stamp->vector = stampVector = psVectorAlloc(numParams, PS_TYPE_F64); + } + assert(stampVector->type.type == PS_TYPE_F64); + assert(stampVector->n == numParams); + psVectorInit(stampVector, 0.0); + + psKernel *reference = stamp->reference; // Reference postage stamp + psKernel *input = stamp->input; // Input postage stamp + psKernel *weight = stamp->weight; // Weight map postage stamp + + // Generate convolutions of the reference + psArray *convolutions = stamp->convolutions; // Convolutions of the reference for each kernel + if (!convolutions) { + stamp->convolutions = convolutions = psArrayAlloc(numKernels); + } + for (int j = 0; j < numKernels; j++) { + if (convolutions->data[j]) { + psFree(convolutions->data[j]); + } + convolutions->data[j] = convolveRef(kernels, j, reference, footprint); +#ifdef TESTING + { + psKernel *conv = convolutions->data[j]; // Convolution of interest + psString filename = NULL; + psStringAppend(&filename, "conv_%03d_%03d.fits", i, j); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, conv->image, 0, NULL); + psFitsClose(fits); + } +#endif + } + + psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms + + // Generate least-squares vector and matrix + for (int j = 0; j < numKernels; j++) { + psKernel *jConv = convolutions->data[j]; // Convolution for j-th element + + // Generate upper diagonals of matrix + for (int k = 0; k < numKernels; k++) { + psKernel *kConv = convolutions->data[k]; // Convolution for k-th element + double sumCC = 0.0; // Sum of the convolution products + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sumCC += jConv->kernel[y][x] * kConv->kernel[y][x] / weight->kernel[y][x]; + } + } + // Generate the pseudo-convolutions from the spatial polynomial terms + for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; + jxOrder++, jIndex += numKernels) { + for (int kyOrder = 0, kIndex = k; kyOrder <= spatialOrder; kyOrder++) { + for (int kxOrder = 0; kxOrder <= spatialOrder - kyOrder; + kxOrder++, kIndex += numKernels) { + stampMatrix->data.F64[jIndex][kIndex] = sumCC * + polyValues->data.F64[jyOrder][jxOrder] * + polyValues->data.F64[kyOrder][kxOrder]; } } } - - // Generate the least-squares matrix and vector - // Upper diagonal only - for (int i = 0; i < bgIndex; i++) { - for (int j = i; j < bgIndex; j++) { - stampMatrix->data.F64[i][j] += convolutions->data.F64[i] * - convolutions->data.F64[j] * invNoise2; + } + } + + // Vector and background term for matrix + double sumC = 0.0; // Sum of the convolution + double sumIC = 0.0; // Sum of the convolution/input product + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double convProduct = jConv->kernel[y][x] / weight->kernel[y][x]; // Convolution / noise^2 + sumC += convProduct; + sumIC += input->kernel[y][x] * convProduct; + } + } + for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; + jxOrder++, jIndex += numKernels) { + stampVector->data.F64[jIndex] = sumIC * polyValues->data.F64[jyOrder][jxOrder]; + stampMatrix->data.F64[jIndex][bgIndex] = sumC * polyValues->data.F64[jyOrder][jxOrder]; + stampMatrix->data.F64[bgIndex][jIndex] = sumC * polyValues->data.F64[jyOrder][jxOrder]; + } + } + + } + psFree(polyValues); + + // Background only terms + double sum1 = 0.0; // Sum of the weighting + double sumI = 0.0; // Sum of the input + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse noise, squared + sum1 += invNoise2; + sumI += input->kernel[y][x] * invNoise2; + } + } + stampMatrix->data.F64[bgIndex][bgIndex] = sum1; + stampVector->data.F64[bgIndex] = sumI; + + // Fill in lower diagonals of symmetric matrix, while checking for bad values + // Note, there are two symmetries going on here --- the matrix is C_i C_j P_i P_j + bool bad = false; // Are there bad values? +#if 0 + for (int j = 0; j < numKernels; j++) { + for (int jSpatial = 0, jIndex = j; jSpatial < numSpatial; jSpatial++, jIndex += numKernels) { + for (int k = 0; k < j; k++) { + for (int kSpatial = 0, kIndex = k; kSpatial < numSpatial; + kSpatial++, kIndex += numKernels) { + double value = stampMatrix->data.F64[kIndex][jIndex]; // Value of matrix + stampMatrix->data.F64[jIndex][kIndex] = value; + if (!isfinite(value)) { + bad = true; } - stampVector->data.F64[i] += input->kernel[y][x] * convolutions->data.F64[i] * - invNoise2; - - // Background term - stampMatrix->data.F64[i][bgIndex] += convolutions->data.F64[i] * invNoise2; } - // Background only terms - stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; - stampVector->data.F64[bgIndex] += input->kernel[y][x] * invNoise2; - } - } - psFree(polyValues); - - // Fill in lower diagonal of symmetric matrix, while checking for bad values - bool bad = false; // Are there bad values? - for (int i = 0; i < bgIndex; i++) { - for (int j = 0; j < i; j++) { - stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i]; - if (!isfinite(stampMatrix->data.F64[j][i])) { - bad = true; - } - } - stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex]; - if (!isfinite(stampMatrix->data.F64[i][bgIndex]) || - !isfinite(stampMatrix->data.F64[i][i]) || - !isfinite(stampVector->data.F64[i])) { + } + + stampMatrix->data.F64[bgIndex][jIndex] = stampMatrix->data.F64[jIndex][bgIndex]; + if (!isfinite(stampMatrix->data.F64[jIndex][bgIndex]) || + !isfinite(stampMatrix->data.F64[jIndex][jIndex]) || + !isfinite(stampVector->data.F64[jIndex])) { bad = true; } } - if (!isfinite(stampVector->data.F64[bgIndex])) { - bad = true; - } - - if (bad) { - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n", - i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); - } else { - stamp->status = PM_SUBTRACTION_STAMP_USED; - } - - if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { - psString matrixName = NULL; - psStringAppend(&matrixName, "matrix%d.fits", i); - psFits *matrixFile = psFitsOpen(matrixName, "w"); - psFree(matrixName); - psFitsWriteImage(matrixFile, NULL, stampMatrix, 0, NULL); - psFitsClose(matrixFile); - } - - } - } - - psFree(convolutions); + } +#endif + if (!isfinite(stampVector->data.F64[bgIndex])) { + bad = true; + } + + if (bad) { + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + } else { + stamp->status = PM_SUBTRACTION_STAMP_USED; + } + + if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { + psString matrixName = NULL; + psStringAppend(&matrixName, "matrix%d.fits", i); + psFits *matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stampMatrix, 0, NULL); + psFitsClose(matrixFile); + } + } return true; @@ -665,10 +765,9 @@ double totalSquareDev = 0.0; // Total square deviation from zero int numStamps = 0; // Number of used stamps + int numKernels = kernels->num; // Number of kernels + int spatialOrder = kernels->spatialOrder; // Order of kernel spatial variations + double devNorm = 1.0 / PS_SQR(2 * footprint + 1); // Normalisation for deviations { - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics - psKernel *convolution = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolution - psKernel *kernelImage = NULL; // The kernel, with which to convolve the stamps float background = solution->data.F64[solution->n-1]; // The difference in background - int size = kernels->size; // Half-size of the kernel for (int i = 0; i < stamps->n; i++) { @@ -681,60 +780,54 @@ stamp->yNorm); // Polynomial terms - psKernel *reference = stamp->reference; // Reference postage stamp - - kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false); + psKernel *input = stamp->input; // Input postage stamp + psKernel *weight = stamp->weight; // Weight postage stamp + psArray *convolutions = stamp->convolutions; // Convolutions of reference image for each kernel +#ifdef TESTING + psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); +#endif + float deviation = 0.0; // Deviation for this stamp for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - convolution->kernel[y][x] = background; - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - convolution->kernel[y][x] += kernelImage->kernel[v][u] * - reference->kernel[y + v][x + u]; + float conv = background; // The value of the convolution + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution of reference + double polynomial = 0.0; // Value of the polynomial + for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; + xOrder++, index += numKernels) { + polynomial += solution->data.F64[index] * + polyValues->data.F64[yOrder][xOrder]; + } } + conv += convolution->kernel[y][x] * polynomial; } + float diff = input->kernel[y][x] - conv; + deviation += PS_SQR(diff) / weight->kernel[y][x]; +#ifdef TESTING + residual->kernel[y][x] = diff; +#endif } } psFree(polyValues); - psImage *convolvedStamp = convolution->image; // Image of the convolution - psImage *input = stamp->input->image; // Input image postage stamp - psImage *weight = stamp->weight->image; // Weight image postage stamp - - // Region of interest - psRegion inRegion = psRegionSet(input->col0 + size, input->col0 + size + 2 * footprint + 1, - input->row0 + size, input->row0 + size + 2 * footprint + 1); - psRegion wtRegion = (input == weight ? inRegion : - psRegionSet(weight->col0 + size, weight->col0 + size + 2 * footprint + 1, - weight->row0 + size, weight->row0 + size + 2 * footprint + 1)); - - psImage *inStamp = psImageSubset(stamp->input->image, inRegion); // Image of stamp - psImage *wtStamp = psImageSubset(stamp->weight->image, wtRegion); // Image of stamp - assert(convolvedStamp->numCols == inStamp->numCols && - convolvedStamp->numRows == inStamp->numRows); - assert(convolvedStamp->numCols == wtStamp->numCols && - convolvedStamp->numRows == wtStamp->numRows); - (void)psBinaryOp(convolvedStamp, inStamp, "-", convolvedStamp); - (void)psBinaryOp(convolvedStamp, convolvedStamp, "*", convolvedStamp); - (void)psBinaryOp(convolvedStamp, convolvedStamp, "/", wtStamp); - psFree(inStamp); - psFree(wtStamp); - if (!psImageStats(stats, convolvedStamp, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, - "Unable to calculate statistics on normalised residual image of stamp."); - psFree(deviations); - psFree(convolution); - psFree(stats); - return -1; - } - deviations->data.F32[i] = sqrt(stats->sampleMean / 2.0); - psTrace("psModules.imcombine", 1, "Deviation for stamp %d (%d,%d): %f\n", - i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); + deviations->data.F32[i] = devNorm * deviation; + psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", + i, (int)stamp->x, (int)stamp->y, deviations->data.F32[i]); totalSquareDev += PS_SQR(deviations->data.F32[i]); numStamps++; - } - - psFree(kernelImage); - psFree(convolution); - psFree(stats); + +#ifdef TESTING + { + psString filename = NULL; + psStringAppend(&filename, "resid_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, residual->image, 0, NULL); + psFitsClose(fits); + } + psFree(residual); +#endif + + } } @@ -823,9 +916,66 @@ } +// 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, psMaskType blank, const psRegion *region, const psVector *solution, - const pmSubtractionKernels *kernels) + const pmSubtractionKernels *kernels, bool useFFT) { PS_ASSERT_IMAGE_NON_NULL(inImage, false); @@ -927,5 +1077,7 @@ for (int j = yMin; j < yMax; j += fullSize) { + int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest for (int i = xMin; i < xMax; i += fullSize) { + int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest // Only generate polynomial values every kernel footprint, since we have already assumed @@ -935,4 +1087,5 @@ 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols, 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows); + kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false); if (inWeight) { @@ -940,7 +1093,26 @@ } - for (int y = j; y < PS_MIN(j + fullSize, yMax); y++) { - for (int x = i; x < PS_MIN(i + fullSize, xMax); x++) { - // Check and propagate the kernel footprint, if required + + if (useFFT) { + // Use Fast Fourier Transform to do the convolution + // This provides a big speed-up for large kernels + convolveFFT(convImage, inImage, kernelImage, psRegionSet(i, xSubMax, j, ySubMax), + background, size); + if (inWeight) { + convolveFFT(convWeight, inWeight, kernelWeight, psRegionSet(i, xSubMax, j, ySubMax), + 0.0, size); + } + } else { + convolveDirect(convImage, inImage, kernelImage, psRegionSet(i, xSubMax, j, ySubMax), + background, size); + if (inWeight) { + convolveDirect(convWeight, inWeight, kernelWeight, psRegionSet(i, xSubMax, j, ySubMax), + 0.0, size); + } + } + + // Propagate the mask + for (int y = j; y < ySubMax; y++) { + for (int x = i; x < xSubMax; x++) { if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) { @@ -950,26 +1122,8 @@ convWeight->data.F32[y][x] = NAN; } - } else { - // Convolve the image - convImage->data.F32[y][x] = background; - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - convImage->data.F32[y][x] += kernelImage->kernel[v][u] * - inImage->data.F32[y + v][x + u]; - } - } - - // Convolve the weight (variance) map - if (inWeight) { - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - convWeight->data.F32[y][x] += kernelWeight->kernel[v][u] * - inWeight->data.F32[y + v][x + u]; - } - } - } } } } + } }