Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 19148) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 19164) @@ -87,5 +87,5 @@ } pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest - if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, false, true)) { + if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, 0, 1.0, region, kernel, false, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); psFree(convRO); 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); Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 19148) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 19164) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-19 00:44:42 $ + * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-08-22 22:45:36 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -26,13 +26,13 @@ /// Mask values for the subtraction mask typedef enum { - 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 + 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 poor or bad + PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be poor or bad + PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad + PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad + PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border + PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected } pmSubtractionMasks; @@ -104,5 +104,7 @@ const pmReadout *ro2, // Input image 2 psImage *subMask, ///< Subtraction mask (or NULL) - psMaskType blank, ///< Mask value for blank regions + psMaskType maskBad, ///< Mask value to give bad pixels + psMaskType maskPoor, ///< Mask value to give poor pixels + float poorFrac, ///< Fraction for "poor" const psRegion *region, ///< Region to convolve (or NULL) const pmSubtractionKernels *kernels, ///< Kernel parameters Index: trunk/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 19148) +++ trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 19164) @@ -138,5 +138,5 @@ if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1, - PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, + PM_SUBTRACTION_MASK_CONVOLVE_1, -size, size, -size, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); @@ -145,21 +145,7 @@ } if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2, - PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, + PM_SUBTRACTION_MASK_CONVOLVE_2, -size, size, -size, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMask(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 (!psImageConvolveMask(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; Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 19148) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 19164) @@ -121,8 +121,8 @@ int inner, int ringsOrder, int binning, float penalty, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, - int iter, float rej, psMaskType maskBad, psMaskType maskBlank, - float badFrac, pmSubtractionMode mode) + int iter, float rej, psMaskType maskVal, psMaskType maskBad, psMaskType maskPoor, + float poorFrac, float badFrac, pmSubtractionMode subMode) { - if (mode != PM_SUBTRACTION_MODE_2) { + if (subMode != PM_SUBTRACTION_MODE_2) { PM_ASSERT_READOUT_NON_NULL(conv1, false); if (conv1->image) { @@ -139,5 +139,5 @@ } } - if (mode != PM_SUBTRACTION_MODE_1) { + if (subMode != PM_SUBTRACTION_MODE_1) { PM_ASSERT_READOUT_NON_NULL(conv2, false); if (conv2->image) { @@ -190,6 +190,9 @@ PS_ASSERT_INT_POSITIVE(iter, false); PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + // Don't care about maskVal // Don't care about maskBad - // Don't care about maskBlank + // Don't care about maskPoor + PS_ASSERT_FLOAT_LARGER_THAN(poorFrac, 0.0, NULL); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(poorFrac, 1.0, NULL); if (isfinite(badFrac)) { PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL); @@ -230,5 +233,5 @@ memCheck("start"); - subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint, + subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint, badFrac, useFFT); if (!subMask) { @@ -267,8 +270,8 @@ if (sources) { stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, footprint, - stampSpacing, mode); + stampSpacing, subMode); } else if (stampsName && strlen(stampsName) > 0) { stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, footprint, - stampSpacing, mode); + stampSpacing, subMode); } @@ -276,14 +279,14 @@ // doesn't matter. if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing, - size, footprint, mode)) { + size, footprint, subMode)) { goto MATCH_ERROR; } - if (mode == PM_SUBTRACTION_MODE_UNSURE) { + if (subMode == PM_SUBTRACTION_MODE_UNSURE) { // Get backgrounds psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator psVector *buffer = NULL;// Buffer for stats - if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskBad, rng)) { + if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 1."); psFree(bgStats); @@ -293,5 +296,5 @@ } float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1 - if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskBad, rng)) { + if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure background of image 2."); psFree(bgStats); @@ -317,5 +320,5 @@ goto MATCH_ERROR; } - mode = newMode; + subMode = newMode; } @@ -323,5 +326,5 @@ if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold, penalty, mode); + stamps, footprint, optThreshold, penalty, subMode); if (!kernels) { psErrorClear(); @@ -332,5 +335,5 @@ // Not an ISIS/GUNK kernel, or the optimum kernel search failed kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder, penalty, mode); + inner, binning, ringsOrder, penalty, subMode); } @@ -344,18 +347,18 @@ } - if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, - PS_META_DUPLICATE_OK, "Subtraction kernels", mode); + PS_META_DUPLICATE_OK, "Subtraction kernels", subMode); psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_DATA_REGION | PS_META_DUPLICATE_OK, "Region over which subtraction was performed", subRegion); } - if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, - PS_META_DUPLICATE_OK, "Subtraction kernels", mode); + PS_META_DUPLICATE_OK, "Subtraction kernels", subMode); psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_DATA_REGION | PS_META_DUPLICATE_OK, @@ -374,5 +377,5 @@ if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing, - size, footprint, mode)) { + size, footprint, subMode)) { goto MATCH_ERROR; } @@ -434,5 +437,5 @@ stamps = NULL; - if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM, PS_META_DUPLICATE_OK, "Number of good stamps", numStamps); @@ -440,5 +443,5 @@ PS_META_DUPLICATE_OK, "RMS deviation of stamps", rmsStamps); } - if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM, PS_META_DUPLICATE_OK, "Number of good stamps", numStamps); @@ -454,7 +457,7 @@ int fullSize = 2 * size + 1 + 1; // Full size of kernel int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize; - psImage *convKernels = psImageAlloc((mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) * + psImage *convKernels = psImageAlloc((subMode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) * imageSize - 1 + - (mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0), + (subMode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0), imageSize - 1, PS_TYPE_F32); psImageInit(convKernels, NAN); @@ -479,5 +482,5 @@ psFree(kernel); - if (mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_DUAL) { kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, (float)j / (float)KERNEL_MOSAIC, @@ -534,6 +537,6 @@ psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, - true, useFFT)) { + if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBad, maskPoor, poorFrac, + region, kernels, true, useFFT)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); goto MATCH_ERROR; @@ -541,5 +544,5 @@ // Set the variance factors - switch (mode) { + switch (subMode) { case PM_SUBTRACTION_MODE_1: { recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false), @@ -565,5 +568,5 @@ // There is data in the readout now - if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { conv1->data_exists = true; if (conv1->parent) { @@ -572,5 +575,5 @@ } } - if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { conv2->data_exists = true; if (conv2->parent) { @@ -590,5 +593,5 @@ weight = NULL; - if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBlank)) { + if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBad)) { psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); goto MATCH_ERROR; @@ -600,5 +603,5 @@ #ifdef TESTING { - if (mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { psFits *fits = psFitsOpen("convolved1.fits", "w"); psFitsWriteImage(fits, NULL, conv1->image, 0, NULL); @@ -606,5 +609,5 @@ } - if (mode == PM_SUBTRACTION_MODE_2 || mode == PM_SUBTRACTION_MODE_DUAL) { + if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { psFits *fits = psFitsOpen("convolved2.fits", "w"); psFitsWriteImage(fits, NULL, conv2->image, 0, NULL); Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 19148) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 19164) @@ -46,6 +46,8 @@ int iter, ///< Rejection iterations float rej, ///< Rejection threshold - psMaskType maskBad, ///< Value to mask - psMaskType maskBlank, ///< Mask for blank region + psMaskType maskVal, ///< Value to mask for input + psMaskType maskBad, ///< Mask for output bad pixels + psMaskType maskPoor, ///< Mask for output poor pixels + float poorFrac, ///< Fraction for "poor" float badFrac, ///< Maximum fraction of bad input pixels to accept pmSubtractionMode mode ///< Mode of subtraction; may be modified Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 19148) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 19164) @@ -73,5 +73,6 @@ static bool checkStampMask(int x, int y, // Coordinates of stamp const psImage *mask, // Mask - pmSubtractionMode mode // Mode for subtraction + pmSubtractionMode mode, // Mode for subtraction + int footprint // Footprint to check for Bad Stuff ) { @@ -79,19 +80,25 @@ return true; } - if (x < 0 || x >= mask->numCols || y < 0 || y >= mask->numRows) { - return false; - } - - psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value + + bool clean = true; // Is the footprint clean? + int numCols = mask->numCols, numRows = mask->numRows; // Size of image + + // Check the footprint bounds + if (x < footprint || x >= numCols - footprint || y < footprint || y >= numRows - footprint) { + clean = false; + } + + // Determine mask value + psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; switch (mode) { case PM_SUBTRACTION_MODE_1: - maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1; + maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1; break; case PM_SUBTRACTION_MODE_2: - maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2; + maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_2; break; case PM_SUBTRACTION_MODE_UNSURE: case PM_SUBTRACTION_MODE_DUAL: - maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2; + maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; break; default: @@ -99,5 +106,35 @@ } - return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true; + // Check the immediate pixel + if (clean && (mask->data.PS_TYPE_MASK_DATA[y][x] & (maskVal | PM_SUBTRACTION_MASK_REJ))) { + clean = false; + } + + int xMin = PS_MAX(x - footprint, 0), xMax = PS_MIN(x + footprint, numCols - 1); // Bounds in x + int yMin = PS_MAX(y - footprint, 0), yMax = PS_MIN(y + footprint, numRows - 1); // Bounds in y + + // Check the footprint + if (clean) { + for (int j = yMin; j <= yMax; j++) { + for (int i = xMin; i <= xMax; i++) { + if (mask->data.PS_TYPE_MASK_DATA[j][i] & maskVal) { + clean = false; + goto CHECK_STAMP_MASK_DONE; + } + } + } + } +CHECK_STAMP_MASK_DONE: + + if (!clean) { + // Mask the footprint, so we don't select something near it again + for (int j = yMin; j <= yMax; j++) { + for (int i = xMin; i <= xMax; i++) { + mask->data.PS_TYPE_MASK_DATA[j][i] |= PM_SUBTRACTION_MASK_REJ; + } + } + } + + return clean; } @@ -259,5 +296,5 @@ for (int y = subRegion->y0; y <= subRegion->y1; y++) { for (int x = subRegion->x0; x <= subRegion->x1; x++) { - if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) { + if (checkStampMask(x, y, subMask, mode, footprint) && image->data.F32[y][x] > fluxStamp) { fluxStamp = image->data.F32[y][x]; xStamp = x; @@ -353,5 +390,5 @@ continue; } - if (!checkStampMask(xPix, yPix, subMask, mode)) { + if (!checkStampMask(xPix, yPix, subMask, mode, footprint)) { // Not a good stamp psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask",