Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 15286) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 15443) @@ -4,4 +4,5 @@ #include +#include #include @@ -45,8 +46,9 @@ psFree(stamp->matrix); psFree(stamp->vector); - psFree(stamp->reference); - psFree(stamp->input); + psFree(stamp->image1); + psFree(stamp->image2); psFree(stamp->weight); - psFree(stamp->convolutions); + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); } @@ -65,5 +67,6 @@ // Is this position unmasked? static bool checkStampMask(int x, int y, // Coordinates of stamp - const psImage *mask + const psImage *mask, // Mask + pmSubtractionMode mode // Mode for subtraction ) { @@ -74,7 +77,23 @@ return false; } - return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER | - PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ? - false : true; + + psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value + switch (mode) { + case PM_SUBTRACTION_MODE_1: + case PM_SUBTRACTION_MODE_TARGET: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1; + break; + case PM_SUBTRACTION_MODE_2: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + case PM_SUBTRACTION_MODE_UNSURE: + case PM_SUBTRACTION_MODE_DUAL: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + + return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true; } @@ -140,8 +159,9 @@ stamp->status = PM_SUBTRACTION_STAMP_INIT; - stamp->reference = NULL; - stamp->input = NULL; + stamp->image1 = NULL; + stamp->image2 = NULL; stamp->weight = NULL; - stamp->convolutions = NULL; + stamp->convolutions1 = NULL; + stamp->convolutions2 = NULL; return stamp; @@ -151,5 +171,5 @@ pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image, const psImage *subMask, const psRegion *region, - float threshold, float spacing) + float threshold, float spacing, pmSubtractionMode mode) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -226,5 +246,5 @@ 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) { + if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) { fluxStamp = image->data.F32[y][x]; xStamp = x; @@ -242,10 +262,13 @@ // Reset the postage stamps since we're making a new stamp - psFree(stamp->reference); - psFree(stamp->input); + psFree(stamp->image1); + psFree(stamp->image2); psFree(stamp->weight); - stamp->reference = stamp->input = stamp->weight = NULL; - - stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + stamp->image1 = stamp->image2 = stamp->weight = NULL; + stamp->convolutions1 = stamp->convolutions2 = NULL; + + stamp->status = PM_SUBTRACTION_STAMP_FOUND; numFound++; psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n", @@ -266,5 +289,6 @@ pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux, const psImage *image, const psImage *subMask, - const psRegion *region, float spacing, int exclusionZone) + const psRegion *region, float spacing, pmSubtractionMode mode) + { PS_ASSERT_VECTOR_NON_NULL(x, NULL); @@ -289,46 +313,6 @@ } PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); - PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL); 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 - 1; 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]; - j < numStars && (jx = x->data.F32[jIndex]) < ix + exclusionZone; - 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 @@ -347,8 +331,4 @@ // 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 @@ -357,5 +337,5 @@ continue; } - if (!checkStampMask(xPix, yPix, subMask)) { + if (!checkStampMask(xPix, yPix, subMask, mode)) { // Not a good stamp continue; @@ -386,5 +366,4 @@ } } - psFree(exclude); // Sort the list by flux, with the brightest last @@ -421,35 +400,27 @@ -bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input, +bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, 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); - if (input) { - PS_ASSERT_IMAGE_NON_NULL(input, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); - PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false); - } - if (weight) { - PS_ASSERT_IMAGE_NON_NULL(weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false); - PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); - } + PS_ASSERT_IMAGE_NON_NULL(image1, false); + PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false); + if (image2) { + PS_ASSERT_IMAGE_NON_NULL(image2, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false); + PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); + } + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); PS_ASSERT_INT_NONNEGATIVE(footprint, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - int numCols = reference->numCols, numRows = reference->numRows; // Size of images + int numCols = image1->numCols, numRows = image1->numRows; // Size of images int size = kernelSize + footprint; // Size of postage stamps - - if (!weight) { - // 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->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { + if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) { continue; } @@ -468,19 +439,26 @@ } + // Catch memory leaks --- these should have been freed and NULLed before + assert(stamp->image1 == NULL); + assert(stamp->image2 == NULL); + assert(stamp->weight == NULL); + psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest - psImage *refSub = psImageSubset(reference, region); // Subimage with stamp - stamp->reference = psKernelAllocFromImage(refSub, size, size); - psFree(refSub); // Drop reference - - if (input) { - psImage *inSub = psImageSubset(input, region); // Subimage with stamp - stamp->input = psKernelAllocFromImage(inSub, size, size); - psFree(inSub); // Drop reference + psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp + stamp->image1 = psKernelAllocFromImage(sub1, size, size); + psFree(sub1); // Drop reference + + if (image2) { + psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp + stamp->image2 = psKernelAllocFromImage(sub2, size, size); + psFree(sub2); // Drop reference } psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp stamp->weight = psKernelAllocFromImage(wtSub, size, size); - psFree(wtSub); // Drop reference + psFree(wtSub); // Drop reference + + stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; } @@ -488,5 +466,5 @@ } - +#if 0 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) { @@ -517,7 +495,7 @@ float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame - psFree(stamp->input); - stamp->input = psKernelAlloc(-size, size, -size, size); - psKernel *input = stamp->input; // Target stamp + psFree(stamp->image2); + stamp->image2 = psKernelAlloc(-size, size, -size, size); + psKernel *target = stamp->image2; // Target stamp // Put in a Waussian, just for fun! @@ -525,5 +503,5 @@ for (int u = -size; u <= size; u++) { float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma)); - input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); + target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); } } @@ -533,9 +511,9 @@ return true; } - +#endif pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask, const psRegion *region, float spacing, - int exclusionZone) + pmSubtractionMode mode) { PS_ASSERT_ARRAY_NON_NULL(sources, NULL); @@ -561,5 +539,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, - spacing, exclusionZone); // Stamps to return + spacing, mode); // Stamps to return psFree(x); psFree(y); @@ -572,2 +550,31 @@ return stamps; } + + +pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask, + const psRegion *region, float spacing, + pmSubtractionMode mode) +{ + PS_ASSERT_STRING_NON_EMPTY(filename, NULL); + // Let pmSubtractionStampsSet take care of the rest of the assertions + + const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file + psArray *data = psVectorsReadFromFile(filename, format); + if (!data) { + psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename); + return NULL; + } + psVector *x = data->data[0], *y = data->data[1]; // Stamp positions + psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes + + // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions + psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing, + mode); + psFree(data); + + return stamps; + +}