Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 14801) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 14802) @@ -47,5 +47,4 @@ static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read - const psArray *stampsData, // Stamp data from a file const pmReadout *reference, // Reference readout const pmReadout *input, // Input readout, or NULL to generate fake stamps @@ -60,35 +59,6 @@ ) { - if (stampsData && *stamps) { - // We've already previously read all the stamps - return true; - } - - psImage *inImage = NULL; // Input image - if (input) { - inImage = input->image; - } - - if (stampsData) { - psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes - if (input) { - // We have x, y because the target is provided by the input image - xStamp = stampsData->data[0]; - yStamp = stampsData->data[1]; - } else { - // We have x, y and flux in order to generate a target - xStamp = stampsData->data[0]; - yStamp = stampsData->data[1]; - fluxStamp = stampsData->data[2]; - } - - psFree(*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); + *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing); if (!*stamps) { psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); @@ -98,11 +68,14 @@ memCheck(" find stamps"); - if (!input && !pmSubtractionGenerateStamps(*stamps, targetWidth, footprint, size)) { + if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps."); return false; } + memCheck(" generate stamps"); + psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionExtractStamps(*stamps, reference->image, inImage, weight, footprint, size)) { + if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL, + weight, footprint, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -225,27 +198,4 @@ pmSubtractionKernels *kernels = NULL; // Kernel basis functions - // Read stamps from file - psArray *stampsData = NULL; // Stamps data read from file - if (stampsName && strlen(stampsName) > 0) { - psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName); - if (input) { - // We have x, y because the target is provided by the input image - stampsData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions - } else { - // We have x, y and flux in order to generate a target - stampsData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions - } - if (!stampsData) { - psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName); - return false; - } - - // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions - psVector *xStamp = stampsData->data[0]; // x positions - psVector *yStamp = stampsData->data[1]; // y positions - psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - } - int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions @@ -281,7 +231,42 @@ } + // Read stamps from file + if (stampsName && strlen(stampsName) > 0) { + psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName); + psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes + if (input) { + // We have x, y because the target is provided by the input image + psArray *stampsData = psVectorsReadFromFile(stampsName, "%f %f"); // Stamp positions + if (!stampsData) { + psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName); + goto ERROR; + } + xStamp = psMemIncrRefCounter(stampsData->data[0]); + yStamp = psMemIncrRefCounter(stampsData->data[1]); + psFree(stampsData); + } else { + // We have x, y and flux in order to generate a target + psArray *stampsData = psVectorsReadFromFile(stampsName, "%f %f %f"); // Stamp positions + if (!stampsData) { + psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName); + goto ERROR; + } + xStamp = psMemIncrRefCounter(stampsData->data[0]); + yStamp = psMemIncrRefCounter(stampsData->data[1]); + fluxStamp = psMemIncrRefCounter(stampsData->data[2]); + psFree(stampsData); + } + + // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions + psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + + stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask, + region, stampSpacing, input ? 0 : 2 * footprint); + } + // Define kernel basis functions if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, NULL, + if (!getStamps(&stamps, reference, input, subMask, weight, NULL, threshold, stampSpacing, targetWidth, size, footprint)) { goto ERROR; @@ -308,5 +293,5 @@ psTrace("psModules.imcombine", 2, "Iteration %d...\n", k); - if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, region, + if (!getStamps(&stamps, reference, input, subMask, weight, region, threshold, stampSpacing, targetWidth, size, footprint)) { goto ERROR; @@ -457,6 +442,4 @@ psFree(subMask); subMask = NULL; - psFree(stampsData); - stampsData = NULL; if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) { @@ -477,5 +460,4 @@ psFree(stamps); psFree(solution); - psFree(stampsData); psFree(stamps); return false; Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14801) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14802) @@ -5,4 +5,16 @@ #include #include + +// All these includes required to get stamps out of an array of pmSources +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmResiduals.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmGrowthCurve.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" + #include "pmSubtraction.h" @@ -133,5 +145,5 @@ -pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, const psImage *image, +pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image, const psImage *subMask, const psRegion *region, float threshold, float spacing) @@ -170,4 +182,5 @@ int numStamps = stamps->num; // Number of stamp regions int numFound = 0; // Number of stamps found + int numValid = 0; // Number of valid regions for (int i = 0; i < numStamps; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -177,4 +190,5 @@ continue; } + numValid++; float xStamp = 0, yStamp = 0; // Coordinates of stamp @@ -238,5 +252,7 @@ } - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound); + if (numValid > 0) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound); + } return stamps; @@ -244,28 +260,30 @@ -pmSubtractionStampList *pmSubtractionSetStamps(const psVector *x, const psVector *y, const psVector *flux, +pmSubtractionStampList *pmSubtractionStampsSet(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); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false); - PS_ASSERT_VECTOR_NON_NULL(y, false); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL); if (flux) { - PS_ASSERT_VECTOR_NON_NULL(flux, false); - PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, false); + PS_ASSERT_VECTOR_NON_NULL(flux, NULL); + PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, NULL); } else { - PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); } if (subMask) { - PS_ASSERT_IMAGE_NON_NULL(subMask, false); - PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); if (image) { - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, false); - } - } + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); + } + } + PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); + PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL); int numStars = x->n; // Number of stars @@ -399,5 +417,5 @@ -bool pmSubtractionExtractStamps(pmSubtractionStampList *stamps, psImage *reference, psImage *input, +bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input, psImage *weight, int footprint, int kernelSize) { @@ -467,5 +485,5 @@ -bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) +bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -510,2 +528,41 @@ return true; } + + +pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask, + const psRegion *region, float spacing, + int exclusionZone) +{ + PS_ASSERT_ARRAY_NON_NULL(sources, NULL); + // Let pmSubtractionStampsSet take care of the rest of the assertions + + int numSources = sources->n; // Number of stars + + psVector *x = psVectorAlloc(numSources, PS_TYPE_F32); // x coordinates + psVector *y = psVectorAlloc(numSources, PS_TYPE_F32); // y coordinates + psVector *flux = psVectorAlloc(numSources, PS_TYPE_F32); // Fluxes + + for (int i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (source->modelPSF) { + x->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + y->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + x->data.F32[i] = source->peak->xf; + y->data.F32[i] = source->peak->yf; + } + flux->data.F32[i] = powf(10.0, -0.4 * source->psfMag); + } + + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, + spacing, exclusionZone); // Stamps to return + psFree(x); + psFree(y); + psFree(flux); + + if (!stamps) { + psError(PS_ERR_UNKNOWN, false, "Unable to set stamps from sources."); + } + + return stamps; +} Index: /trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14801) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14802) @@ -67,5 +67,5 @@ /// Find stamps on an image -pmSubtractionStampList *pmSubtractionFindStamps(pmSubtractionStampList *stamps, ///< Output stamps, or NULL +pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL const psImage *image, ///< Image for which to find stamps const psImage *mask, ///< Mask, or NULL @@ -79,5 +79,5 @@ /// 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 +pmSubtractionStampList *pmSubtractionStampsSet(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 @@ -86,9 +86,18 @@ const psRegion *region, ///< Region to search, or NULL float spacing, ///< Rough spacing for stamps - int exclusionZone // Exclusion zone around each star + int exclusionZone ///< Exclusion zone around each star + ); + +/// +pmSubtractionStampList *pmSubtractionStampsSetFromSources( + const psArray *sources, ///< Sources for each stamp + const psImage *subMask, ///< 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(pmSubtractionStampList *stamps, ///< Stamps +bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, ///< Stamps psImage *reference, ///< Reference image psImage *input, ///< Input image (or NULL) @@ -99,5 +108,5 @@ /// Generate target for stamps based on list -bool pmSubtractionGenerateStamps(pmSubtractionStampList *stamps, ///< Stamps +bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, ///< Stamps float fwhm, ///< FWHM for each stamp int footprint, ///< Stamp footprint size