Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 19307) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 19335) @@ -223,9 +223,17 @@ if (threaded) { psMutexUnlock(image); - psMutexLock(mask); - } - psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask - if (threaded) { - psMutexUnlock(mask); + } + + psImage *subMask; // Subimage mask + if (mask) { + if (threaded) { + psMutexLock(mask); + } + subMask = psImageSubset(mask, border); + if (threaded) { + psMutexUnlock(mask); + } + } else { + subMask = NULL; } @@ -238,15 +246,17 @@ if (threaded) { psMutexUnlock(image); - psMutexLock(mask); - } - psFree(subMask); - if (threaded) { - psMutexUnlock(mask); - } - - // Now, we have to chop off the borders, and stick it in where it belongs - // No locking because we own this one - psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges - + } + + if (mask) { + if (threaded) { + psMutexLock(mask); + } + psFree(subMask); + if (threaded) { + psMutexUnlock(mask); + } + } + + // Now, we have to stick it in where it belongs int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region if (background != 0.0) { @@ -262,5 +272,4 @@ } } - psFree(subConv); psFree(convolved); @@ -303,4 +312,6 @@ float background, // Background value to apply psRegion region, // Region to convolve + psMaskType maskBad, // Value to give bad pixels + psMaskType maskPoor, // Value to give poor pixels float poorFrac, // Fraction for "poor" bool useFFT, // Use FFT to convolve? @@ -313,12 +324,15 @@ } - psMaskType maskSource; // Mask these values when - psMaskType maskTarget; // Set this mask value when convolving + psMaskType subBad; // Bad pixels in subtraction mask + psMaskType subConvBad; // Bad pixels in subtraction mask when convolving + psMaskType subConvPoor; // Poor pixels in subtraction mask 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; + subBad = PM_SUBTRACTION_MASK_BAD_1; + subConvBad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_1; + subConvPoor = PM_SUBTRACTION_MASK_CONVOLVE_1; } else { - maskSource = PM_SUBTRACTION_MASK_BAD_2; - maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2; + subBad = PM_SUBTRACTION_MASK_BAD_2; + subConvBad = PM_SUBTRACTION_MASK_CONVOLVE_BAD_2; + subConvPoor = PM_SUBTRACTION_MASK_CONVOLVE_2; } @@ -327,9 +341,10 @@ // Use Fast Fourier Transform to do the convolution // This provides a big speed-up for large kernels - convolveFFT(convImage, image, subMask, maskSource, *kernelImage, region, background, kernels->size); + convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size); if (weight) { - convolveFFT(convWeight, weight, subMask, maskSource, *kernelWeight, region, 0.0, kernels->size); + convolveFFT(convWeight, weight, subMask, subBad, *kernelWeight, region, 0.0, kernels->size); } } else { + // XXX Direct convolution doesn't account for bad pixels yet convolveDirect(convImage, image, *kernelImage, region, background, kernels->size); if (weight) { @@ -344,16 +359,16 @@ if (box > 0) { int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds - bool threaded = pmSubtractionThreaded(); // Are we running threaded? + if (threaded) { psMutexLock(subMask); } psImage *image = psImageSubset(subMask, psRegionSet(colMin - box, colMax + box, - rowMin - box, rowMax + box)); + rowMin - box, rowMax + box)); // Mask if (threaded) { psMutexUnlock(subMask); } - psImage *convolved = psImageConvolveMask(NULL, image, maskSource, maskTarget, + psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad, -box, box, -box, box); // Convolved subtraction mask @@ -366,18 +381,18 @@ } - int numBytes = (colMax - colMin) * PSELEMTYPE_SIZEOF(PS_TYPE_MASK); // Number of bytes to copy psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns"); psAssert(convolved->numRows - 2 * box == rowMax - rowMin, "Bad number of rows"); - // Since we're copying back into the source, we need to lock - if (threaded) { - psMutexLock(subMask); - } for (int yTarget = rowMin, ySource = box; yTarget < rowMax; yTarget++, ySource++) { - memcpy(&subMask->data.PS_TYPE_MASK_DATA[yTarget][colMin], - &convolved->data.PS_TYPE_MASK_DATA[ySource][box], numBytes); - } - if (threaded) { - psMutexUnlock(subMask); + // Dereference images + psMaskType *target = &convMask->data.PS_TYPE_MASK_DATA[yTarget][colMin]; // Target values + psMaskType *source = &convolved->data.PS_TYPE_MASK_DATA[ySource][box]; // Source values + for (int xTarget = colMin; xTarget < colMax; xTarget++, target++, source++) { + if (*source & subConvBad) { + *target |= maskBad; + } else if (*source & subConvPoor) { + *target |= maskPoor; + } + } } @@ -896,11 +911,11 @@ if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight, - ro1->image, ro1->weight, subMask, kernels, polyValues, background, - *region, poorFrac, useFFT, false); + ro1->image, ro1->weight, subMask, kernels, polyValues, background, *region, + maskBad, maskPoor, poorFrac, useFFT, false); } if (kernels->mode == PM_SUBTRACTION_MODE_2 || 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); + ro2->image, ro2->weight, subMask, kernels, polyValues, background, *region, + maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL); } @@ -909,39 +924,4 @@ psFree(polyValues); - // 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++) { - // 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; - } - } - } - } if ((kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) && ro1->mask) { psMaskType **target = convMask->data.PS_TYPE_MASK_DATA; // Target mask @@ -1160,5 +1140,8 @@ psArrayAdd(args, 1, (pmReadout*)ro1); // Casting away const psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const + // Since adding to the array can impact the reference count, we need to lock + psMutexLock(subMask); psArrayAdd(args, 1, (psImage*)subMask); // Casting away const + psMutexUnlock(subMask); PS_ARRAY_ADD_SCALAR(args, maskBad, PS_TYPE_U8); PS_ARRAY_ADD_SCALAR(args, maskPoor, PS_TYPE_U8);