Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 19281) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 19282) @@ -4,4 +4,5 @@ #include +#include #include @@ -14,6 +15,6 @@ -psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold, - const psArray *subRegions, const psArray *kernels) +psPixels *pmStackReject(const psPixels *in, int numCols, int numRows, float threshold, float poorFrac, + const psArray *subRegions, const psArray *subKernels) { PS_ASSERT_PIXELS_NON_NULL(in, NULL); @@ -21,6 +22,6 @@ PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL); PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL); - PS_ASSERT_ARRAY_NON_NULL(kernels, NULL); - PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL); + PS_ASSERT_ARRAY_NON_NULL(subKernels, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, subKernels, NULL); // Trivial case @@ -29,45 +30,19 @@ } - // Get the original image size - int numRegions = subRegions->n; // Number of regions - int numCols = 0, numRows = 0; // Size of original image - int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image + // Check consistency of kernels + int numRegions = subRegions->n; // Number of regions int size = 0; // Size of kernel for (int i = 0; i < numRegions; i++) { - psRegion *subRegion = subRegions->data[i]; // Region of interest - if (subRegion->x0 < minCols) { - minCols = subRegion->x0; - } - if (subRegion->y0 < minRows) { - minRows = subRegion->y0; - } - if (subRegion->x1 > numCols) { - numCols = subRegion->x1; - } - if (subRegion->y1 > numRows) { - numRows = subRegion->y1; - } - - pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest + pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest if (size == 0) { - size = kernel->size; - } else if (kernel->size != size) { + size = kernels->size; + } else if (kernels->size != size) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d", - size, kernel->size); + size, kernels->size); return NULL; } } - // Adjust the size for the size of the subimage - if (valid) { - minCols = PS_MAX(valid->x0, minCols); - minRows = PS_MAX(valid->y0, minRows); - numCols = PS_MIN(valid->x1, numCols); - numRows = PS_MIN(valid->y1, numRows); - } - psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows); - - psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1), - 0x01); // Mask + psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols - 1, 0, numRows - 1), 0x01); // Mask psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve psFree(mask); @@ -77,15 +52,8 @@ pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image inRO->image = image; - inRO->col0 = minCols; - inRO->row0 = minRows; for (int i = 0; i < numRegions; i++) { psRegion *region = subRegions->data[i]; // Region of interest - if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 || - region->y0 > valid->y1 || region->y1 < valid->y0)) { - // Region is outside of our sub-image - continue; - } - pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest - if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernel, false, true)) { + pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest + if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernels, false, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); psFree(convRO); @@ -98,9 +66,9 @@ // Image of the kernel at the centre of the region - float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) / - (float)kernel->numCols; - float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) / - (float)kernel->numRows; - psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false); + float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) / + (float)kernels->numCols; + float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) / + (float)kernels->numRows; + psImage *image = pmSubtractionKernelImage(kernels, xNorm, yNorm, false); if (!image) { psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); @@ -118,13 +86,12 @@ // Range for normalisation - int yMin = PS_MAX(minRows, region->y0) - inRO->row0; - int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0; - int xMin = PS_MAX(minCols, region->x0) - inRO->col0; - int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0; + int xMin = PS_MAX(0, region->x0), xMax = PS_MIN(numCols - 1, region->x1); + int yMin = PS_MAX(0, region->y0), yMax = PS_MIN(numRows - 1, region->y1); psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n", sum, xMin, xMax, yMin, yMax); + sum = 1.0 / sum; for (int y = yMin; y <= yMax; y++) { for (int x = xMin; x <= xMax; x++) { - convRO->image->data.F32[y][x] /= sum; + convRO->image->data.F32[y][x] *= sum; } } @@ -152,6 +119,5 @@ for (int x = size; x < convolved->numCols - size; x++) { if (convolved->data.F32[y][x] > threshold) { - // Pixel coordinates in "bad" correspond to the full image - bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x + minCols, y + minRows); + bad = psPixelsAdd(bad, bad->nalloc, x, y); } } @@ -160,35 +126,52 @@ // Now, grow the mask to include everything that touches a bad pixel in the convolution - int x0 = minCols, y0 = minRows; // Offset for mask image - mask = psPixelsToMask(NULL, bad, psRegionSet(x0, numCols - 1, y0, numRows - 1), 0xff); - for (int i = 0; i < bad->n; i++) { - int xPix = bad->data[i].x - x0, yPix = bad->data[i].y - y0; // Coordinates in frame of mask image - // Convolution limits - int xMin = PS_MAX(xPix - size, 0); - int xMax = PS_MIN(xPix + size, mask->numCols - 1); - int yMin = PS_MAX(yPix - size, 0); - int yMax = PS_MIN(yPix + size, mask->numRows - 1); - for (int y = yMin; y <= yMax; y++) { - for (int x = xMin; x <= xMax; x++) { - assert(x < mask->numCols && y < mask->numRows); - mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff; + // Mask values: + // 0x0f: we think this is bad + // 0xf0: this is within the convolution kernel of a bad pixel + mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0x0f); + for (int i = 0; i < subRegions->n; i++) { + psRegion *region = subRegions->data[i]; // Subtraction region + pmSubtractionKernels *kernels = subKernels->data[i]; // Subtraction kernel + + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1; // Full size of kernel + + // Get region for convolution: [xMin:xMax,yMin:yMax] + int xMin = PS_MAX(region->x0, size), xMax = PS_MIN(region->x1, numCols - size); + int yMin = PS_MAX(region->y0, size), yMax = PS_MIN(region->y1, numRows - size); + + psImage *polyValues = NULL; // Pre-calculated polynomial values + for (int j = yMin; j < yMax; j += fullSize) { + int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest + for (int i = xMin; i < xMax; i += fullSize) { + int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest + + polyValues = p_pmSubtractionPolynomialFromCoords(polyValues, kernels, numCols, numRows, + i + size + 1, j + size + 1); + int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues, + false, poorFrac); // Radius of bad box + if (box > 0) { + // Convolve a subimage, then stick it in the original + psImage *subMask = psImageSubset(mask, psRegionSet(i - box, xSubMax + box, + j - box, ySubMax + box)); // Subimage + psImage *convolved = psImageConvolveMask(NULL, subMask, 0x0f, 0xf0, + -box, box, -box, box); // Convolved mask + psFree(subMask); + + int numBytes = (xSubMax - i) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy + psAssert(convolved->numCols - 2 * box == xSubMax - i, "Bad number of columns"); + psAssert(convolved->numRows - 2 * box == ySubMax - j, "Bad number of rows"); + + for (int yTarget = j, ySource = box; yTarget < ySubMax; yTarget++, ySource++) { + memcpy(&mask->data.PS_TYPE_MASK_DATA[yTarget][i], + &convolved->data.PS_TYPE_MASK_DATA[ySource][box], numBytes); + } + psFree(convolved); + } } } + psFree(polyValues); } bad = psPixelsFromMask(bad, mask, 0xff); - psFree(mask); - - // Convert coordinates to frame of original image - for (int i = 0; i < bad->n; i++) { - int x = bad->data[i].x + x0; - int y = bad->data[i].y + y0; - if (x < 0 || x >= numCols || y < 0 || y >= numRows) { - psWarning("Bad pixel coordinate %d: %d,%d --- ignored.", - i, x, y); - continue; - } - bad->data[i].x = x; - bad->data[i].y = y; - } return bad; Index: trunk/psModules/src/imcombine/pmStackReject.h =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.h (revision 19281) +++ trunk/psModules/src/imcombine/pmStackReject.h (revision 19282) @@ -10,6 +10,7 @@ /// We apply a matched filter to the corresponding mask image, and threshold to find the original pixels psPixels *pmStackReject(const psPixels *in, ///< List of pixels in the convolved image - const psRegion *valid, ///< Valid region to consider + int numCols, int numRows, ///< Size of image of interest float threshold, ///< Threshold on convolved image, 0..1 + float poorFrac, ///< Fraction for "poor" const psArray *regions, ///< Array of image regions for image const psArray *kernels ///< Array of kernel parameters for each region Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 19281) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 19282) @@ -300,5 +300,5 @@ psImage *subMask, // Subtraction mask const pmSubtractionKernels *kernels, // Kernels - psImage *polyValues, // Polynomial values + const psImage *polyValues, // Polynomial values float background, // Background value to apply psRegion region, // Region to convolve @@ -340,32 +340,6 @@ // Convolve the mask for bad pixels if (subMask && convMask) { - psKernel *kernel = *kernelWeight; // Kernel of interest - int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds - - // Determine the thresholds - double sumKernel2 = 0.0; // Sum of the kernel-squared - for (int y = yMin; y <= yMax; y++) { - for (int x = xMin; x <= xMax; x++) { - sumKernel2 += kernel->kernel[y][x]; - } - } - float threshold = sumKernel2 * poorFrac; // Threshold between poor and bad - - // Get bounds of threshold region - // Start with the entire kernel, and keep reducing the size of the box until it drops below threshold - int box = kernels->size; // Size of box with bad pixels - for (double sumBox = sumKernel2; box > 0; box--) { - for (int x = -box; x <= box; x++) { - sumBox -= kernel->kernel[-box][x] + kernel->kernel[box][x]; - } - for (int y = -box + 1; y <= box - 1; y++) { - // Note: not doing corners - sumBox -= kernel->kernel[y][-box] + kernel->kernel[y][box]; - } - if (sumBox < threshold) { - break; - } - } - + int box = p_pmSubtractionBadRadius(*kernelImage, kernels, polyValues, + wantDual, poorFrac); // Size of bad box if (box > 0) { int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds @@ -415,4 +389,6 @@ return; } + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -430,4 +406,65 @@ psFree(conv); return convolved; +} + +int p_pmSubtractionBadRadius(psKernel *preKernel, const pmSubtractionKernels *kernels, + const psImage *polyValues, bool wantDual, float poorFrac) +{ + psKernel *kernel; // Kernel to use + if (!preKernel) { + kernel = solvedKernel(NULL, kernels, polyValues, wantDual); + } else { + kernel = psMemIncrRefCounter(preKernel); + } + PS_ASSERT_IMAGE_NON_NULL(polyValues, -1); + + int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds + + // Determine the threshold between bad and poor + double sumKernel2 = 0.0; // Sum of the kernel-squared + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + sumKernel2 += PS_SQR(kernel->kernel[y][x]); + } + } + float threshold = sumKernel2 * poorFrac; // Threshold between poor and bad + + // Get bounds of threshold region + // Start with the entire kernel, and keep reducing the size of the box until it drops below threshold + int box = kernels->size; // Size of box with bad pixels + for (double sumBox = sumKernel2; box > 0; box--) { + for (int x = -box; x <= box; x++) { + sumBox -= PS_SQR(kernel->kernel[-box][x]) + PS_SQR(kernel->kernel[box][x]); + } + for (int y = -box + 1; y <= box - 1; y++) { + // Note: not doing corners + sumBox -= PS_SQR(kernel->kernel[y][-box]) + PS_SQR(kernel->kernel[y][box]); + } + if (sumBox < threshold) { + break; + } + } + + psFree(kernel); + + return box; +} + +psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels, + int numCols, int numRows, int x, int y) +{ + assert(kernels); + assert(numCols > 0 && numRows > 0); + + // Size to use when calculating normalised coordinates (different from actual size when convolving + // subimage) + int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols); + int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows); + + // Normalised coordinates + float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize; + float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize; + + return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm); } @@ -847,13 +884,4 @@ int xMin = region->x0, xMax = region->x1, yMin = region->y0, yMax = region->y1; // Bounds of patch - // Size to use when calculating normalised coordinates (different from actual size when convolving - // subimage) - int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols); - int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows); - - // Normalised coordinates - float yNorm = 2.0 * (float)(yMin + y0 + size + 1 - yNormSize/2.0) / (float)yNormSize; // Normalised coord - float xNorm = 2.0 * (float)(xMin + x0 + size + 1 - xNormSize/2.0) / (float)xNormSize; // Normalised coord - psKernel *kernelImage = NULL; // Kernel for the images psKernel *kernelWeight = NULL; // Kernel for the weight maps @@ -861,6 +889,7 @@ // Only generate polynomial values every kernel footprint, since we have already assumed // (with the stamps) that it does not vary rapidly on this scale. - psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, - xNorm, yNorm); // Pre-calculated polynomial values + psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows, + xMin + x0 + size + 1, + yMin + y0 + size + 1); float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 19281) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 19282) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-22 22:45:36 $ + * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-29 04:04:19 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -126,4 +126,21 @@ ); +/// Given pixel coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j +/// +/// Same as p_pmSubtractionPolynomial except that the normalisation is applied +psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, ///< Output matrix, or NULL + const pmSubtractionKernels *kernels, ///< Kernel parameters + int numCols, int numRows, ///< Size of image of interest + int x, int y ///< Position of interest + ); + +/// Return the radius from the centre of the convolution kernel that distinguishes "bad" and "poor" pixels +int p_pmSubtractionBadRadius(psKernel *preKernel, ///< Pre-calculated convolution kernel + const pmSubtractionKernels *kernels, ///< Kernel parameters + const psImage *polyValues, ///< Polynomial values + bool wantDual, ///< Calculate for the dual kernel? + float poorFrac ///< Fraction for "poor" + ); + /// @} #endif