Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 19148) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 19164) @@ -249,28 +249,19 @@ psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges - if (threaded) { - psMutexLock(target); - } - psImage *subTarget = psImageSubset(target, region); // Target for this subregion - if (threaded) { - psMutexUnlock(target); - } + int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region if (background != 0.0) { - psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32)); + for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { + for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) { + target->data.F32[yTarget][xTarget] = convolved->data.F32[ySource][xSource] + background; + } + } } 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); + int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy + for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) { + memcpy(&target->data.F32[yTarget][xMin], &convolved->data.F32[ySource][size], numBytes); } } psFree(subConv); psFree(convolved); - if (threaded) { - psMutexLock(target); - } - psFree(subTarget); - if (threaded) { - psMutexUnlock(target); - } return; @@ -291,6 +282,5 @@ 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]; + target->data.F32[y][x] += kernel->kernel[v][u] * image->data.F32[y - v][x - u]; } } @@ -303,4 +293,5 @@ static inline void convolveRegion(psImage *convImage, // Convolved image (output) psImage *convWeight, // Convolved weight map (output), or NULL + psImage *convMask, // Convolve mask (output), or NULL psKernel **kernelImage, // Convolution kernel for the image psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL @@ -308,9 +299,9 @@ psImage *weight, // Weight map to convolve, or NULL psImage *subMask, // Subtraction mask - psMaskType maskVal, // Mask value to apply in convolution const pmSubtractionKernels *kernels, // Kernels psImage *polyValues, // Polynomial values float background, // Background value to apply psRegion region, // Region to convolve + float poorFrac, // Fraction for "poor" bool useFFT, // Use FFT to convolve? bool wantDual // Want the dual convolution? @@ -318,15 +309,25 @@ { *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual); - if (weight) { + if (weight || subMask) { *kernelWeight = varianceKernel(*kernelWeight, *kernelImage); } + psMaskType maskSource; // Mask these values when + psMaskType maskTarget; // Set this mask value when convolving + if (kernels->mode == PM_SUBTRACTION_MODE_1 || (kernels->mode == PM_SUBTRACTION_MODE_DUAL && !wantDual)) { + maskSource = PM_SUBTRACTION_MASK_BAD_1; + maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1; + } else { + maskSource = PM_SUBTRACTION_MASK_BAD_2; + maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2; + } + + // Convolve the image and weight if (useFFT) { // Use Fast Fourier Transform to do the convolution // This provides a big speed-up for large kernels - - convolveFFT(convImage, image, subMask, maskVal, *kernelImage, region, background, kernels->size); + convolveFFT(convImage, image, subMask, maskSource, *kernelImage, region, background, kernels->size); if (weight) { - convolveFFT(convWeight, weight, subMask, maskVal, *kernelWeight, region, 0.0, kernels->size); + convolveFFT(convWeight, weight, subMask, maskSource, *kernelWeight, region, 0.0, kernels->size); } } else { @@ -334,4 +335,76 @@ if (weight) { convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size); + } + } + + // 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; + } + } + + if (box > 0) { + bool threaded = pmSubtractionThreaded(); // Are we running threaded? + if (threaded) { + psMutexLock(subMask); + } + psImage *image = psImageSubset(subMask, psRegionSet(region.x0 + xMin, region.x1 + xMax, + region.y0 + yMin, region.y1 + yMax)); + if (threaded) { + psMutexUnlock(subMask); + } + + psImage *convolved = psImageConvolveMask(NULL, image, maskSource, maskTarget, + -box, box, -box, box); // Convolved subtraction mask + + if (threaded) { + psMutexLock(subMask); + } + psFree(image); + if (threaded) { + psMutexUnlock(subMask); + } + + int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds + int numBytes = (colMax - colMin) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy + + // Since we're copying back into the source, we need to lock + if (threaded) { + psMutexLock(subMask); + } + for (int yTarget = rowMin, ySource = -yMin; yTarget < rowMax; yTarget++, ySource++) { + memcpy(&subMask->data.PS_TYPE_MASK_DATA[yTarget][colMin], + &convolved->data.PS_TYPE_MASK_DATA[ySource][-xMin], numBytes); + } + if (threaded) { + psMutexUnlock(subMask); + } + + // No need to lock: we own this + psFree(convolved); } } @@ -753,10 +826,18 @@ // XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data -static bool subtractionConvolvePatch(int numCols, int numRows, int x0, int y0, - pmReadout *out1, pmReadout *out2, psImage *convMask, - const pmReadout *ro1, const pmReadout *ro2, psImage *subMask, - psMaskType blank, psMaskType maskSource, psMaskType maskTarget, - const psRegion *region, const pmSubtractionKernels *kernels, - bool doBG, bool useFFT) +static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image + int x0, int y0, // Offsets for image + pmReadout *out1, pmReadout *out2, // Output readouts + psImage *convMask, // Output convolved mask + const pmReadout *ro1, const pmReadout *ro2, // Input readouts + psImage *subMask, // Input subtraction mask + psMaskType maskBad, // Mask value to give bad pixels + psMaskType maskPoor, // Mask value to give poor pixels + float poorFrac, // Fraction for "poor" + const psRegion *region, // Patch to convolve + const pmSubtractionKernels *kernels, // Kernels + bool doBG, // Add in background when convolving? + bool useFFT // Use FFT to do the convolution? + ) { int size = kernels->size; // Half-size of kernel @@ -782,12 +863,12 @@ if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - convolveRegion(out1->image, out1->weight, &kernelImage, &kernelWeight, ro1->image, ro1->weight, - subMask, maskSource, kernels, polyValues, background, *region, useFFT, - false); + convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight, + ro1->image, ro1->weight, subMask, kernels, polyValues, background, + *region, poorFrac, useFFT, false); } if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - convolveRegion(out2->image, out2->weight, &kernelImage, &kernelWeight, ro2->image, ro2->weight, - subMask, maskSource, kernels, polyValues, background, *region, useFFT, - kernels->mode == PM_SUBTRACTION_MODE_DUAL); + convolveRegion(out2->image, out2->weight, convMask, &kernelImage, &kernelWeight, + ro2->image, ro2->weight, subMask, kernels, polyValues, background, + *region, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL); } @@ -796,13 +877,35 @@ psFree(polyValues); - // Propagate the mask + // Propagate the subtraction mask if (subMask) { psMaskType **target = convMask->data.PS_TYPE_MASK_DATA; // Target mask psMaskType **source = subMask->data.PS_TYPE_MASK_DATA; // Source mask + psMaskType poor, bad; // Mask value for "poor" and "bad" pixels in subtraction mask + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + poor = PM_SUBTRACTION_MASK_CONVOLVE_1; + bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; + break; + case PM_SUBTRACTION_MODE_2: + poor = PM_SUBTRACTION_MASK_CONVOLVE_2; + bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; + break; + case PM_SUBTRACTION_MODE_DUAL: + poor = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_1; + bad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 | + PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; + break; + default: + psAbort("Unrecognised subtraction mode: %x", kernels->mode); + } + for (int y = yMin; y < yMax; y++) { for (int x = xMin; x < xMax; x++) { - if (source[y][x] & maskTarget) { - target[y][x] |= blank; + // Pixels marked as "bad" shouldn't also be marked as "poor" + if (source[y][x] & bad) { + target[y][x] |= maskBad; + } else if (source[y][x] & poor){ + target[y][x] |= maskPoor; } } @@ -848,7 +951,7 @@ const pmReadout *ro2 = args->data[8]; // Input readout 2 psImage *subMask = args->data[9]; // Subtraction mask - psMaskType blank = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value - psMaskType maskSource = PS_SCALAR_VALUE(args->data[11], U8); // Mask value corresponding to source image - psMaskType maskTarget = PS_SCALAR_VALUE(args->data[12], U8); // Mask value corresponding to target image + psMaskType maskBad = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value for bad pixels + psMaskType maskPoor = PS_SCALAR_VALUE(args->data[11], U8); // Output mask value for poor pixels + float poorFrac = PS_SCALAR_VALUE(args->data[12], F32); // Fraction for "poor" const psRegion *region = args->data[13]; // Region to convolve const pmSubtractionKernels *kernels = args->data[14]; // Kernels @@ -856,12 +959,12 @@ bool useFFT = PS_SCALAR_VALUE(args->data[16], U8); // Use FFT for convolution? - return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, - ro1, ro2, subMask, blank, maskSource, maskTarget, - region, kernels, doBG, useFFT); + return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask, + maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT); } bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, - psImage *subMask, psMaskType blank, const psRegion *region, - const pmSubtractionKernels *kernels, bool doBG, bool useFFT) + psImage *subMask, psMaskType maskBad, psMaskType maskPoor, float poorFrac, + const psRegion *region, const pmSubtractionKernels *kernels, + bool doBG, bool useFFT) { int numCols = 0, numRows = 0; // Image dimensions @@ -998,23 +1101,4 @@ yMin = PS_MAX(region->y0, yMin); yMax = PS_MIN(region->y1, yMax); - } - - psMaskType maskSource; // Mask these pixels when convolving - psMaskType maskTarget; // Mark these pixels as bad when propagating the subtractionMask - switch (kernels->mode) { - 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; - case PM_SUBTRACTION_MODE_DUAL: - maskSource = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; - maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; - break; - default: - psAbort("Unsupported subtraction mode: %x", kernels->mode); } @@ -1045,7 +1129,7 @@ psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const psArrayAdd(args, 1, (psImage*)subMask); // Casting away const - PS_ARRAY_ADD_SCALAR(args, blank, PS_TYPE_U8); - PS_ARRAY_ADD_SCALAR(args, maskSource, PS_TYPE_U8); - PS_ARRAY_ADD_SCALAR(args, maskTarget, PS_TYPE_U8); + PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_U8); + PS_ARRAY_ADD_SCALAR(args, maskPoor, PS_TYPE_U8); + PS_ARRAY_ADD_SCALAR(args, poorFrac, PS_TYPE_F32); psArrayAdd(args, 1, subRegion); psArrayAdd(args, 1, (pmSubtractionKernels*)kernels); // Casting away const @@ -1060,5 +1144,5 @@ } else { subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask, - blank, maskSource, maskTarget, subRegion, kernels, doBG, useFFT); + maskBad, maskPoor, poorFrac, subRegion, kernels, doBG, useFFT); } psFree(subRegion);