Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 23837) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 23839) @@ -850,4 +850,12 @@ psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit); + + psString ds9name = NULL; // Filename for ds9 region file + static int ds9num = 0; // File number for ds9 region file + psStringAppend(&ds9name, "stamps_reject_%d.ds9", ds9num); + FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "rejected stamps"); + psFree(ds9name); + ds9num++; + int numRejected = 0; // Number of stamps rejected int numGood = 0; // Number of good stamps @@ -872,4 +880,5 @@ } } + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); // Set stamp for replacement @@ -897,8 +906,13 @@ numGood++; newMean += deviations->data.F32[i]; + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); } } } newMean /= numGood; + + if (ds9) { + fclose(ds9); + } if (numRejected > 0) { Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 23837) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 23839) @@ -749,4 +749,11 @@ } + psString ds9name = NULL; // Filename for ds9 region file + static int ds9num = 0; // File number for ds9 region file + psStringAppend(&ds9name, "stamps_solution_%d.ds9", ds9num); + FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "solution stamps"); + psFree(ds9name); + ds9num++; + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { // Accumulate the least-squares matricies and vectors @@ -761,5 +768,8 @@ (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; + } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); } } @@ -830,4 +840,5 @@ (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1); (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; } @@ -990,4 +1001,8 @@ // XXXXX Free temporary matrices and vectors + } + + if (ds9) { + fclose(ds9); } Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 23837) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 23839) @@ -149,4 +149,53 @@ return clean; } + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Functions for generating ds9 region files +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +static bool ds9regions = false; // Save ds9 region files? + +void pmSubtractionRegions(bool state) +{ + ds9regions = state; +} + +FILE *pmSubtractionStampsFile(const pmSubtractionStampList *stamps, const char *filename, + const char *description) +{ + if (!ds9regions || !stamps || !filename) { + return NULL; + } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Writing %s to ds9 region file: %s", + description, filename); + + FILE *file = fopen(filename, "w"); + + // Outline the stamps + for (int i = 0; i < stamps->num; i++) { + psRegion *region = stamps->regions->data[i]; // Region of interest + float xCentre = 0.5 * (region->x0 + region->x1), yCentre = 0.5 * (region->y0 + region->y1); + fprintf(file, "image;box(%f,%f,%f,%f,0) # color=blue\nimage;text(%f,%f,{%d}) # color=blue\n", + xCentre, yCentre, region->x1 - region->x0, region->y1 - region->y0, + xCentre, yCentre, i); + } + + return file; +} + +void pmSubtractionStampPrint(FILE *ds9, float x, float y, float size, const char *color) +{ + if (!ds9regions || !ds9) { + return; + } + fprintf(ds9, "image;circle(%f,%f,%f)", x, y, size); + if (color && strlen(color) > 0) { + fprintf(ds9, " # color=%s", color); + } + fprintf(ds9, "\n"); + return; +} + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -366,34 +415,5 @@ } -bool pmSubtractionStampsRegions(pmSubtractionStampList *stamps, const char *filename) -{ - PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); - PS_ASSERT_STRING_NON_EMPTY(filename, false); - - FILE *outFile = fopen(filename, "w"); // File to write - if (!outFile) { - psError(PS_ERR_IO, false, "Unable to open %s for writing", filename); - return false; - } - - for (int i = 0; i < stamps->num; i++) { - // Outline the stamp - psRegion *region = stamps->regions->data[i]; // Region of interest - fprintf(outFile, "image;box(%f,%f,%f,%f,0)\n", - 0.5 * (region->x0 + region->x1), 0.5 * (region->y0 + region->y1), - region->x1 - region->x0, region->y1 - region->y0); - - psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists - - for (int j = 0; j < xList->n; j++) { - fprintf(outFile, "image;circle(%f,%f,%d)\n", - xList->data.F32[j], yList->data.F32[j], stamps->footprint); - } - } - - fclose(outFile); - - return true; -} + pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, @@ -424,4 +444,10 @@ int numStamps = stamps->num; // Number of stamps + psString ds9name = NULL; // Filename for ds9 region file + static int ds9num = 0; // File number for ds9 region file + psStringAppend(&ds9name, "stamps_all_%d.ds9", ds9num); + FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "all stamps"); // ds9 region file + ds9num++; + // Initialise the lists stamps->x = psArrayAlloc(numStamps); @@ -442,4 +468,5 @@ psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region", xPix, yPix); + pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "magenta"); continue; } @@ -448,4 +475,5 @@ psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask", xPix, yPix); + pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red"); continue; } @@ -471,4 +499,5 @@ psTrace("psModules.imcombine", 9, "Putting input stamp (%d,%d) into subregion %d", xPix, yPix, j); + pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "green"); } } @@ -477,5 +506,10 @@ psTrace("psModules.imcombine", 9, "Unable to find subregion for stamp (%d,%d)", xPix, yPix); - } + pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "yellow"); + } + } + + if (ds9) { + fclose(ds9); } @@ -508,8 +542,4 @@ } - if (psTraceGetLevel("psModules.imcombine") >= 9) { - pmSubtractionStampsRegions(stamps, "stamps_all.ds9"); - } - return stamps; } @@ -581,50 +611,4 @@ } -#if 0 -bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int kernelSize) -{ - PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); - PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false); - PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - - int size = kernelSize + stamps->footprint; // Size of postage 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->stamps->data[i]; // Stamp of interest - if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) { - continue; - } - - float x = stamp->x, y = stamp->y; // Coordinates of stamp - float flux = stamp->flux; // Flux of star - if (!isfinite(flux)) { - psWarning("Unable to generate PSF for stamp %d --- bad flux.", i); - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - continue; - } - - float xStamp = x - (int)(x + 0.5); // x coordinate of star in stamp frame - float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame - - psFree(stamp->image2); - stamp->image2 = psKernelAlloc(-size, size, -size, size); - psKernel *target = stamp->image2; // Target stamp - - // Put in a Waussian, just for fun! - for (int v = -size; v <= size; v++) { - for (int u = -size; u <= size; u++) { - float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma)); - target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); - } - } - - } - - return true; -} -#endif - pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image, const psImage *subMask, const psRegion *region, Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 23837) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 23839) @@ -125,11 +125,26 @@ -/// Write a ds9 region file with stamp positions +/// Turn on/off generation of ds9 region files /// /// Intended for debugging -bool pmSubtractionStampsRegions(pmSubtractionStampList *stamps, ///< Stamps - const char *filename ///< Filename to which to write regions +void pmSubtractionRegions(bool state ///< Generate ds9 region files? ); +/// Open a file for ds9 regions +/// +/// Intended for debugging +FILE *pmSubtractionStampsFile(const pmSubtractionStampList *stamps, ///< List of stamps, for outlines + const char *filename, ///< Filename to write + const char *description ///< Description of file + ); + +/// Print a stamp position to ds9 region file +/// +/// Intended for debugging +void pmSubtractionStampPrint(FILE *ds9, ///< ds9 region file + float x, float y, ///< Position of stamp + float size,///< Size of circle + const char *color ///< Colour + ); #endif