Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 15427) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 15443) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.71 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-01 01:11:03 $ + * @version $Revision: 1.72 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-03 02:28:24 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -76,5 +76,6 @@ // Generate the kernel to apply to the variance from the normal kernel -static psKernel *varianceKernel(psKernel *normalKernel // Normal kernel +static psKernel *varianceKernel(psKernel *out, // Output kernel + psKernel *normalKernel // Normal kernel ) { @@ -83,5 +84,7 @@ int yMin = normalKernel->yMin, yMax = normalKernel->yMax; - psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax); // Kernel to return + if (!out) { + out = psKernelAlloc(xMin, xMax, yMin, yMax); + } // Take the square of the normal kernel @@ -93,5 +96,5 @@ sumNormal += value; sumVariance += value2; - kernel->kernel[v][u] = value2; + out->kernel[v][u] = value2; } } @@ -99,8 +102,7 @@ // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel // This is required to keep the relative scaling between the image and the weight map - psBinaryOp(kernel->image, kernel->image, "*", - psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32)); - - return kernel; + psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32)); + + return out; } @@ -264,10 +266,275 @@ } -// 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 - ) + +// Convolve an image using FFT +static void convolveFFT(psImage *target,// Place the result in here + const psImage *image, // Image to convolve + const psImage *mask, // Mask image + psMaskType maskVal, // Value to mask + 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 *subMask = mask ? psImageSubset((psImage*)mask, border) : NULL; // Subimage mask + psImage *convolved = psImageConvolveFFT(subImage, subMask, maskVal, kernel, 0.0); // Convolution + psFree(subImage); + psFree(subMask); + + // 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 + psImage *mask, // Mask to mark as blank (or NULL) + psImage *weight, // Weight map to mark as blank (or NULL) + int x, int y, // Coordinates to mark blank + psMaskType blank // Blank mask value + ) +{ + image->data.F32[y][x] = NAN; + if (mask) { + mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + } + if (weight) { + weight->data.F32[y][x] = NAN; + } + 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, NULL, 0, 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; +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psMaskType maskVal, + int size, int footprint, float badFrac, bool useFFT) +{ + PS_ASSERT_IMAGE_NON_NULL(mask1, NULL); + PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_MASK, NULL); + if (mask2) { + PS_ASSERT_IMAGE_NON_NULL(mask2, NULL); + PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_MASK, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL); + } + PS_ASSERT_INT_NONNEGATIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); + if (isfinite(badFrac)) { + PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL); + } + + int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images + + // Dereference inputs for convenience + psMaskType **data1 = mask1->data.PS_TYPE_MASK_DATA; + psMaskType **data2 = NULL; + if (mask2) { + data2 = mask2->data.PS_TYPE_MASK_DATA; + } + + // First, a pass through to determine the fraction of bad pixels + if (isfinite(badFrac) && badFrac != 1.0) { + int numBad = 0; // Number of bad pixels + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (data1[y][x] & maskVal) { + numBad++; + continue; + } + if (data2 && data2[y][x] & maskVal) { + numBad++; + } + } + } + if (numBad > badFrac * numCols * numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", + numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac); + return NULL; + } + } + + // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask + psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask + psImageInit(mask, 0); + psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience + + // Block out a border around the edge of the image + + // Bottom stripe + for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Either side + for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { + for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Top stripe + for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + + // XXX Could do something smarter here --- we will get images that are predominantly masked (where the + // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad + // pixel is wasting time. As a first cut, I've put in a check on the fraction of bad pixels, but we could + // imagine looking for the edge of big regions and convolving just at the edge. As a second cut, allow + // use of FFT convolution. + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (data1[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1; + } + if (data2 && data2[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2; + } + } + } + + // Block out the entire stamp footprint around bad input pixels. + + // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad + // reference pixel (within 'size'). Then we want to block out with the FOOTPRINT mask everything within a + // footprint's distance of those (within 'footprint'). + + if (useFFT) { + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, + -size, size, -size, size, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, + -size, size, -size, size, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, + PM_SUBTRACTION_MASK_FOOTPRINT_1, + -footprint, footprint, -footprint, footprint, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, + PM_SUBTRACTION_MASK_FOOTPRINT_2, + -footprint, footprint, -footprint, footprint, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + } else { + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_1, + PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, + -size, size, -size, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, + -size, size, -size, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, + PM_SUBTRACTION_MASK_FOOTPRINT_1, + -footprint, footprint, -footprint, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, + PM_SUBTRACTION_MASK_FOOTPRINT_2, + -footprint, footprint, -footprint, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + } + + return mask; +} + +// Convolve a stamp by a single kernel basis function +static psKernel *convolveStampSingle(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) { @@ -326,11 +593,10 @@ } case PM_SUBTRACTION_KERNEL_RINGS: { - psKernel *convolved = psKernelAlloc(-footprint, footprint, - -footprint, footprint); // Convolved image if (index == kernels->subIndex) { // Simply copying over the image return convolveOffset(image, 0, 0, footprint); } - + psKernel *convolved = psKernelAlloc(-footprint, footprint, + -footprint, footprint); // Convolved image psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data psVector *uCoords = preCalc->data[0]; // u coordinates @@ -357,255 +623,35 @@ } -// Convolve an image using FFT -static void convolveFFT(psImage *target,// Place the result in here - const psImage *image, // Image to convolve - const psImage *mask, // Mask image - 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 *subMask = mask ? psImageSubset((psImage*)mask, border) : NULL; // Subimage mask - psImage *convolved = psImageConvolveFFT(subImage, subMask, PM_SUBTRACTION_MASK_REF, - kernel, 0.0); // Convolution - psFree(subImage); - psFree(subMask); - - // 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 - psImage *mask, // Mask to mark as blank (or NULL) - psImage *weight, // Weight map to mark as blank (or NULL) - int x, int y, // Coordinates to mark blank - psMaskType blank // Blank mask value +// Convolve the stamp by each of the kernel basis functions +static psArray *convolveStamp(psArray *convolutions, // The convolutions + const psKernel *image, // Image to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int footprint // Stamp half-size ) { - image->data.F32[y][x] = NAN; - if (mask) { - mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; - } - if (weight) { - weight->data.F32[y][x] = NAN; - } - 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, NULL, 0, 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; -} - - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Public functions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -psImage *pmSubtractionMask(const psImage *refMask, const psImage *inMask, psMaskType maskVal, - int size, int footprint, float badFrac, bool useFFT) -{ - PS_ASSERT_IMAGE_NON_NULL(refMask, NULL); - PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL); - if (inMask) { - PS_ASSERT_IMAGE_NON_NULL(inMask, NULL); - PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL); - } - PS_ASSERT_INT_NONNEGATIVE(size, NULL); - PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); - if (isfinite(badFrac)) { - PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL); - PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL); - } - - // Size of the images - int numCols = refMask->numCols; - int numRows = refMask->numRows; - - // Dereference inputs for convenience - psMaskType **inData = NULL; - if (inMask) { - inData = inMask->data.PS_TYPE_MASK_DATA; - } - psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA; - - // First, a pass through to determine the fraction of bad pixels - if (isfinite(badFrac) && badFrac != 1.0) { - int numBad = 0; // Number of bad pixels - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - if (inData && inData[y][x] & maskVal) { - numBad++; - continue; - } - if (refData[y][x] & maskVal) { - numBad++; - } - } - } - if (numBad > badFrac * numCols * numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", - numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac); - return NULL; - } - } - - // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask - psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask - psImageInit(mask, 0); - psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience - - // Block out a border around the edge of the image - - // Bottom stripe - for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { - for (int x = 0; x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - // Either side - for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { - for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - // Top stripe - for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - - // XXX Could do something smarter here --- we will get images that are predominantly masked (where the - // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad - // pixel is wasting time. As a first cut, I've put in a check on the fraction of bad pixels, but we could - // imagine looking for the edge of big regions and convolving just at the edge. As a second cut, allow - // use of FFT convolution. - - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - if (inData && inData[y][x] & maskVal) { - maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT; - } - if (refData[y][x] & maskVal) { - maskData[y][x] |= PM_SUBTRACTION_MASK_REF; - } - } - } - - // Block out the entire stamp footprint around bad input pixels. - - // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad - // reference pixel (within 'size'). Then we want to block out with the FOOTPRINT mask everything within a - // footprint's distance of those (within 'footprint'). - - if (useFFT) { - if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_INPUT, PM_SUBTRACTION_MASK_FOOTPRINT, - -footprint, footprint, -footprint, footprint, 0.5)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad input pixels."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_REF, PM_SUBTRACTION_MASK_CONVOLVE, - -size, size, -size, size, 0.5)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad reference pixels."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE, PM_SUBTRACTION_MASK_FOOTPRINT, - -footprint, footprint, -footprint, footprint, 0.5)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad reference pixels."); - psFree(mask); - return NULL; - } - } else { - if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_INPUT, PM_SUBTRACTION_MASK_FOOTPRINT, - -footprint, footprint, -footprint, footprint)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad input pixels."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_REF, PM_SUBTRACTION_MASK_CONVOLVE, - -size, size, -size, size)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad reference pixels."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE, - PM_SUBTRACTION_MASK_FOOTPRINT, - -footprint, footprint, -footprint, footprint)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad reference pixels."); - psFree(mask); - return NULL; - } - } - - return mask; -} - -bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint) + assert(image); + assert(kernels); + assert(footprint >= 0); + + if (convolutions) { + // Already done + return convolutions; + } + + int numKernels = kernels->num; // Number of kernels + convolutions = psArrayAlloc(numKernels); + + for (int i = 0; i < numKernels; i++) { + convolutions->data[i] = convolveStampSingle(kernels, i, image, footprint); + } + + return convolutions; +} + + +bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint, + pmSubtractionMode mode) { PS_ASSERT_PTR_NON_NULL(stamp, false); - PS_ASSERT_KERNEL_NON_NULL(stamp->reference, false); PS_ASSERT_PTR_NON_NULL(kernels, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); @@ -617,14 +663,19 @@ } - if (stamp->convolutions) { - // Already done - return true; - } - - stamp->convolutions = psArrayAlloc(kernels->num); - psKernel *reference = stamp->reference; // Reference postage stamp - - for (int i = 0; i < kernels->num; i++) { - stamp->convolutions->data[i] = convolveRef(kernels, i, reference, footprint); + switch (mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); + break; + case PM_SUBTRACTION_MODE_2: + stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint); + break; + case PM_SUBTRACTION_MODE_UNSURE: + case PM_SUBTRACTION_MODE_DUAL: + stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); + stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint); + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); } @@ -634,5 +685,5 @@ bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, - int footprint) + int footprint, pmSubtractionMode mode) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -676,14 +727,28 @@ psVectorInit(stampVector, 0.0); - psKernel *input = stamp->input; // Input postage stamp - psKernel *weight = stamp->weight; // Weight map postage stamp - - // Generate convolutions of the reference - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + // Generate convolutions + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); return NULL; } - psArray *convolutions = stamp->convolutions; // Convolutions of the reference for each kernel + psKernel *weight = stamp->weight; // Weight map postage stamp + psKernel *target; // Target postage stamp + psArray *convolutions; // Convolutions for each kernel + + switch (mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + convolutions = stamp->convolutions1; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + convolutions = stamp->convolutions2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms @@ -739,10 +804,10 @@ // Vector and background term for matrix double sumC = 0.0; // Sum of the convolution - double sumIC = 0.0; // Sum of the convolution/input product + double sumTC = 0.0; // Sum of the convolution/target 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; + sumTC += target->kernel[y][x] * convProduct; } } @@ -752,5 +817,5 @@ continue; } - if (!isfinite(sumIC)) { + if (!isfinite(sumTC)) { psStringAppend(&badString, "\nBad sumIC detected at %d", j); bad = true; @@ -761,5 +826,5 @@ for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex += numKernels) { - stampVector->data.F64[jIndex] = sumIC * polyValues->data.F64[jyOrder][jxOrder]; + stampVector->data.F64[jIndex] = sumTC * polyValues->data.F64[jyOrder][jxOrder]; double convPoly = sumC * polyValues->data.F64[jyOrder][jxOrder]; // Value @@ -774,10 +839,10 @@ // Background only terms double sum1 = 0.0; // Sum of the weighting - double sumI = 0.0; // Sum of the input + double sumT = 0.0; // Sum of the target 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; + sumT += target->kernel[y][x] * invNoise2; } } @@ -785,5 +850,5 @@ psStringAppend(&badString, "\nBad sum1 detected"); bad = true; - } else if (!isfinite(sumI)) { + } else if (!isfinite(sumT)) { psStringAppend(&badString, "\nBad sumI detected"); bad = true; @@ -791,5 +856,5 @@ stampMatrix->data.F64[bgIndex][bgIndex] = sum1; - stampVector->data.F64[bgIndex] = sumI; + stampVector->data.F64[bgIndex] = sumT; if (bad) { @@ -890,119 +955,150 @@ } - -int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, psImage *subMask, const psVector *solution, - int footprint, float sigmaRej, const pmSubtractionKernels *kernels) +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, const psVector *solution, + int footprint, const pmSubtractionKernels *kernels, + pmSubtractionMode mode) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); + PS_ASSERT_VECTOR_NON_NULL(solution, NULL); + PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); + PS_ASSERT_PTR_NON_NULL(kernels, NULL); + PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + + polyTerms(kernels->bgOrder), NULL); + PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); + + psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps + double devNorm = 1.0 / PS_SQR(2 * footprint + 1); // Normalisation for deviations + int numKernels = kernels->num; // Number of kernels + int spatialOrder = kernels->spatialOrder; // Order of kernel spatial variations + float background = solution->data.F64[solution->n-1]; // The difference in background + + psVector *coeff = psVectorAlloc(numKernels, PS_TYPE_F64); // Coefficients for the kernel basis functions + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + deviations->data.F32[i] = NAN; + continue; + } + + // Calculate coefficients of the kernel basis functions + psImage *polyValues = spatialPolyValues(kernels->spatialOrder, + stamp->xNorm, stamp->yNorm); // Polynomial terms + for (int j = 0; j < numKernels; j++) { + 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]; + } + } + coeff->data.F64[j] = polynomial; + } + psFree(polyValues); + + // Calculate residuals + psKernel *weight = stamp->weight; // Weight postage stamp + psKernel *target; // Target postage stamp + psArray *convolutions; // Convolution postage stamps for each kernel basis function + switch (mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + convolutions = stamp->convolutions1; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + convolutions = stamp->convolutions2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + +#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++) { + float conv = background; // The value of the convolution + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution + conv += convolution->kernel[y][x] * coeff->data.F64[j]; + } + float diff = target->kernel[y][x] - conv; + deviation += PS_SQR(diff) / weight->kernel[y][x]; +#ifdef TESTING + residual->kernel[y][x] = diff; +#endif + } + } + + deviations->data.F32[i] = sqrtf(devNorm * deviation); + psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); + if (!isfinite(deviations->data.F32[i])) { + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + psTrace("psModules.imcombine", 5, + "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + continue; + } + +#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 + + } + psFree(coeff); + + return deviations; +} + + +int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, const psVector *deviations, + psImage *subMask, float sigmaRej, int footprint) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1); + PS_ASSERT_VECTOR_NON_NULL(deviations, -1); + PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, -1); PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); - PS_ASSERT_VECTOR_NON_NULL(solution, -1); - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1); - PS_ASSERT_PTR_NON_NULL(kernels, -1); - PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + - polyTerms(kernels->bgOrder), -1); - PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, -1); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, -1); - - // Measure deviations - psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps + 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 - { - float background = solution->data.F64[solution->n-1]; // The difference in background - - for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest - if (stamp->status != PM_SUBTRACTION_STAMP_USED) { - continue; - } - - psImage *polyValues = spatialPolyValues(kernels->spatialOrder, stamp->xNorm, - stamp->yNorm); // Polynomial terms - - 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++) { - float conv = background; // The value of the convolution - for (int j = 0; j < numKernels; j++) { - psKernel *convolution = convolutions->data[j]; // Convolution of reference - - // XXX Precalculate these values, so that we're not calculating the same thing for - // every pixel. - 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); - - deviations->data.F32[i] = devNorm * deviation; - if (!isfinite(deviations->data.F32[i])) { - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - psTrace("psModules.imcombine", 5, - "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", - i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); - continue; - } - psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", - i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); - totalSquareDev += PS_SQR(deviations->data.F32[i]); - numStamps++; - -#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 - - } + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + totalSquareDev += PS_SQR(deviations->data.F32[i]); + numStamps++; } if (numStamps == 0) { psError(PS_ERR_UNKNOWN, true, "No good stamps found."); - psFree(deviations); return -1; } if (!isfinite(sigmaRej) || sigmaRej <= 0.0) { + // User just wanted to calculate and record the RMS for posterity psLogMsg("psModules.imcombine", PS_LOG_INFO, "RMS deviation: %f", sqrt(totalSquareDev / numStamps)); - psFree(deviations); return 0; } + + float limit = sigmaRej * sqrt(totalSquareDev / (double)numStamps); // Limit on maximum deviation + psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit); int numRejected = 0; // Number of stamps rejected int numGood = 0; // Number of good stamps double newSquareDev = 0.0; // New square deviation - - float limit = sigmaRej * sqrt(totalSquareDev / numStamps); // Limit on maximum deviation - psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit); - for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -1024,6 +1120,11 @@ stamp->status = PM_SUBTRACTION_STAMP_REJECTED; // Recalculate convolutions - psFree(stamp->convolutions); - stamp->convolutions = NULL; + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + stamp->convolutions1 = stamp->convolutions2 = NULL; + psFree(stamp->image1); + psFree(stamp->image2); + psFree(stamp->weight); + stamp->image1 = stamp->image2 = stamp->weight = NULL; } else { numGood++; @@ -1033,8 +1134,14 @@ } - psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n", - numGood, numRejected, sqrt(totalSquareDev / numStamps), sqrt(newSquareDev / numGood)); - - psFree(deviations); + if (numRejected > 0) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, + "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n", + numGood, numRejected, sqrt(totalSquareDev / numStamps), + sqrt(newSquareDev / (double)numGood)); + } else { + psLogMsg("psModules.imcombine", PS_LOG_INFO, + "%d good stamps; 0 rejected.\nRMS deviation: %f\n", + numGood, sqrt(totalSquareDev / numStamps)); + } return numRejected; @@ -1096,38 +1203,46 @@ } -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, bool useFFT) -{ - PS_ASSERT_IMAGE_NON_NULL(inImage, false); - PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false); +bool pmSubtractionConvolve(pmReadout *out, const pmReadout *ro1, const pmReadout *ro2, + const psImage *subMask, psMaskType blank, const psRegion *region, + const psVector *solution, const pmSubtractionKernels *kernels, + pmSubtractionMode mode, bool useFFT) +{ + PS_ASSERT_PTR_NON_NULL(out, false); + PS_ASSERT_PTR_NON_NULL(ro1, false); + PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); + PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false); + if (ro1->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false); + PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); + } + PS_ASSERT_PTR_NON_NULL(ro2, false); + PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); + PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); + if (ro2->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false); + PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false); + } if (subMask) { PS_ASSERT_IMAGE_NON_NULL(subMask, false); PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false); - } - if (inWeight) { - PS_ASSERT_IMAGE_NON_NULL(inWeight, false); - PS_ASSERT_IMAGE_TYPE(inWeight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(inWeight, inImage, false); - } - PS_ASSERT_PTR_NON_NULL(outImage, false); - if (*outImage) { - PS_ASSERT_IMAGE_NON_NULL(*outImage, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(*outImage, inImage, false); - PS_ASSERT_IMAGE_TYPE(*outImage, PS_TYPE_F32, false); - } - if (outMask && *outMask) { - PS_ASSERT_IMAGE_NON_NULL(*outMask, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false); - PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, ro1->image, false); + } + if (out->image) { + PS_ASSERT_IMAGE_NON_NULL(out->image, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out->image, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out->image, PS_TYPE_F32, false); + } + if (out->mask) { + PS_ASSERT_IMAGE_NON_NULL(out->mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out->mask, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out->mask, PS_TYPE_MASK, false); PS_ASSERT_IMAGE_NON_NULL(subMask, false); } - if (outWeight && *outWeight) { - PS_ASSERT_IMAGE_NON_NULL(*outWeight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(*outWeight, inImage, false); - PS_ASSERT_IMAGE_TYPE(*outWeight, PS_TYPE_F32, false); - PS_ASSERT_IMAGE_NON_NULL(inWeight, false); + if (out->weight) { + PS_ASSERT_IMAGE_NON_NULL(out->weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out->weight, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out->weight, PS_TYPE_F32, false); } PS_ASSERT_VECTOR_NON_NULL(solution, false); @@ -1145,9 +1260,9 @@ return false; } - if (region->x0 < 0 || region->x1 > inImage->numCols || - region->y0 < 0 || region->y1 > inImage->numRows) { + if (region->x0 < 0 || region->x1 > ro1->image->numCols || + region->y0 < 0 || region->y1 > ro1->image->numRows) { psString string = psRegionToString(*region); psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)", - string, inImage->numCols, inImage->numRows); + string, ro1->image->numCols, ro1->image->numRows); psFree(string); return false; @@ -1155,27 +1270,42 @@ } + psImage *image, *weight; // Image and weight map to convolve + switch (mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + image = ro1->image; + weight = ro1->weight; + break; + case PM_SUBTRACTION_MODE_2: + image = ro2->image; + weight = ro2->image; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + int numBackground = polyTerms(kernels->bgOrder); // Number of background terms float background = solution->data.F64[solution->n - numBackground]; // The difference in background - int numCols = inImage->numCols, numRows = inImage->numRows; // Image dimensions - - psImage *convImage = *outImage; // Convolved image + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions + + psImage *convImage = out->image; // Convolved image if (!convImage) { - *outImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); - convImage = *outImage; + out->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + convImage = out->image; } psImage *convMask = NULL; // Convolved mask image - if (outMask && subMask) { - if (!*outMask) { - *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - } - convMask = *outMask; + if (subMask) { + if (!out->mask) { + out->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + } + convMask = out->mask; psImageInit(convMask, 0); } psImage *convWeight = NULL; // Convolved weight (variance) image - if (outWeight && inWeight) { - if (!*outWeight) { - *outWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - convWeight = *outWeight; + if (weight) { + if (!out->weight) { + out->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight = out->weight; psImageInit(convWeight, 0.0); } @@ -1198,4 +1328,19 @@ } + psMaskType maskSource, maskTarget; // Mask values for source and target + switch (mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + maskSource = PM_SUBTRACTION_MASK_BAD_1; + maskTarget = PM_SUBTRACTION_MASK_BAD_2 | PM_SUBTRACTION_MASK_CONVOLVE_1; + break; + case PM_SUBTRACTION_MODE_2: + maskSource = PM_SUBTRACTION_MASK_BAD_2; + maskTarget = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + for (int j = yMin; j < yMax; j += fullSize) { int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest @@ -1211,36 +1356,35 @@ kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues); - if (inWeight) { - kernelWeight = varianceKernel(kernelImage); - } - - + if (weight) { + kernelWeight = varianceKernel(kernelWeight, kernelImage); + } + + psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve if (useFFT) { // Use Fast Fourier Transform to do the convolution // This provides a big speed-up for large kernels - convolveFFT(convImage, inImage, subMask, kernelImage, psRegionSet(i, xSubMax, j, ySubMax), + convolveFFT(convImage, image, subMask, maskSource, kernelImage, subRegion, background, size); - if (inWeight) { - convolveFFT(convWeight, inWeight, subMask, kernelWeight, - psRegionSet(i, xSubMax, j, ySubMax), 0.0, size); + if (weight) { + convolveFFT(convWeight, weight, subMask, maskSource, kernelWeight, subRegion, + 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); + convolveDirect(convImage, image, kernelImage, subRegion, background, size); + if (weight) { + convolveDirect(convWeight, weight, kernelWeight, subRegion, 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))) { - convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; - convImage->data.F32[y][x] = NAN; - if (inWeight) { - convWeight->data.F32[y][x] = NAN; + if (subMask) { + for (int y = j; y < ySubMax; y++) { + for (int x = i; x < xSubMax; x++) { + if (subMask->data.PS_TYPE_MASK_DATA[y][x] & maskTarget) { + convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + convImage->data.F32[y][x] = NAN; + if (weight) { + convWeight->data.F32[y][x] = NAN; + } } }