Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 13734) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 13735) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-05-24 21:17:01 $ + * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-06-09 01:04:02 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -333,4 +333,109 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// +psImage *pmSubtractionMask(const psImage *inMask, const psImage *refMask, psMaskType maskVal, + int size, int footprint) +{ + PS_ASSERT_IMAGE_NON_NULL(inMask, NULL); + PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL); + PS_ASSERT_IMAGE_NON_NULL(refMask, NULL); + PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL); + PS_ASSERT_INT_NONNEGATIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); + + // Size of the images + int numCols = inMask->numCols; + int numRows = inMask->numRows; + + // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask + psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask + psImageInit(mask, 0); + + // Dereference for convenience + psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; + psMaskType **inData = inMask->data.PS_TYPE_MASK_DATA; + psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA; + + // Block out a border around the edge of the image + + // Bottom stripe + for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Either side + for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { + for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Top stripe + for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + + // Mask the bad pixels, and around them + for (int y = footprint + size; y < numRows - footprint - size; y++) { + for (int x = 0; x < numCols; x++) { + if (inData[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT; + // Block out the entire stamp footprint around this pixel + for (int v = y - footprint; v <= y + footprint; v++) { + for (int u = x - footprint; u <= x + footprint; u++) { + maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; + } + } + } + if (refData[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_REF; + + // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with + // this bad pixel (within 'size'). Then we want to block out with the FOOTPRINT mask + // everything within a footprint's distance of those (within 'footprint'). + + // Bottom stripe + for (int v = PS_MAX(y - footprint - size, 0); v < PS_MAX(y - size, 0); v++) { + for (int u = PS_MAX(x - footprint - size, 0); + u <= PS_MIN(x + footprint + size, numCols - 1); u++) { + maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; + } + } + // Middle + for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) { + // Left side + for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) { + maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; + } + // Centre + for (int u = PS_MAX(x - size, 0); u <= PS_MIN(x + size, numCols - 1); u++) { + maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_CONVOLVE; + } + // Right side + for (int u = PS_MIN(x + size + 1, numCols - 1); + u <= PS_MIN(x + footprint + size, numCols - 1); u++) { + maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; + } + } + // Top stripe + for (int v = PS_MIN(y + size + 1, numRows - 1); + v <= PS_MIN(y + footprint + size, numRows - 1); v++) { + for (int u = PS_MAX(x - footprint - size, 0); + u <= PS_MIN(x + footprint + size, numCols - 1); u++) { + maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; + } + } + } + } + } + + return mask; +} + bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input, const psImage *weight, const pmSubtractionKernels *kernels, int footprint) @@ -550,6 +655,6 @@ int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage, - psImage *mask, psMaskType badStampMaskVal, const psVector *solution, - int footprint, float sigmaRej, const pmSubtractionKernels *kernels) + psImage *subMask, const psVector *solution, int footprint, float sigmaRej, + const pmSubtractionKernels *kernels) { PS_ASSERT_ARRAY_NON_NULL(stamps, -1); @@ -559,7 +664,7 @@ PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1); PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -1); - PS_ASSERT_IMAGE_NON_EMPTY(mask, -1); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, -1); - PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, -1); + PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); + PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, subMask, -1); PS_ASSERT_VECTOR_NON_NULL(solution, -1); PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1); @@ -654,5 +759,5 @@ for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { - mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal; + subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; } } @@ -707,14 +812,13 @@ bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, - const psImage *inImage, const psImage *inWeight, const psImage *inMask, - psMaskType maskVal, psMaskType blank, - const psVector *solution, const pmSubtractionKernels *kernels) + const psImage *inImage, const psImage *inWeight, const psImage *subMask, + psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels) { PS_ASSERT_IMAGE_NON_NULL(inImage, false); PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false); - if (inMask) { - PS_ASSERT_IMAGE_NON_NULL(inMask, false); - PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, inImage, false); + if (subMask) { + PS_ASSERT_IMAGE_NON_NULL(subMask, false); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false); } if (inWeight) { @@ -733,5 +837,5 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false); PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGE_NON_NULL(inMask, false); + PS_ASSERT_IMAGE_NON_NULL(subMask, false); } if (outWeight && *outWeight) { @@ -763,5 +867,5 @@ } psImage *convMask = NULL; // Convolved mask image - if (outMask && inMask) { + if (outMask && subMask) { if (!*outMask) { *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); @@ -812,18 +916,12 @@ // Check and propagate the kernel footprint, if required - if (inMask) { - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - convMask->data.PS_TYPE_MASK_DATA[y][x] |= - inMask->data.PS_TYPE_MASK_DATA[y + v][x + u]; - } + if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] & + (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) { + convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + convImage->data.F32[y][x] = NAN; + if (inWeight) { + convWeight->data.F32[y][x] = NAN; } - if (convMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) { - convImage->data.F32[y][x] = NAN; - if (inWeight) { - convWeight->data.F32[y][x] = NAN; - } - continue; - } + continue; } Index: /trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.h (revision 13734) +++ /trunk/psModules/src/imcombine/pmSubtraction.h (revision 13735) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-05-11 02:21:01 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-06-09 01:04:02 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -19,4 +19,22 @@ /// @addtogroup imcombine Image Combinations /// @{ + +typedef enum { + PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking + PM_SUBTRACTION_MASK_REF = 0x01, // Reference image is bad + PM_SUBTRACTION_MASK_INPUT = 0x02, // Input image is bad + PM_SUBTRACTION_MASK_CONVOLVE = 0x04, // If convolved, would be bad + PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint + PM_SUBTRACTION_MASK_BORDER = 0x10, // Image border + PM_SUBTRACTION_MASK_REJ = 0x20, // Previously tried as a stamp, and rejected +} pmSubtractionMasks; + +/// Generate a mask for use in the subtraction process +psImage *pmSubtractionMask(const psImage *inMask, ///< Mask for the input image + const psImage *refMask, ///< Mask for the reference image (will be convolved) + psMaskType maskVal, ///< Value to mask out + int size, ///< Half-size of the kernel (pmSubtractionKernels.size) + int footprint ///< Half-size of the kernel footprint + ); /// Calculate the least-squares equation to match the image quality @@ -38,6 +56,5 @@ const psImage *refImage, ///< Reference image psImage *inImage, ///< Input image - psImage *mask, ///< Mask image - psMaskType badStampMaskVal, ///< Value to use in mask for bad stamp + psImage *subMask, ///< Subtraction mask const psVector *solution, ///< Solution vector int footprint, ///< Region to mask if stamp is bad @@ -58,6 +75,5 @@ const psImage *inImage, ///< Input image const psImage *inWeight, ///< Input weight map (or NULL) - const psImage *inMask, ///< Input mask (or NULL) - psMaskType maskVal, ///< Value to mask + const psImage *subMask, ///< Subtraction mask (or NULL) psMaskType blank, ///< Mask value for blank regions const psVector *solution, ///< The solution vector Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 13734) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 13735) @@ -6,4 +6,5 @@ #include +#include "pmSubtraction.h" #include "pmSubtractionStamps.h" @@ -39,28 +40,23 @@ } - -psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *mask, - psMaskType maskVal, psMaskType used, float threshold, - float spacing, int border) +psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, + float threshold, float spacing) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - if (mask) { - PS_ASSERT_IMAGE_NON_NULL(mask, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL); - } - if (used == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "The mask value for used stamps cannot be zero."); - return NULL; + if (subMask) { + PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); } PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); - PS_ASSERT_INT_NONNEGATIVE(border, NULL); - PS_ASSERT_INT_LARGER_THAN(image->numCols, 2 * border, NULL); - PS_ASSERT_INT_LARGER_THAN(image->numRows, 2 * border, NULL); - psMaskType badCentral = maskVal | used; // For central pixel of stamp, don't get a stamp already used - int xNumStamps = (image->numCols - 2 * border) / spacing + 1; // Number of stamps in x - int yNumStamps = (image->numRows - 2 * border) / spacing + 1; // Number of stamps in y + // Size of image + int numRows = image->numRows; + int numCols = image->numCols; + + // Number of stamps + int xNumStamps = numCols / spacing + 1; + int yNumStamps = numRows / spacing + 1; if (stamps) { @@ -78,7 +74,4 @@ } } - // Footprint of image - int numRows = image->numRows; - int numCols = image->numCols; for (int j = 0, index = 0; j < yNumStamps; j++) { @@ -96,39 +89,20 @@ // Bounds of region to search for stamp - int yMin = border + j * (numRows - border) / (yNumStamps + 1); - int yMax = border + (j + 1) * (numRows - border) / (yNumStamps + 1) - 1; - int xMin = border + i * (numCols - border) / (xNumStamps + 1); - int xMax = border + (i + 1) * (numCols - border) / (xNumStamps + 1) - 1; - assert(yMax < image->numRows - border && xMax < image->numCols - border && - yMin >= border && xMin >= border); + int yMin = j * numRows / (yNumStamps + 1); + int yMax = (j + 1) * numRows / (yNumStamps + 1) - 1; + int xMin = i * numCols / (xNumStamps + 1); + int xMax = (i + 1) * numCols / (xNumStamps + 1) - 1; + assert(yMax < image->numRows && xMax < image->numCols && + yMin >= 0 && xMin >= 0); for (int y = yMin; y <= yMax ; y++) { for (int x = xMin; x <= xMax ; x++) { - if (image->data.F32[y][x] > fluxBest) { - bool ok = true; - // Check kernel footprint for bad pixels - if (mask && (mask->data.PS_TYPE_MASK_DATA[y][x] & badCentral)) { - ok = false; - } else { - for (int v = -border; v <= border && ok; v++) { - for (int u = -border; u <= border && ok; u++) { - if ((mask && - (mask->data.PS_TYPE_MASK_DATA[y + v][x + u] & maskVal)) || - !isfinite(image->data.F32[y + v][x + u])) { - ok = false; - if (mask) { - // Mark it so we don't have to look so hard next time - mask->data.PS_TYPE_MASK_DATA[y][x] |= used; - } - } - } - } - } - - if (ok) { - fluxBest = image->data.F32[y][x]; - xBest = x; - yBest = y; - } + if ((!subMask || !(subMask->data.PS_TYPE_MASK_DATA[y][x] & + (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT | + PM_SUBTRACTION_MASK_REJ))) && + image->data.F32[y][x] > fluxBest) { + fluxBest = image->data.F32[y][x]; + xBest = x; + yBest = y; } } Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 13734) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 13735) @@ -3,5 +3,4 @@ #include - /// Status of stamp @@ -26,9 +25,6 @@ const psImage *image, ///< Image for which to find stamps const psImage *mask, ///< Mask - psMaskType maskVal, ///< Value for mask - psMaskType bad, ///< Mask value for bad stamps float threshold, ///< Threshold for stamps in the image - float spacing, ///< Rough spacing for stamps - int border ///< Size of border around image; no stamps in border + float spacing ///< Rough spacing for stamps );