Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14800) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14801) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.58 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-05 21:20:20 $ + * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-10 20:10:05 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -582,7 +582,8 @@ -bool pmSubtractionCalculateEquation(psArray *stamps, const pmSubtractionKernels *kernels, int footprint) -{ - PS_ASSERT_ARRAY_NON_NULL(stamps, false); +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, + int footprint) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); PS_ASSERT_PTR_NON_NULL(kernels, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); @@ -601,6 +602,6 @@ // We iterate over each stamp, allocate the matrix and vectors if // necessary, and then calculate those matrix/vectors. - for (int i = 0; i < stamps->n; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { continue; @@ -736,12 +737,12 @@ -psVector *pmSubtractionSolveEquation(psVector *solution, const psArray *stamps) -{ - PS_ASSERT_ARRAY_NON_NULL(stamps, NULL); +psVector *pmSubtractionSolveEquation(psVector *solution, const pmSubtractionStampList *stamps) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); // Check inputs, while summing the stamp matrices and vectors long numParams = -1; // Number of parameters - for (int i = 0; i < stamps->n; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest PS_ASSERT_PTR_NON_NULL(stamp, NULL); if (stamp->status != PM_SUBTRACTION_STAMP_USED) { @@ -772,6 +773,6 @@ psVectorInit(sumVector, 0.0); psImageInit(sumMatrix, 0.0); - for (int i = 0; i < stamps->n; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (stamp->status == PM_SUBTRACTION_STAMP_USED) { (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); @@ -810,8 +811,8 @@ -int pmSubtractionRejectStamps(psArray *stamps, psImage *subMask, const psVector *solution, +int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, psImage *subMask, const psVector *solution, int footprint, float sigmaRej, const pmSubtractionKernels *kernels) { - PS_ASSERT_ARRAY_NON_NULL(stamps, -1); + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1); PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); @@ -825,5 +826,5 @@ // Measure deviations - psVector *deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); // Mean deviation for stamps + psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps double totalSquareDev = 0.0; // Total square deviation from zero int numStamps = 0; // Number of used stamps @@ -834,6 +835,6 @@ float background = solution->data.F64[solution->n-1]; // The difference in background - for (int i = 0; i < stamps->n; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // The stamp of interest + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest if (stamp->status != PM_SUBTRACTION_STAMP_USED) { continue; @@ -899,6 +900,6 @@ int numRejected = 0; - for (int i = 0; i < stamps->n; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) { // Mask out the stamp in the image so you it's not found again Index: /trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.h (revision 14800) +++ /trunk/psModules/src/imcombine/pmSubtraction.h (revision 14801) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-05 00:15:28 $ + * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-10 20:10:05 $ * * Copyright 2004-207 Institute for Astronomy, University of Hawaii @@ -47,5 +47,5 @@ /// Calculate the least-squares equation to match the image quality -bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps const pmSubtractionKernels *kernels, ///< Kernel parameters int footprint ///< Half-size of region over which to calculate equation @@ -54,9 +54,9 @@ /// Solve the least-squares equation to match the image quality psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL - const psArray *stamps ///< Array of stamps + const pmSubtractionStampList *stamps ///< Stamps ); /// Reject stamps -int pmSubtractionRejectStamps(psArray *stamps, ///< Array of stamps to check for rejection +int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, ///< Stamps psImage *subMask, ///< Subtraction mask const psVector *solution, ///< Solution vector Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 14800) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 14801) @@ -46,5 +46,5 @@ -static bool getStamps(psArray **stamps, // Stamps to read +static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read const psArray *stampsData, // Stamp data from a file const pmReadout *reference, // Reference readout @@ -84,11 +84,12 @@ psFree(*stamps); - *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, subMask, region); - } else { - psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region, - threshold, stampSpacing); - } - if (!stamps) { + // Apply exclusion zone if we're matching to a nominated PSF; otherwise don't care + *stamps = pmSubtractionSetStamps(xStamp, yStamp, fluxStamp, reference->image, subMask, + region, stampSpacing, input ? 0 : footprint); + } + psTrace("psModules.imcombine", 3, "Finding stamps...\n"); + *stamps = pmSubtractionFindStamps(*stamps, reference->image, subMask, region, + threshold, stampSpacing); + if (!*stamps) { psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); return false; @@ -220,5 +221,5 @@ psRegion *region = NULL; // Iso-kernel region psString regionString = NULL; // String for region - psArray *stamps = NULL; // Stamps for matching PSF + pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF psVector *solution = NULL; // Solution to match PSF pmSubtractionKernels *kernels = NULL; // Kernel basis functions Index: /trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 14800) +++ /trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 14801) @@ -171,6 +171,7 @@ pmSubtractionKernels *pmSubtractionKernelsOptimumISIS(pmSubtractionKernelsType type, int size, int inner, - int spatialOrder, psVector *fwhms, int maxOrder, - const psArray *stamps, int footprint, float tolerance) + int spatialOrder, psVector *fwhms, int maxOrder, + const pmSubtractionStampList *stamps, int footprint, + float tolerance) { if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { @@ -184,5 +185,5 @@ PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); PS_ASSERT_INT_NONNEGATIVE(maxOrder, NULL); - PS_ASSERT_ARRAY_NON_NULL(stamps, NULL); + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, NULL); @@ -208,9 +209,12 @@ // Need to save the stamp inputs --- we're changing the values! - int numStamps = stamps->n; // Number of stamps + int numStamps = stamps->num; // Number of stamps psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs + psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps + psVectorInit(badStamps, 0); for (int i = 0; i < numStamps; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || stamp->status == PM_SUBTRACTION_STAMP_NONE) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { + badStamps->data.U8[i] = 0xff; continue; } @@ -230,6 +234,6 @@ int numPixels = 0; // Number of pixels contributing to chi^2 for (int i = 0; i < numStamps; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || stamp->status == PM_SUBTRACTION_STAMP_NONE) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (badStamps->data.U8[i]) { continue; } @@ -238,4 +242,5 @@ psFree(inputs); psFree(kernels); + psFree(badStamps); return NULL; } @@ -255,4 +260,5 @@ psFree(inputs); psFree(kernels); + psFree(badStamps); return NULL; } @@ -287,9 +293,8 @@ for (int j = 0; j < numStamps; j++) { - pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || - stamp->status == PM_SUBTRACTION_STAMP_NONE) { + if (badStamps->data.U8[j]) { continue; } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint); } @@ -301,9 +306,8 @@ double chi2 = 0.0; // Chi^2 for (int j = 0; j < numStamps; j++) { - pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || - stamp->status == PM_SUBTRACTION_STAMP_NONE) { + if (badStamps->data.U8[j]) { continue; } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint); } @@ -329,4 +333,5 @@ psFree(ranking); psFree(kernels); + psFree(badStamps); return NULL; } @@ -336,9 +341,8 @@ // Remove its contribution, and don't include it in the future. for (int j = 0; j < numStamps; j++) { - pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || - stamp->status == PM_SUBTRACTION_STAMP_NONE) { + if (badStamps->data.U8[j]) { continue; } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint); } @@ -365,4 +369,5 @@ psFree(ranking); psFree(kernels); + psFree(badStamps); return NULL; } @@ -376,4 +381,7 @@ psArray *convNew = psArrayAlloc(numStamps); for (int i = 0; i < numStamps; i++) { + if (badStamps->data.U8[i]) { + continue; + } convNew->data[i] = psArrayAlloc(newSize); } @@ -388,6 +396,9 @@ for (int j = 0; j < numStamps; j++) { + if (badStamps->data.U8[j]) { + continue; + } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest psArray *convolutions = convNew->data[j]; // Convolutions for this stamp - pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]); } @@ -405,9 +416,13 @@ for (int i = 0; i < numStamps; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + if (badStamps->data.U8[i]) { + continue; + } + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest psFree(stamp->convolutions); stamp->convolutions = convNew->data[i]; } + psFree(badStamps); psFree(ranking); Index: /trunk/psModules/src/imcombine/pmSubtractionParams.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionParams.h (revision 14800) +++ /trunk/psModules/src/imcombine/pmSubtractionParams.h (revision 14801) @@ -3,5 +3,6 @@ #include -#include "pmSubtractionKernels.h" +#include +#include /// Generate a set of optimum kernels for ISIS (or GUNK) @@ -12,5 +13,5 @@ psVector *fwhms, ///< Gaussian FWHMs to try int maxOrder, ///< Maximum polynomial order - const psArray *stamps, ///< Stamps + const pmSubtractionStampList *stamps, ///< Stamps int footprint, ///< Convolution footprint for stamps float tolerance ///< Maximum difference in chi^2 Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14800) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14801) @@ -9,11 +9,24 @@ #include "pmSubtractionStamps.h" +#define STAMP_LIST_BUFFER 20 // Number of stamps to add to list at a time + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Private (file-static) functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Free function for pmSubtractionStampList +static void subtractionStampListFree(pmSubtractionStampList *list // Stamp list to free + ) +{ + psFree(list->stamps); + psFree(list->regions); + psFree(list->x); + psFree(list->y); + psFree(list->flux); +} + // Free function for pmSubtractionStamp -void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free - ) +static void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free + ) { psFree(stamp->matrix); @@ -28,5 +41,5 @@ static bool checkStampRegion(int x, int y, // Coordinates of stamp const psRegion *region // Region of interest - ) + ) { if (!region) { @@ -40,5 +53,5 @@ static bool checkStampMask(int x, int y, // Coordinates of stamp const psImage *mask - ) + ) { if (!mask) { @@ -54,5 +67,48 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -pmSubtractionStamp *pmSubtractionStampAlloc(pmSubtractionStampStatus status) +pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, + float spacing) +{ + pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return + psMemSetDeallocator(list, (psFreeFunc)subtractionStampListFree); + + // Get region in which to find stamps: [xMin:xMax,yMin:yMax] + int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows; + if (region) { + xMin = PS_MAX(region->x0, xMin); + xMax = PS_MIN(region->x1, xMax); + yMin = PS_MAX(region->y0, yMin); + yMax = PS_MIN(region->y1, yMax); + } + int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest + int xStamps = xSize / spacing + 1, yStamps = ySize / spacing + 1; // Number of stamps in x and y + + list->num = xStamps * yStamps; + list->stamps = psArrayAlloc(list->num); + list->regions = psArrayAlloc(list->num); + + for (int y = 0, index = 0; y < yStamps; y++) { + int yStart = yMin + y * ySize / (yStamps + 1); // Subregion starts here + int yStop = yMin + (y + 1) * ySize / (yStamps + 1) - 1; // Subregion stops here + assert(yStart >= yMin && yStop < yMax); + + for (int x = 0; x < xStamps; x++, index++) { + int xStart = xMin + x * xSize / (xStamps + 1); // Subregion starts here + int xStop = xMin + (x + 1) * xSize / (xStamps + 1) - 1; // Subregion stops here + assert(xStart >= xMin && xStop < xMax); + + list->stamps->data[index] = pmSubtractionStampAlloc(); + list->regions->data[index] = psRegionAlloc(xStart, xStop, yStart, yStop); + } + } + + list->x = NULL; + list->y = NULL; + list->flux = NULL; + + return list; +} + +pmSubtractionStamp *pmSubtractionStampAlloc(void) { pmSubtractionStamp *stamp = psAlloc(sizeof(pmSubtractionStamp)); // Stamp to return @@ -66,5 +122,5 @@ stamp->matrix = NULL; stamp->vector = NULL; - stamp->status = status; + stamp->status = PM_SUBTRACTION_STAMP_INIT; stamp->reference = NULL; @@ -77,6 +133,7 @@ -psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, - const psRegion *region, float threshold, float spacing) +pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, const psImage *image, + const psImage *subMask, const psRegion *region, + float threshold, float spacing) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -107,89 +164,79 @@ int numRows = image->numRows, numCols = image->numCols; // Size of image - // Get region in which to find stamps: [xMin:xMax,yMin:yMax] - int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows; - if (region) { - xMin = PS_MAX(region->x0, xMin); - xMax = PS_MIN(region->x1, xMax); - yMin = PS_MAX(region->y0, yMin); - yMax = PS_MIN(region->y1, yMax); - } - int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest - - // Number of stamps - int xNumStamps = xSize / spacing + 1; - int yNumStamps = ySize / spacing + 1; - + if (!stamps) { + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, spacing); + } + + int numStamps = stamps->num; // Number of stamp regions int numFound = 0; // Number of stamps found - - if (stamps) { - PS_ASSERT_INT_EQUAL(stamps->n, xNumStamps * yNumStamps, NULL); - // Just in case, check for NULL stamps - for (int i = 0; i < xNumStamps * yNumStamps; i++) { - if (!stamps->data[i]) { - stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_REJECTED); + for (int i = 0; i < numStamps; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + + // Only find a new stamp if we need to + if (stamp->status != PM_SUBTRACTION_STAMP_REJECTED && stamp->status != PM_SUBTRACTION_STAMP_INIT) { + continue; + } + + float xStamp = 0, yStamp = 0; // Coordinates of stamp + float fluxStamp = NAN; // Flux of stamp + bool goodStamp = false; // Found a good stamp? + + // A couple different ways of finding stamps: + if (stamps->x && stamps->y) { + // Get the next stamp from the list + psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists + 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]; + + // Chop off the top of the list + xList->n = index; + yList->n = index; + fluxList->n = index; + + goodStamp = true; } - } - } else { - stamps = psArrayAlloc(xNumStamps * yNumStamps); - for (int i = 0; i < xNumStamps * yNumStamps ; i++) { - stamps->data[i] = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_INIT); - } - } - - for (int j = 0, index = 0; j < yNumStamps; j++) { - for (int i = 0; i < xNumStamps; i++, index++) { - pmSubtractionStamp *stamp = stamps->data[index]; // Stamp of interest - - // Only find a new stamp if we need to - if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || - stamp->status == PM_SUBTRACTION_STAMP_INIT) { - - // Find maximum non-masked value in the image section, - // but don't include a footprint around the edge - float fluxBest = threshold; // Flux of best stamp pixel - int xBest = 0, yBest = 0; // Coordinates of best stamp - - // Bounds of region to search for stamp - int yStart = yMin + j * ySize / (yNumStamps + 1); - int yStop = yMin + (j + 1) * ySize / (yNumStamps + 1) - 1; - int xStart = xMin + i * xSize / (xNumStamps + 1); - int xStop = xMin + (i + 1) * xSize / (xNumStamps + 1) - 1; - assert(yStop < numRows && xStop < numCols && yStart >= 0 && xStart >= 0); - - for (int y = yStart; y <= yStop ; y++) { - for (int x = xStart; x <= xStop ; x++) { - if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxBest) { - fluxBest = image->data.F32[y][x]; - xBest = x; - yBest = y; - } + } else { + // Use a simple method of automatically finding stamps --- take the highest pixel in the subregion + fluxStamp = threshold; + psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest + for (int y = subRegion->y0; y <= subRegion->y1 ; y++) { + for (int x = subRegion->x0; x <= subRegion->y1 ; x++) { + if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) { + fluxStamp = image->data.F32[y][x]; + xStamp = x; + yStamp = y; + goodStamp = true; } } - - stamp->x = xBest; - stamp->y = yBest; - stamp->flux = fluxBest; - - // Reset the postage stamps since we're making a new stamp - psFree(stamp->reference); - psFree(stamp->input); - psFree(stamp->weight); - stamp->reference = stamp->input = stamp->weight = NULL; - - if (fluxBest > threshold) { - stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; - numFound++; - psTrace("psModules.imcombine", 5, "Found stamp in region (%d,%d): %d,%d\n", - i, j, (int)stamp->x, (int)stamp->y); - } else { - stamp->status = PM_SUBTRACTION_STAMP_NONE; - } } } - } - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps in %dx%d regions", - numFound, xNumStamps, yNumStamps); + + if (goodStamp) { + stamp->x = xStamp; + stamp->y = yStamp; + stamp->flux = fluxStamp; + + // Reset the postage stamps since we're making a new stamp + psFree(stamp->reference); + psFree(stamp->input); + psFree(stamp->weight); + stamp->reference = stamp->input = stamp->weight = NULL; + + stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; + numFound++; + psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n", + i, (int)stamp->x, (int)stamp->y); + } else { + stamp->status = PM_SUBTRACTION_STAMP_NONE; + } + } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound); return stamps; @@ -197,6 +244,7 @@ -psArray *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux, - const psImage *subMask, const psRegion *region) +pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux, + const psImage *image, const psImage *subMask, + const psRegion *region, float spacing, int exclusionZone) { PS_ASSERT_VECTOR_NON_NULL(x, false); @@ -209,30 +257,141 @@ PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false); PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, false); + } else { + PS_ASSERT_IMAGE_NON_NULL(image, false); } if (subMask) { PS_ASSERT_IMAGE_NON_NULL(subMask, false); PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); - } - - int num = x->n; // Number of stamps - psArray *stamps = psArrayAllocEmpty(num); // Stamps, to return - - for (int i = 0; i < num; i++) { - int xStamp = x->data.F32[i] + 0.5, yStamp = y->data.F32[i] + 0.5; // Pixel coords of stamp - if (!checkStampRegion(xStamp, yStamp, region)) { - continue; - } - if (!checkStampMask(xStamp, yStamp, subMask)) { - continue; - } - - pmSubtractionStamp *stamp = pmSubtractionStampAlloc(PM_SUBTRACTION_STAMP_CALCULATE); // Stamp - stamp->x = x->data.F32[i]; - stamp->y = y->data.F32[i]; - if (flux) { - stamp->flux = flux->data.F32[i]; - } - stamps->data[stamps->n++] = stamp; - } + if (image) { + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, false); + } + } + + int numStars = x->n; // Number of stars + psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star? + psVectorInit(exclude, 0); + + // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so + // when we're trying to get a good subtraction. + if (exclusionZone > 0) { + psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone); + // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then + // to exclude stars within a nominated distance, we need only examine (i.e., calculate the + // 2-dimensional distance to) other stars in the list that are within that distance of the x + // coordinate of the star of interest. By marking both stars when we find a violation of the + // exclusion zone, we need only search upwards from the x coordinate of the star of interest. + + int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars + psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x + for (int i = 0; i < numStars; i++) { + int iIndex = indexes->data.S32[i]; // Index for star i + if (exclude->data.U8[iIndex]) { + continue; + } + float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i + float jx, jy; // Coordinates for star j + for (int j = i + 1, jIndex = indexes->data.S32[j]; + (jx = x->data.F32[jIndex]) < ix + exclusionZone && j < numStars; + j++, jIndex = indexes->data.S32[j]) { + jy = y->data.F32[jIndex]; + if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) { + exclude->data.U8[iIndex] = 0xff; + exclude->data.U8[jIndex] = 0xff; + psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n", + (int)x->data.F32[iIndex], (int)y->data.F32[iIndex], + (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]); + // Would 'break' here, but there might be more stamps within the exclusion zone. + } + } + } + psFree(indexes); + } + + pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, + region, spacing); // Stamp list + int numStamps = stamps->num; // Number of stamps + + // Initialise the lists + stamps->x = psArrayAlloc(numStamps); + stamps->y = psArrayAlloc(numStamps); + stamps->flux = psArrayAlloc(numStamps); + for (int i = 0; i < numStamps; i++) { + stamps->x->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); + stamps->y->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); + stamps->flux->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); + } + + // Put the stars into their appropriate subregions + for (int i = 0; i < numStars; i++) { + if (exclude->data.U8[i]) { + // Star has been excluded + continue; + } + float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp + int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp + if (!checkStampRegion(xPix, yPix, region)) { + // It's not in the big region + continue; + } + if (!checkStampMask(xPix, yPix, subMask)) { + // Not a good stamp + continue; + } + + for (int j = 0; j < numStamps; j++) { + psRegion *subRegion = stamps->regions->data[j]; // Subregion of interest + if (checkStampRegion(xPix, yPix, subRegion)) { + psVector *xList = stamps->x->data[j], *yList = stamps->y->data[j]; // Pixel lists + psVector *fluxList = stamps->flux->data[j]; // Flux list + + int index = xList->n; // Index of new stamp candidate + xList->data.F32[index] = xStamp; + yList->data.F32[index] = yStamp; + + if (flux) { + fluxList->data.F32[index] = flux->data.F32[i]; + } else { + fluxList->data.F32[index] = image->data.F32[yPix][xPix]; + } + + psVectorExtend(xList, STAMP_LIST_BUFFER, 1); + psVectorExtend(yList, STAMP_LIST_BUFFER, 1); + psVectorExtend(fluxList, STAMP_LIST_BUFFER, 1); + + break; + } + } + } + psFree(exclude); + + // Sort the list by flux, with the brightest last + for (int i = 0; i < numStamps; i++) { + psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Pixel lists + psVector *fluxList = stamps->flux->data[i]; // Flux list + + psVector *indexes = psVectorSortIndex(NULL, fluxList); // Indices to sort flux + int num = indexes->n; // Number of candidate stamps in this subregion + + psVector *xSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of x list + psVector *ySorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of y list + psVector *fluxSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of flux list + for (int j = 0; j < num; j++) { + int k = indexes->data.S32[j]; // Sorted index + xSorted->data.F32[j] = xList->data.F32[k]; + ySorted->data.F32[j] = yList->data.F32[k]; + fluxSorted->data.F32[j] = fluxList->data.F32[k]; + } + psFree(indexes); + + psFree(stamps->x->data[i]); + psFree(stamps->y->data[i]); + psFree(stamps->flux->data[i]); + + stamps->x->data[i] = xSorted; + stamps->y->data[i] = ySorted; + stamps->flux->data[i] = fluxSorted; + } + return stamps; @@ -240,8 +399,8 @@ -bool pmSubtractionExtractStamps(psArray *stamps, psImage *reference, psImage *input, psImage *weight, - int footprint, int kernelSize) -{ - PS_ASSERT_ARRAY_NON_NULL(stamps, false); +bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, psImage *reference, psImage *input, + psImage *weight, int footprint, int kernelSize) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); PS_ASSERT_IMAGE_NON_NULL(reference, false); PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false); @@ -263,10 +422,11 @@ if (!weight) { - // Use the input as a rough approximation to the variance map, and HOPE that it's positive. + // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE + // that it's positive. weight = input ? input : reference; } - for (int i = 0; i < stamps->n; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { continue; @@ -307,16 +467,17 @@ -bool pmSubtractionGenerateStamps(psArray *stamps, float fwhm, int footprint, int kernelSize) -{ - PS_ASSERT_ARRAY_NON_NULL(stamps, false); +bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false); PS_ASSERT_INT_NONNEGATIVE(footprint, false); + PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); int size = kernelSize + footprint; // Size of postage stamps - int num = stamps->n; // Number of stamps + int num = stamps->num; // Number of stamps float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma for (int i = 0; i < num; i++) { - pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) { continue; @@ -347,4 +508,4 @@ } - return stamps; -} + return true; +} Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14800) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14801) @@ -15,4 +15,38 @@ } pmSubtractionStampStatus; +/// A list of stamps +typedef struct { + long num; ///< Number of stamps + psArray *stamps; ///< The stamps + psArray *regions; ///< Regions for each stamp + psArray *x, *y; ///< Coordinates for possible stamps (or NULL) + psArray *flux; ///< Fluxes for possible stamps (or NULL) +} pmSubtractionStampList; + +/// Allocate a list of stamps +pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image + int numRows, // Number of rows in image + const psRegion *region, // Region for stamps, or NULL + float spacing // Rough average spacing between stamps + ); + +/// Assertion for stamp list to be valid +#define PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(LIST, RETURNVALUE) { \ + PS_ASSERT_PTR_NON_NULL(LIST, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->stamps, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->regions, RETURNVALUE); \ + PS_ASSERT_INT_POSITIVE((LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->stamps, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->regions, (LIST)->num, RETURNVALUE); \ + if ((LIST)->x || (LIST)->y || (LIST)->flux) { \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->x, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->y, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->flux, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->x, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->y, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->flux, (LIST)->num, RETURNVALUE); \ + } \ +} + /// A stamp for image subtraction typedef struct { @@ -26,26 +60,35 @@ psImage *matrix; ///< Associated matrix psVector *vector; ///< Assoicated vector - pmSubtractionStampStatus status; ///< Status ofstamp + pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp; +/// Allocate a stamp +pmSubtractionStamp *pmSubtractionStampAlloc(void); + /// Find stamps on an image -psArray *pmSubtractionFindStamps(psArray *stamps, ///< Output stamps, or NULL - const psImage *image, ///< Image for which to find stamps - const psImage *mask, ///< Mask, or NULL - const psRegion *region, ///< Region in which to search for stamps, or NULL - float threshold, ///< Threshold for stamps in the image - float spacing ///< Rough spacing for stamps +pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, ///< Output stamps, or NULL + const psImage *image, ///< Image for which to find stamps + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float threshold, ///< Threshold for stamps in the image + float spacing ///< Rough spacing for stamps ); /// Set stamps based on a list of x,y -psArray *pmSubtractionSetStamps(const psVector *x, ///< x coordinates for each stamp - const psVector *y, ///< y coordinates for each stamp - const psVector *flux, ///< Flux for each stamp, or NULL - const psImage *mask, ///< Mask, or NULL - const psRegion *region ///< Region in which to search for stamps, or NULL +/// +/// We may optionally apply an exclusion zone around each star --- this is important when we're convolving to +/// a specified PSF; less so when we're only trying to get a good subtraction. +pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, ///< x coordinates for each stamp + const psVector *y, ///< y coordinates for each stamp + const psVector *flux, ///< Flux for each stamp, or NULL + const psImage *image, ///< Image for flux of stamp + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float spacing, ///< Rough spacing for stamps + int exclusionZone // Exclusion zone around each star ); /// Extract stamps from the images -bool pmSubtractionExtractStamps(psArray *stamps, ///< Stamps +bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, ///< Stamps psImage *reference, ///< Reference image psImage *input, ///< Input image (or NULL) @@ -56,5 +99,5 @@ /// Generate target for stamps based on list -bool pmSubtractionGenerateStamps(psArray *stamps, ///< Stamps +bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, ///< Stamps float fwhm, ///< FWHM for each stamp int footprint, ///< Stamp footprint size