Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 24058) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 24066) @@ -78,29 +78,50 @@ } -// Is this position unmasked? -static bool checkStampMask(int x, int y, // Coordinates of stamp - int numCols, int numRows, // Size of image - const psImage *mask, // Mask - pmSubtractionMode mode, // Mode for subtraction - int footprint, // Size of stamp - int border // Size of border - ) -{ - // Check the footprint bounds - if (x < border || x >= numCols - border || y < border || y >= numRows - border) { - return false; - } - - if (!mask) { - return true; - } - - // Check the central pixel - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; - return false; - } - - return true; + +// Search a region for a suitable stamp +bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return + float *fluxStamp, // Flux of stamp, to return + const psImage *image1, const psImage *image2, // Images to search + float thresh1, float thresh2, // Thresholds for images + const psImage *subMask, // Subtraction mask + int xMin, int xMax, int yMin, int yMax, // Bounds of search + int numCols, int numRows, // Size of images + int border // Border around image + ) +{ + bool found = false; // Found a suitable stamp? + *fluxStamp = -INFINITY; // Flux of best stamp + + // Ensure we're not going to go outside the bounds of the image + xMin = PS_MAX(border, xMin); + xMax = PS_MIN(numCols - border - 1, xMax); + yMin = PS_MAX(border, yMin); + yMax = PS_MIN(numRows - border - 1, yMax); + + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + if ((image1 && image1->data.F32[y][x] < thresh1) || + (image2 && image2->data.F32[y][x] < thresh2)) { + continue; + } + + if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & + (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) { + return false; + } + + // We take the MIN to attempt to avoid transients in both images + float flux = (image1 && image2) ? PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]) : + ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel + if (flux > *fluxStamp) { + *xStamp = x; + *yStamp = y; + *fluxStamp = flux; + found = true; + } + } + } + + return found; } @@ -314,7 +335,7 @@ numSearch++; - float xStamp = 0, yStamp = 0; // Coordinates of stamp - float fluxStamp = NAN; // Flux of stamp - bool goodStamp = false; // Found a good stamp? + int xStamp = 0, yStamp = 0; // Coordinates of stamp + float fluxStamp = -INFINITY; // Flux of stamp + bool goodStamp = false; // Found a good stamp? // A couple different ways of finding stamps: @@ -324,66 +345,35 @@ psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes - // Take stamp off the top of the (sorted) list - if (xList->n > 0) { - int index = xList->n - 1; // Index of new stamp - xStamp = xList->data.F32[index]; - yStamp = yList->data.F32[index]; - fluxStamp = fluxList->data.F32[index]; + // Take stamps off the top of the (sorted) list + for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) { + int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre // Chop off the top of the list - xList->n = index; - yList->n = index; - fluxList->n = index; - - goodStamp = true; - int x = xStamp + 0.5, y = yStamp + 0.5; // Pixel location - if (image1 && image1->data.F32[y][x] < thresh1) { - goodStamp = false; - } else if (image2 && image2->data.F32[y][x] < thresh2) { - goodStamp = false; - } - } else { - psTrace("psModules.imcombine", 9, "No sources in subregion %d", i); + xList->n = j; + yList->n = j; + fluxList->n = j; + + // Fish around a bit to see if we can find a pixel that isn't masked + psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n", + i, xCentre, yCentre); + + // Search bounds + int search = footprint - size; // Search radius + int xMin = PS_MAX(border, xCentre - search); + int xMax = PS_MIN(numCols - border -1, xCentre + search); + int yMin = PS_MAX(border, yCentre - search); + int yMax = PS_MIN(numRows - border - 1, yCentre + search); + + goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, + subMask, xMin, xMax, yMin, yMax, numCols, numRows, border); } } else { // Use a simple method of automatically finding stamps --- take the highest pixel in the // subregion - fluxStamp = 0.0; psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest - const psImage *image = NULL; // Image for flux of stamp - switch (mode) { - case PM_SUBTRACTION_MODE_1: - image = image2 ? image2 : image1; - break; - case PM_SUBTRACTION_MODE_2: - image = image1 ? image1 : image2; - break; - case PM_SUBTRACTION_MODE_UNSURE: - case PM_SUBTRACTION_MODE_DUAL: - // If both are available, we'll take the MIN value below in the hope of rejecting - // transients in both - image = (image1 && image2) ? NULL : (image1 ? image1 : image2); - default: - psAbort("Unrecognised mode: %x", mode); - } - - for (int y = subRegion->y0; y <= subRegion->y1; y++) { - for (int x = subRegion->x0; x <= subRegion->x1; x++) { - if ((image1 && image1->data.F32[y][x] < thresh1) || - (image2 && image2->data.F32[y][x] < thresh2)) { - continue; - } - float value = image ? image->data.F32[y][x] : - PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]); // Value of image - if (value > fluxStamp && - checkStampMask(x, y, numCols, numRows, subMask, mode, footprint, border)) { - fluxStamp = image->data.F32[y][x]; - xStamp = x; - yStamp = y; - goodStamp = true; - } - } - } + goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, + subMask, subRegion->x0, subRegion->x1, subRegion->y0, subRegion->y1, + numCols, numRows, border); } @@ -453,6 +443,4 @@ region, footprint, spacing); // Stamp list int numStamps = stamps->num; // Number of stamps - int numCols = image->numCols, numRows = image->numRows; // Size of image - int border = footprint + size; // Border size psString ds9name = NULL; // Filename for ds9 region file @@ -480,11 +468,4 @@ psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region", xPix, yPix); - pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "magenta"); - continue; - } - if (!checkStampMask(xPix, yPix, numCols, numRows, subMask, mode, footprint, border)) { - // Not a good stamp - psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask", - xPix, yPix); pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red"); continue;