Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 19164) +++ 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;