Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 15433) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 15443) @@ -52,13 +52,15 @@ // Convolve the image with the kernel --- we're basically applying a matched filter and then thresholding - psImage *convolved = NULL; // Convolved image + pmReadout *convRO = pmReadoutAlloc(NULL); // Readout with convolved image + pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image + inRO->image = image; for (int i = 0; i < numRegions; i++) { psRegion *region = regions->data[i]; // Region of interest psVector *solution = solutions->data[i]; // Solution of interest - if (!pmSubtractionConvolve(&convolved, NULL, NULL, image, NULL, NULL, 0, - region, solution, kernels, true)) { + if (!pmSubtractionConvolve(convRO, inRO, NULL, NULL, 0, region, solution, kernels, true, + PM_SUBTRACTION_MODE_1)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); - psFree(convolved); - psFree(image); + psFree(convRO); + psFree(inRO); return NULL; } @@ -73,6 +75,6 @@ if (!kernel) { psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); - psFree(convolved); - psFree(image); + psFree(convRO); + psFree(inRO); return NULL; } @@ -85,8 +87,10 @@ psFree(kernel); - psImage *subConv = psImageSubset(convolved, *region); // Sub-image of convolved image + psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32)); } - psFree(image); + psFree(inRO); + psImage *convolved = psMemIncrRefCounter(convRO->image); + psFree(convRO); // Threshold the convolved image Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 15433) +++ 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; + } } } Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 15443) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-10-17 02:45:40 $ + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-11-03 02:28:24 $ * * Copyright 2004-207 Institute for Astronomy, University of Hawaii @@ -16,18 +16,24 @@ #include -#include "pmSubtractionKernels.h" -#include "pmSubtractionStamps.h" + +#include +#include +#include +#include /// @addtogroup imcombine Image Combinations /// @{ +/// Mask values for the subtraction mask typedef enum { - PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking - PM_SUBTRACTION_MASK_REF = 0x01, // Reference image is bad - PM_SUBTRACTION_MASK_INPUT = 0x02, // Input image is bad - PM_SUBTRACTION_MASK_CONVOLVE = 0x04, // If convolved, would be bad - PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint - PM_SUBTRACTION_MASK_BORDER = 0x10, // Image border - PM_SUBTRACTION_MASK_REJ = 0x20, // Previously tried as a stamp, and rejected + PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking + PM_SUBTRACTION_MASK_BAD_1 = 0x01, // Image 1 is bad + PM_SUBTRACTION_MASK_BAD_2 = 0x02, // Image 2 is bad + PM_SUBTRACTION_MASK_CONVOLVE_1 = 0x04, // If image 1 is convolved, would be bad + PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be bad + PM_SUBTRACTION_MASK_FOOTPRINT_1 = 0x10, // Bad pixel within the stamp footprint of image 1 + PM_SUBTRACTION_MASK_FOOTPRINT_2 = 0x20, // Bad pixel within the stamp footprint of image 2 + PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border + PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected } pmSubtractionMasks; @@ -45,5 +51,6 @@ bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve const pmSubtractionKernels *kernels, ///< Kernel parameters - int footprint ///< Half-size of region over which to calculate equation + int footprint, ///< Half-size of region over which to calculate equation + pmSubtractionMode mode ///< Mode for subtraction (which to convolve) ); @@ -51,5 +58,6 @@ bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps const pmSubtractionKernels *kernels, ///< Kernel parameters - int footprint ///< Half-size of region over which to calculate equation + int footprint, ///< Half-size of region over which to calculate equation + pmSubtractionMode mode ///< Mode for subtraction (which to convolve) ); @@ -59,11 +67,18 @@ ); +/// Calculate deviations +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps + const psVector *solution, ///< Solution vector + int footprint, ///< Half-size of stamp + const pmSubtractionKernels *kernels, ///< Kernel parameters + pmSubtractionMode mode ///< Mode for subtraction + ); + /// Reject stamps int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, ///< Stamps + const psVector *deviations, ///< Deviations for each stamp psImage *subMask, ///< Subtraction mask - const psVector *solution, ///< Solution vector - int footprint, ///< Region to mask if stamp is bad float sigmaRej, ///< Number of RMS deviations above zero at which to reject - const pmSubtractionKernels *kernels ///< Kernel parameters + int footprint ///< Half-size of stamp ); @@ -81,9 +96,7 @@ /// Convolve image in preparation for subtraction -bool pmSubtractionConvolve(psImage **outImage, ///< Output image - psImage **outWeight, ///< Output weight map (or NULL) - psImage **outMask, ///< Output mask (or NULL) - const psImage *inImage, ///< Input image - const psImage *inWeight, ///< Input weight map (or NULL) +bool pmSubtractionConvolve(pmReadout *out, ///< Output image + const pmReadout *ro1, // Input image 1 + const pmReadout *ro2, // Input image 2 const psImage *subMask, ///< Subtraction mask (or NULL) psMaskType blank, ///< Mask value for blank regions @@ -91,4 +104,5 @@ const psVector *solution, ///< The solution vector const pmSubtractionKernels *kernels, ///< Kernel parameters + pmSubtractionMode mode, ///< Mode for subtraction bool useFFT ///< Use Fast Fourier Transform for the convolution? ); Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 15443) @@ -14,4 +14,14 @@ PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction } pmSubtractionKernelsType; + +/// Modes --- specifies which image to convolve +typedef enum { + PM_SUBTRACTION_MODE_ERR, // Error in the mode + PM_SUBTRACTION_MODE_TARGET, // Convolve image 1 to match target PSF + PM_SUBTRACTION_MODE_1, // Convolve image 1 + PM_SUBTRACTION_MODE_2, // Convolve image 2 + PM_SUBTRACTION_MODE_UNSURE, // Not sure yet which image to convolve so try to satisfy both + PM_SUBTRACTION_MODE_DUAL, // Dual convolution +} pmSubtractionMode; /// Kernels specification Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15443) @@ -50,6 +50,6 @@ static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read - const pmReadout *reference, // Reference readout - const pmReadout *input, // Input readout, or NULL to generate fake stamps + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 const psImage *subMask, // Mask for subtraction, or NULL psImage *weight, // Weight map @@ -57,11 +57,11 @@ float threshold, // Threshold for stamp finding float stampSpacing, // Spacing between stamps - float targetWidth,// Target width for fake stamps int size, // Kernel half-size - int footprint // Convolution footprint for stamps + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction ) { psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing); + *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode); if (!*stamps) { psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); @@ -71,14 +71,6 @@ memCheck(" find stamps"); - if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps."); - return false; - } - - memCheck(" generate stamps"); - psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL, - weight, footprint, size)) { + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -94,40 +86,40 @@ -bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input, +bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2, int footprint, float regionSize, float stampSpacing, float threshold, - const psArray *sources, const char *stampsName, float targetWidth, + const psArray *sources, const char *stampsName, pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *isisWidths, const psVector *isisOrders, int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad, - psMaskType maskBlank, float badFrac) + psMaskType maskBlank, float badFrac, pmSubtractionMode mode) { PS_ASSERT_PTR_NON_NULL(convolved, false); - PS_ASSERT_PTR_NON_NULL(reference, false); - PS_ASSERT_IMAGE_NON_NULL(reference->image, false); - PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false); - if (reference->mask) { - PS_ASSERT_IMAGE_NON_NULL(reference->mask, false); - PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false); - } - if (reference->weight) { - PS_ASSERT_IMAGE_NON_NULL(reference->weight, false); - PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false); - } - if (input) { - PS_ASSERT_IMAGE_NON_NULL(input->image, false); - PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false); - if (input->mask) { - PS_ASSERT_IMAGE_NON_NULL(input->mask, false); - PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false); - } - if (input->weight) { - PS_ASSERT_IMAGE_NON_NULL(input->weight, false); - PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, 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->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false); + PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, 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); + } + if (ro2) { + PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); + PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false); + if (ro2->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false); + PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, 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); } } else if (!stampsName && !sources) { @@ -144,5 +136,4 @@ } // stampsName may be anything - // targetWidth can be just about anything (except maybe negative, but it can be NAN) // We'll check kernel type when we allocate the kernels PS_ASSERT_INT_POSITIVE(size, false); @@ -196,20 +187,18 @@ } - psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight - if (input) { - inImage = input->image; - inMask = input->mask; - } - // Where does our weight map come from? psImage *weight = NULL; // Weight image to use - if (input && input->weight) { - weight = input->weight; - } else if (reference->weight) { - weight = reference->weight; - } else if (input) { - weight = input->image; + if (ro1->weight && ro2 && ro2->weight) { + weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight); + } else if (ro1->weight) { + weight = psMemIncrRefCounter(ro1->weight); + } else if (ro2) { + if (ro2->weight) { + weight = psMemIncrRefCounter(ro2->weight); + } else { + weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); + } } else { - weight = reference->image; + weight = psMemIncrRefCounter(ro1->image); } @@ -221,12 +210,13 @@ psVector *solution = NULL; // Solution to match PSF pmSubtractionKernels *kernels = NULL; // Kernel basis functions - - int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions memCheck("start"); - subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint, badFrac, useFFT); + subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint, + badFrac, useFFT); if (!subMask) { psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psFree(weight); return false; } @@ -260,39 +250,41 @@ if (sources) { - stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, - input ? 0 : 2 * footprint); + stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode); } else if (stampsName && strlen(stampsName) > 0) { - // Read stamps from file - psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName); - const char *stampFormat = input ? "%f %f" : "%f %f %f"; // Format for reading stamp file - psArray *stampsData = stampsData = psVectorsReadFromFile(stampsName, stampFormat); - psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes - if (!stampsData) { - psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName); - goto ERROR; - } - xStamp = stampsData->data[0]; - yStamp = stampsData->data[1]; - if (!input) { - fluxStamp = stampsData->data[2]; - } - - // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions - psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - - stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask, - region, stampSpacing, input ? 0 : footprint + size); - psFree(stampsData); + stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode); + } + + // We get the stamps here; we will also attempt to get stamps at the first iteration, but it + // doesn't matter. + if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; + } + + if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) { + pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode + switch (newMode) { + case PM_SUBTRACTION_MODE_1: + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2."); + break; + case PM_SUBTRACTION_MODE_2: + if (mode == PM_SUBTRACTION_MODE_TARGET) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Input PSF is larger than target PSF --- can't match image."); + goto MATCH_ERROR; + } + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); + break; + default: + psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order."); + goto MATCH_ERROR; + } + mode = newMode; } // Define kernel basis functions if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - if (!getStamps(&stamps, reference, input, subMask, weight, NULL, - threshold, stampSpacing, targetWidth, size, footprint)) { - goto ERROR; - } kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold); + stamps, footprint, optThreshold, mode); if (!kernels) { psErrorClear(); @@ -314,13 +306,13 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - if (!getStamps(&stamps, reference, input, subMask, weight, region, - threshold, stampSpacing, targetWidth, size, footprint)) { - goto ERROR; + if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; } psTrace("psModules.imcombine", 3, "Calculating equation...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) { + if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; + goto MATCH_ERROR; } @@ -331,15 +323,27 @@ if (!solution) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; + goto MATCH_ERROR; } memCheck(" solve equation"); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, + mode); // Deviations for each stamp + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + + memCheck(" calculate deviations"); + psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels); + numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); - goto ERROR; - } + psFree(deviations); + goto MATCH_ERROR; + } + psFree(deviations); + memCheck(" reject stamps"); } @@ -350,7 +354,14 @@ if (!solution) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - (void)pmSubtractionRejectStamps(stamps, subMask, solution, footprint, NAN, kernels); + goto MATCH_ERROR; + } + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, + mode); // Deviations for each stamp + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + (void)pmSubtractionRejectStamps(stamps, deviations, subMask, footprint, NAN); + psFree(deviations); } psFree(stamps); @@ -372,5 +383,5 @@ psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); psFree(convKernels); - goto ERROR; + goto MATCH_ERROR; } @@ -380,5 +391,5 @@ psFree(kernel); psFree(convKernels); - goto ERROR; + goto MATCH_ERROR; } psFree(kernel); @@ -416,9 +427,8 @@ psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask, - reference->image, reference->weight, subMask, maskBlank, region, - solution, kernels, useFFT)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image."); - goto ERROR; + if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region, + solution, kernels, mode, useFFT)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); + goto MATCH_ERROR; } psFree(kernels); @@ -461,8 +471,10 @@ psFree(subMask); subMask = NULL; + psFree(weight); + weight = NULL; if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) { psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); - goto ERROR; + goto MATCH_ERROR; } @@ -472,5 +484,5 @@ return true; -ERROR: +MATCH_ERROR: psFree(region); psFree(regionString); @@ -479,4 +491,131 @@ psFree(stamps); psFree(solution); + psFree(weight); return false; } + +// Calculate the second order moments for an image +static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments + int radius // Maximum radius + ) +{ + assert(kernel && kernel->kernel); + + int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x + int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y + + float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment) + float sum = 0.0; // Sum (zero-th moment) + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + xCentroid += kernel->kernel[y][x] * x; + yCentroid += kernel->kernel[y][x] * y; + sum += kernel->kernel[y][x]; + } + } + xCentroid /= sum; + yCentroid /= sum; + + float eta20 = 0.0, eta02 = 0.0; // Second moments + for (int y = yMin; y <= yMax; y++) { + float yDiff = y - yCentroid; + for (int x = xMin; x <= xMax; x++) { + float xDiff = x - xCentroid; + eta20 += PS_SQR(xDiff) * kernel->kernel[y][x]; + eta02 += PS_SQR(yDiff) * kernel->kernel[y][x]; + } + } + + // Normalise to calculate the scale-invariant + float sum2 = PS_SQR(sum); + eta20 /= sum2; + eta02 /= sum2; + // eta11 /= sum2; + + return eta20 + eta02; +} + +#if 0 +// Calculate the deviations for a particular subtraction order +static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel + pmSubtractionStampList *stamps, // Stamps to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int footprint, // Stamp footprint + pmSubtractionMode mode // Mode of subtraction + ) +{ + assert(stamps); + assert(footprint >= 0); + assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2); + + if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + psVector *solution = pmSubtractionSolveEquation(NULL, stamps); + if (!solution) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + if (sumKernel) { + float sum = 0.0; // Sum of the kernel + psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + sum += image->data.F32[y][x]; + } + } + psFree(image); + *sumKernel = sum; + } + + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode); + psFree(solution); + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return NULL; + } + + return deviations; +} +#endif + +pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius) +{ + PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR); + + psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps + psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { + mask->data.PS_TYPE_MASK_DATA[i] = 0xff; + continue; + } + mask->data.PS_TYPE_MASK_DATA[i] = 0; + moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) / + subtractionOrderMoment(stamp->image2, radius); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + if (!psVectorStats(stats, moments, NULL, mask, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio."); + psFree(mask); + psFree(moments); + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(moments); + psFree(mask); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian); + pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2); + psFree(stats); + + return mode; +} + + + Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 15443) @@ -4,12 +4,13 @@ #include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmSubtractionKernels.h" +#include +#include +#include +#include - +/// Match two images bool pmSubtractionMatch(pmReadout *convolved, ///< Output convolved data - const pmReadout *reference, ///< Reference data - const pmReadout *input, ///< Input data + const pmReadout *ro1, ///< Image 1 + const pmReadout *ro2, ///< Image 2 // Stamp parameters int footprint, ///< Stamp half-size @@ -19,5 +20,4 @@ const psArray *sources, ///< Sources for stamps const char *stampsName, ///< Filename for stamps - float targetWidth, ///< Width of PSF for simulated target // Kernel parameters pmSubtractionKernelsType type, ///< Kernel type @@ -38,6 +38,13 @@ psMaskType maskBad, ///< Value to mask psMaskType maskBlank, ///< Mask for blank region - float badFrac ///< Maximum fraction of bad input pixels to accept + float badFrac, ///< Maximum fraction of bad input pixels to accept + pmSubtractionMode mode ///< Mode of subtraction ); +/// Determine which image to convolve +pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, ///< Stamps that have been extracted + int footprint ///< Stamp half-size + ); + + #endif Index: trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 15443) @@ -50,4 +50,21 @@ #endif +/// Select the appropriate convolution, given the kernel basis function and subtraction mode +static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp + int kernelIndex, // Index for kernel component + pmSubtractionMode mode // Mode of subtraction + ) +{ + switch (mode) { + case PM_SUBTRACTION_MODE_1: + return stamp->convolutions1->data[kernelIndex]; + case PM_SUBTRACTION_MODE_2: + return stamp->convolutions2->data[kernelIndex]; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + return NULL; // Unreached +} + // Accumulate cross-term sums for a stamp static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2 @@ -55,14 +72,15 @@ double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 const pmSubtractionStamp *stamp, // Stamp with weight - const psKernel *input, // Input image, I(x) + const psKernel *target, // Target stamp int kernelIndex, // Index for kernel component - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution @@ -82,9 +100,10 @@ const pmSubtractionStamp *stamp, // Stamp with input and weight int kernelIndex, // Index for kernel component - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { @@ -100,18 +119,19 @@ } -static double accumulateChi2(psKernel *input, // Input stamp +static double accumulateChi2(const psKernel *target, // Target stamp pmSubtractionStamp *stamp, // Stamp with weight int kernelIndex, // Index for kernel component double coeff, // Coefficient of convolution double bg, // Background term - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { double chi2 = 0.0; psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution @@ -125,17 +145,28 @@ // Return the initial value of chi^2 -static double initialChi2(psKernel *input, // Input stamp +static double initialChi2(const psKernel *target, // Target stamp const pmSubtractionStamp *stamp, // Stamp with weight - int footprint // Size of convolution + int footprint, // Size of convolution + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight map - psKernel *reference = stamp->reference; // Reference stamp + psKernel *source; // Source stamp + switch (mode) { + case PM_SUBTRACTION_MODE_1: + source = stamp->image1; + break; + case PM_SUBTRACTION_MODE_2: + source = stamp->image2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } double chi2 = 0.0; // Chi^2 for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight - psF32 *ref = &reference->kernel[y][-footprint]; // Derference reference + psF32 *ref = &source->kernel[y][-footprint]; // Derference reference for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { float diff = *in - *ref; // Temporary value @@ -148,16 +179,16 @@ // Subtract a convolution from the input -static void subtractConvolution(psKernel *input, // Input stamp +static void subtractConvolution(psKernel *target, // Target stamp const pmSubtractionStamp *stamp, // Stamp with weight int kernelIndex, // Index for kernel component float coeff, // Coefficient of subtraction float bg, // Background term - int footprint // Size of region of interest - ) -{ - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest - + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction + ) +{ + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, conv++) { @@ -173,5 +204,5 @@ int spatialOrder, const psVector *fwhms, int maxOrder, const pmSubtractionStampList *stamps, int footprint, - float tolerance) + float tolerance, pmSubtractionMode mode) { if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { @@ -210,5 +241,5 @@ // Need to save the stamp inputs --- we're changing the values! int numStamps = stamps->num; // Number of stamps - psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs + psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps psVectorInit(badStamps, 0); @@ -219,7 +250,17 @@ continue; } - 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, size + footprint); + psKernel *target; // Target image of interest + switch (mode) { + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image + targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); psFree(copy); // Drop reference } @@ -238,7 +279,7 @@ continue; } - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); - psFree(inputs); + psFree(targets); psFree(kernels); psFree(badStamps); @@ -258,5 +299,5 @@ "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n", i, (int)stamp->x, (int)stamp->y); - psFree(inputs); + psFree(targets); psFree(kernels); psFree(badStamps); @@ -265,8 +306,8 @@ for (int j = 0; j < numKernels; j++) { - accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint); - } - - lastChi2 += initialChi2(inputs->data[i], stamp, footprint); + accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode); + } + + lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode); numPixels += PS_SQR(2 * footprint + 1); } @@ -297,5 +338,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint); + accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode); } @@ -310,5 +351,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint); + chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode); } @@ -328,5 +369,5 @@ if (bestIndex == -1) { psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter); - psFree(inputs); + psFree(targets); psFree(sumC); psFree(sumCC); @@ -345,5 +386,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint); + subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode); } @@ -361,5 +402,5 @@ lastChi2 = bestChi2; } - psFree(inputs); + psFree(targets); psFree(sumC); psFree(sumCC); @@ -401,5 +442,5 @@ pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest psArray *convolutions = convNew->data[j]; // Convolutions for this stamp - convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]); + convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode)); } } @@ -420,6 +461,18 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - psFree(stamp->convolutions); - stamp->convolutions = convNew->data[i]; + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + switch (mode) { + case PM_SUBTRACTION_MODE_1: + stamp->convolutions1 = convNew->data[i]; + stamp->convolutions2 = NULL; + break; + case PM_SUBTRACTION_MODE_2: + stamp->convolutions1 = NULL; + stamp->convolutions2 = convNew->data[i]; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } } Index: trunk/psModules/src/imcombine/pmSubtractionParams.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.h (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtractionParams.h (revision 15443) @@ -15,5 +15,6 @@ const pmSubtractionStampList *stamps, ///< Stamps int footprint, ///< Convolution footprint for stamps - float tolerance ///< Maximum difference in chi^2 + float tolerance, ///< Maximum difference in chi^2 + pmSubtractionMode mode // Mode for subtraction ); Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 15443) @@ -4,4 +4,5 @@ #include +#include #include @@ -45,8 +46,9 @@ psFree(stamp->matrix); psFree(stamp->vector); - psFree(stamp->reference); - psFree(stamp->input); + psFree(stamp->image1); + psFree(stamp->image2); psFree(stamp->weight); - psFree(stamp->convolutions); + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); } @@ -65,5 +67,6 @@ // Is this position unmasked? static bool checkStampMask(int x, int y, // Coordinates of stamp - const psImage *mask + const psImage *mask, // Mask + pmSubtractionMode mode // Mode for subtraction ) { @@ -74,7 +77,23 @@ return false; } - return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | - PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ? - false : true; + + psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value + switch (mode) { + case PM_SUBTRACTION_MODE_1: + case PM_SUBTRACTION_MODE_TARGET: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1; + break; + case PM_SUBTRACTION_MODE_2: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + case PM_SUBTRACTION_MODE_UNSURE: + case PM_SUBTRACTION_MODE_DUAL: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + + return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true; } @@ -140,8 +159,9 @@ stamp->status = PM_SUBTRACTION_STAMP_INIT; - stamp->reference = NULL; - stamp->input = NULL; + stamp->image1 = NULL; + stamp->image2 = NULL; stamp->weight = NULL; - stamp->convolutions = NULL; + stamp->convolutions1 = NULL; + stamp->convolutions2 = NULL; return stamp; @@ -151,5 +171,5 @@ pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image, const psImage *subMask, const psRegion *region, - float threshold, float spacing) + float threshold, float spacing, pmSubtractionMode mode) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -226,5 +246,5 @@ for (int y = subRegion->y0; y <= subRegion->y1 ; y++) { for (int x = subRegion->x0; x <= subRegion->y1 ; x++) { - if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) { + if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) { fluxStamp = image->data.F32[y][x]; xStamp = x; @@ -242,10 +262,13 @@ // Reset the postage stamps since we're making a new stamp - psFree(stamp->reference); - psFree(stamp->input); + psFree(stamp->image1); + psFree(stamp->image2); psFree(stamp->weight); - stamp->reference = stamp->input = stamp->weight = NULL; - - stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + stamp->image1 = stamp->image2 = stamp->weight = NULL; + stamp->convolutions1 = stamp->convolutions2 = NULL; + + stamp->status = PM_SUBTRACTION_STAMP_FOUND; numFound++; psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n", @@ -266,5 +289,6 @@ pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux, const psImage *image, const psImage *subMask, - const psRegion *region, float spacing, int exclusionZone) + const psRegion *region, float spacing, pmSubtractionMode mode) + { PS_ASSERT_VECTOR_NON_NULL(x, NULL); @@ -289,46 +313,6 @@ } PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); - PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL); int numStars = x->n; // Number of stars - psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star? - psVectorInit(exclude, 0); - - // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so - // when we're trying to get a good subtraction. - if (exclusionZone > 0) { - psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone); - // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then - // to exclude stars within a nominated distance, we need only examine (i.e., calculate the - // 2-dimensional distance to) other stars in the list that are within that distance of the x - // coordinate of the star of interest. By marking both stars when we find a violation of the - // exclusion zone, we need only search upwards from the x coordinate of the star of interest. - - int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars - psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x - for (int i = 0; i < numStars - 1; i++) { - int iIndex = indexes->data.S32[i]; // Index for star i - if (exclude->data.U8[iIndex]) { - continue; - } - float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i - float jx, jy; // Coordinates for star j - for (int j = i + 1, jIndex = indexes->data.S32[j]; - j < numStars && (jx = x->data.F32[jIndex]) < ix + exclusionZone; - j++, jIndex = indexes->data.S32[j]) { - jy = y->data.F32[jIndex]; - if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) { - exclude->data.U8[iIndex] = 0xff; - exclude->data.U8[jIndex] = 0xff; - psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n", - (int)x->data.F32[iIndex], (int)y->data.F32[iIndex], - (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]); - // Would 'break' here, but there might be more stamps within the exclusion zone. - } - } - } - psFree(indexes); - } - pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, region, spacing); // Stamp list @@ -347,8 +331,4 @@ // Put the stars into their appropriate subregions for (int i = 0; i < numStars; i++) { - if (exclude->data.U8[i]) { - // Star has been excluded - continue; - } float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp @@ -357,5 +337,5 @@ continue; } - if (!checkStampMask(xPix, yPix, subMask)) { + if (!checkStampMask(xPix, yPix, subMask, mode)) { // Not a good stamp continue; @@ -386,5 +366,4 @@ } } - psFree(exclude); // Sort the list by flux, with the brightest last @@ -421,35 +400,27 @@ -bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input, +bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, psImage *weight, int footprint, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); - PS_ASSERT_IMAGE_NON_NULL(reference, false); - PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false); - if (input) { - PS_ASSERT_IMAGE_NON_NULL(input, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); - PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false); - } - if (weight) { - PS_ASSERT_IMAGE_NON_NULL(weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false); - PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); - } + PS_ASSERT_IMAGE_NON_NULL(image1, false); + PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false); + if (image2) { + PS_ASSERT_IMAGE_NON_NULL(image2, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false); + PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); + } + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); PS_ASSERT_INT_NONNEGATIVE(footprint, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - int numCols = reference->numCols, numRows = reference->numRows; // Size of images + int numCols = image1->numCols, numRows = image1->numRows; // Size of images int size = kernelSize + footprint; // Size of postage stamps - - if (!weight) { - // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE - // that it's positive. - weight = input ? input : reference; - } for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { + if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) { continue; } @@ -468,19 +439,26 @@ } + // Catch memory leaks --- these should have been freed and NULLed before + assert(stamp->image1 == NULL); + assert(stamp->image2 == NULL); + assert(stamp->weight == NULL); + psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest - psImage *refSub = psImageSubset(reference, region); // Subimage with stamp - stamp->reference = psKernelAllocFromImage(refSub, size, size); - psFree(refSub); // Drop reference - - if (input) { - psImage *inSub = psImageSubset(input, region); // Subimage with stamp - stamp->input = psKernelAllocFromImage(inSub, size, size); - psFree(inSub); // Drop reference + psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp + stamp->image1 = psKernelAllocFromImage(sub1, size, size); + psFree(sub1); // Drop reference + + if (image2) { + psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp + stamp->image2 = psKernelAllocFromImage(sub2, size, size); + psFree(sub2); // Drop reference } psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp stamp->weight = psKernelAllocFromImage(wtSub, size, size); - psFree(wtSub); // Drop reference + psFree(wtSub); // Drop reference + + stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; } @@ -488,5 +466,5 @@ } - +#if 0 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) { @@ -517,7 +495,7 @@ float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame - psFree(stamp->input); - stamp->input = psKernelAlloc(-size, size, -size, size); - psKernel *input = stamp->input; // Target stamp + psFree(stamp->image2); + stamp->image2 = psKernelAlloc(-size, size, -size, size); + psKernel *target = stamp->image2; // Target stamp // Put in a Waussian, just for fun! @@ -525,5 +503,5 @@ for (int u = -size; u <= size; u++) { float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma)); - input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); + target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); } } @@ -533,9 +511,9 @@ return true; } - +#endif pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask, const psRegion *region, float spacing, - int exclusionZone) + pmSubtractionMode mode) { PS_ASSERT_ARRAY_NON_NULL(sources, NULL); @@ -561,5 +539,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, - spacing, exclusionZone); // Stamps to return + spacing, mode); // Stamps to return psFree(x); psFree(y); @@ -572,2 +550,31 @@ return stamps; } + + +pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask, + const psRegion *region, float spacing, + pmSubtractionMode mode) +{ + PS_ASSERT_STRING_NON_EMPTY(filename, NULL); + // Let pmSubtractionStampsSet take care of the rest of the assertions + + const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file + psArray *data = psVectorsReadFromFile(filename, format); + if (!data) { + psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename); + return NULL; + } + psVector *x = data->data[0], *y = data->data[1]; // Stamp positions + psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes + + // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions + psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing, + mode); + psFree(data); + + return stamps; + +} Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 15433) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 15443) @@ -9,7 +9,8 @@ typedef enum { PM_SUBTRACTION_STAMP_INIT, ///< Initial state + PM_SUBTRACTION_STAMP_FOUND, ///< Found a suitable source for this stamp + PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp PM_SUBTRACTION_STAMP_USED, ///< Use this stamp PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected - PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region } pmSubtractionStampStatus; @@ -54,8 +55,9 @@ float flux; ///< Flux float xNorm, yNorm; ///< Normalised position - psKernel *reference; ///< Reference image postage stamp - psKernel *input; ///< Input image postage stamp + psKernel *image1; ///< Reference image postage stamp + psKernel *image2; ///< Input image postage stamp psKernel *weight; ///< Weight image postage stamp - psArray *convolutions; ///< Convolutions of the reference for each kernel component + psArray *convolutions1; ///< Convolutions of image 1 for each kernel component + psArray *convolutions2; ///< Convolutions of image 2 for each kernel component psImage *matrix; ///< Associated matrix psVector *vector; ///< Assoicated vector @@ -72,11 +74,9 @@ const psRegion *region, ///< Region to search, or NULL float threshold, ///< Threshold for stamps in the image - float spacing ///< Rough spacing for stamps + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction ); /// Set stamps based on a list of x,y -/// -/// We may optionally apply an exclusion zone around each star --- this is important when we're convolving to -/// a specified PSF; less so when we're only trying to get a good subtraction. pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp const psVector *y, ///< y coordinates for each stamp @@ -86,14 +86,23 @@ const psRegion *region, ///< Region to search, or NULL float spacing, ///< Rough spacing for stamps - int exclusionZone ///< Exclusion zone around each star + pmSubtractionMode mode ///< Mode for subtraction ); -/// +/// Set stamps based on a list of sources pmSubtractionStampList *pmSubtractionStampsSetFromSources( - const psArray *sources, ///< Sources for each stamp - const psImage *subMask, ///< Mask, or NULL - const psRegion *region, ///< Region to search, or NULL - float spacing, ///< Rough spacing for stamps - int exclusionZone ///< Exclusion zone around each star + const psArray *sources, ///< Sources for each stamp + const psImage *subMask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Set stamps based on values in a file +pmSubtractionStampList *pmSubtractionStampsSetFromFile( + const char *filename, ///< Filename of file containing x,y (or x,y,flux) on each line + const psImage *subMask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction ); @@ -107,10 +116,3 @@ ); -/// Generate target for stamps based on list -bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, ///< Stamps - float fwhm, ///< FWHM for each stamp - int footprint, ///< Stamp footprint size - int size ///< Kernel half-size - ); - #endif