Index: branches/eam_branches/20090522/ippconfig/recipes/psphot.config =================================================================== --- branches/eam_branches/20090522/ippconfig/recipes/psphot.config (revision 24238) +++ branches/eam_branches/20090522/ippconfig/recipes/psphot.config (revision 24529) @@ -11,10 +11,11 @@ SAVE.PLOTS BOOL FALSE +# format for output CMF files +OUTPUT.FORMAT STR PS1_V1 + # the zero point is used to set a basic scale for DVO # XXX it may not currently be read : double check this (EAM) ZERO_POINT F32 25.000 # zero point used by DVO ZERO_PT F32 25.000 # zero point used by DVO - -OUTPUT.FORMAT STR PS1_V1 # these parameter govern how the background is measured @@ -53,4 +54,6 @@ MOMENTS_SY_MAX F32 50.0 MOMENTS_AR_MAX F32 1.5 # maximum axial ratio: 1 / AR < (sx / sy) < AR +MOMENTS_MIN_PIXEL_SN F32 1.0 # XXX caution on this: too high and we will excessively clip +MOMENTS_GAUSS_SIGMA F32 2.0 # XXX TEST THIS: sigma in pixels (should be >> seeing equiv) # basic object statistics Index: branches/eam_branches/20090522/ppSim/src/ppSimArguments.c =================================================================== --- branches/eam_branches/20090522/ppSim/src/ppSimArguments.c (revision 24238) +++ branches/eam_branches/20090522/ppSim/src/ppSimArguments.c (revision 24529) @@ -169,4 +169,5 @@ } ppSimType type = ppSimTypeFromString (typeStr); + // XXX handle error ppSimArgToRecipeStr(&status, options, "FILTER", arguments, "-filter"); // Filter name Index: branches/eam_branches/20090522/ppSim/src/ppSimUtils.c =================================================================== --- branches/eam_branches/20090522/ppSim/src/ppSimUtils.c (revision 24238) +++ branches/eam_branches/20090522/ppSim/src/ppSimUtils.c (revision 24529) @@ -117,4 +117,5 @@ if (!strcasecmp (typeStr, "OBJECT")) return PPSIM_TYPE_OBJECT; psAbort("Should never get here."); + // XXX raise error instead of aborting } Index: branches/eam_branches/20090522/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- branches/eam_branches/20090522/psModules/src/imcombine/pmPSFEnvelope.c (revision 24238) +++ branches/eam_branches/20090522/psModules/src/imcombine/pmPSFEnvelope.c (revision 24529) @@ -298,5 +298,6 @@ } - if (!pmSourceMoments(source, maxRadius)) { + // measure the source moments: tophat windowing, pixel S/N > 1.0 + if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) { // Can't do anything about it; limp along as best we can psErrorClear(); Index: branches/eam_branches/20090522/psModules/src/objects/pmSource.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSource.c (revision 24238) +++ branches/eam_branches/20090522/psModules/src/objects/pmSource.c (revision 24529) @@ -291,7 +291,7 @@ bool status = true; // Status of MD lookup - float PSF_CLUMP_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_SN_LIM"); + float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); if (!status) { - PSF_CLUMP_SN_LIM = 0; + PSF_SN_LIM = 0; } float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); @@ -342,5 +342,5 @@ } - if (src->moments->SN < PSF_CLUMP_SN_LIM) { + if (src->moments->SN < PSF_SN_LIM) { psTrace("psModules.objects", 10, "Rejecting source from clump because of low S/N (%f)\n", src->moments->SN); @@ -450,5 +450,5 @@ if (tmpSrc->moments == NULL) continue; - if (tmpSrc->moments->SN < PSF_CLUMP_SN_LIM) + if (tmpSrc->moments->SN < PSF_SN_LIM) continue; @@ -534,10 +534,9 @@ bool status; float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM"); - if (!status) - PSF_SN_LIM = 20.0; + if (!status) PSF_SN_LIM = 20.0; float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); - if (!status) - PSF_CLUMP_NSIGMA = 1.5; - float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); + if (!status) PSF_CLUMP_NSIGMA = 1.5; + + // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; @@ -582,6 +581,7 @@ source->type = PM_SOURCE_TYPE_STAR; source->mode |= PM_SOURCE_MODE_SATSTAR; - // recalculate moments here with larger box? - pmSourceMoments (source, INNER_RADIUS); + // why do we recalculate moments here? + // we already attempt to do this in psphotSourceStats + // pmSourceMoments (source, INNER_RADIUS); Nsatstar ++; continue; Index: branches/eam_branches/20090522/psModules/src/objects/pmSource.h =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSource.h (revision 24238) +++ branches/eam_branches/20090522/psModules/src/objects/pmSource.h (revision 24529) @@ -213,6 +213,8 @@ */ bool pmSourceMoments( - pmSource *source, ///< The input pmSource for which moments will be computed - float radius ///< Use a circle of pixels around the peak + pmSource *source, ///< The input pmSource for which moments will be computed + float radius, ///< Use a circle of pixels around the peak + float sigma, ///< size of Gaussian window function (<= 0.0 -> skip window) + float minSN ///< minimum pixel significance ); Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceGroup.h =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceGroup.h (revision 24529) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceGroup.h (revision 24529) @@ -0,0 +1,44 @@ +/* @file pmSourceGroup.h + * + * @author EAM, IfA + * + * @version $Revision: $ + * @date $Date: 2009-02-16 22:30:50 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# ifndef PM_SOURCE_GROUP_H +# define PM_SOURCE_GROUP_H + +#include +#include "pmPeaks.h" +#include "pmModel.h" +#include "pmMoments.h" +#include "pmSourceExtendedPars.h" + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +#include + +/** pmSourceGroup data structure + * + * A source group is a connected set of source measurements with a + * common connection. Each source in the group represents the + * detection of some astronomical object on a single each. The group + * represents the related collection of measurements. The fits are + * coupled in some way. For example, they may all have the same + * position, but independent fluxes. Or, they may have a common set + * of positions and shape parameters. Or the position in each image + * may be related by a function. + * + * XXX is thre any info that is neaded for each source group beyond that carried by the sources (besides ID)? + */ +struct pmSource { + int seq; ///< ID for output (generated on write OR set on read) + psArray *sources; +} pmSourceGroup; + + +/// @} +# endif /* PM_SOURCE_GROUP_H */ Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_RAW.c (revision 24238) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_RAW.c (revision 24529) @@ -254,9 +254,12 @@ } + fprintf (f, "# %5s %5s %8s %7s %7s %6s %6s %10s %7s %7s %7s %4s %4s %5s\n", + "x", "y", "peak", "Mx", "My", "Mxx", "Myy", "Sum", "Peak", "Sky", "SN", "nPix", "type", "mode"); + for (i = 0; i < sources->n; i++) { source = sources->data[i]; if (source->moments == NULL) continue; - fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", + fprintf (f, "%5d %5d %8.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", source->peak->x, source->peak->y, source->peak->value, source->moments->Mx, source->moments->My, Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceMoments.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceMoments.c (revision 24238) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceMoments.c (revision 24529) @@ -57,5 +57,5 @@ # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) -bool pmSourceMoments(pmSource *source, psF32 radius) +bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -84,4 +84,6 @@ psF32 Y1 = 0.0; psF32 R2 = PS_SQR(radius); + psF32 minSN2 = PS_SQR(minSN); + psF32 rsigma2 = 0.5 / PS_SQR(sigma); // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -97,4 +99,7 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 yDiff = row - yPeak; + if (fabs(yDiff) > radius) continue; psF32 *vPix = source->pixels->data.F32[row]; @@ -113,5 +118,5 @@ psF32 xDiff = col - xPeak; - psF32 yDiff = row - yPeak; + if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) @@ -121,15 +126,33 @@ psF32 wDiff = *vWgt; - // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) continue; - - Var += wDiff; - Sum += pDiff; - - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; - - X1 += xWght; - Y1 += yWght; + // skip pixels below specified significance level + if (PS_SQR(pDiff) < minSN2*wDiff) continue; + if (pDiff < 0) continue; + + if (sigma > 0.0) { + // apply a pseudo-gaussian weight + + // XXX a lot of extra flops; can we do pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); + psF32 weight = 1.0 / t; + + // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); + + wDiff *= weight; + pDiff *= weight; + } + + Var += wDiff; + Sum += pDiff; + + psF32 xWght = xDiff * pDiff; + psF32 yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + + // fprintf (stderr, " : %6.1f %6.1f %8.1f %8.1f\n", X1, Y1, Sum, Var); peakPixel = PS_MAX (*vPix, peakPixel); @@ -188,4 +211,7 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { + psF32 yDiff = row - yCM; + if (fabs(yDiff) > radius) continue; + psF32 *vPix = source->pixels->data.F32[row]; psF32 *vWgt = source->variance->data.F32[row]; @@ -203,5 +229,5 @@ psF32 xDiff = col - xCM; - psF32 yDiff = row - yCM; + if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) @@ -215,6 +241,21 @@ // XXX EAM : should this limit be user-defined? - if (PS_SQR(pDiff) < wDiff) continue; + if (PS_SQR(pDiff) < minSN2*wDiff) continue; if (pDiff < 0) continue; + + if (sigma > 0.0) { + // apply a pseudo-gaussian weight + + // XXX a lot of extra flops; can we do pre-calculate? + psF32 z = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2; + assert (z >= 0.0); + psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0)); + psF32 weight = 1.0 / t; + + // fprintf (stderr, "%6.1f %6.1f %8.1f %8.1f %5.3f ", xDiff, yDiff, pDiff, wDiff, weight); + + wDiff *= weight; + pDiff *= weight; + } Sum += pDiff; Index: branches/eam_branches/20090522/psModules/src/objects/pmSourceUtils.c =================================================================== --- branches/eam_branches/20090522/psModules/src/objects/pmSourceUtils.c (revision 24238) +++ branches/eam_branches/20090522/psModules/src/objects/pmSourceUtils.c (revision 24529) @@ -94,24 +94,35 @@ source->peak = pmPeakAlloc (xChip, yChip, Io, PM_PEAK_LONE); - int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); - int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); - int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); - int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); + float xReadout, yReadout; - // XXX fix the binning : currently not selected from concepts - // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y - // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y - int xBin = 1; - int yBin = 1; + // if we have information about the chip & cell, adjust the coordinates chip->cell->readout + // otherwise, assume 0,0 offset and 1,1 parity + if (cell) { + int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0"); + int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0"); + int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY"); + int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY"); - // Position on the cell - float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin); - float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin); + // XXX fix the binning : currently not selected from concepts + // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y + // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y + int xBin = 1; + int yBin = 1; - // Position on the readout - // float xReadout = CELL_TO_READOUT(xCell, x0Readout); - // float yReadout = CELL_TO_READOUT(yCell, y0Readout); + // Position on the cell + float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin); + float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin); + + // Position on the readout + // float xReadout = CELL_TO_READOUT(xCell, x0Readout); + // float yReadout = CELL_TO_READOUT(yCell, y0Readout); + xReadout = xCell; + yReadout = yCell; + } else { + xReadout = xChip; + yReadout = yChip; + } - pmSourceDefinePixels (source, readout, xCell, yCell, radius); + pmSourceDefinePixels (source, readout, xReadout, yReadout, radius); return (source); Index: branches/eam_branches/20090522/psModules/test/objects/tap_pmSource.c =================================================================== --- branches/eam_branches/20090522/psModules/test/objects/tap_pmSource.c (revision 24238) +++ branches/eam_branches/20090522/psModules/test/objects/tap_pmSource.c (revision 24529) @@ -432,5 +432,5 @@ psArray *sources = psArrayAlloc(NUM_SOURCES); psMetadata *recipe = psMetadataAlloc(); - bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0); + bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0); @@ -452,5 +452,5 @@ psArray *sources = psArrayAlloc(NUM_SOURCES); psMetadata *recipe = psMetadataAlloc(); - bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0); + bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0); @@ -490,5 +490,5 @@ } psMetadata *recipe = psMetadataAlloc(); - bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0); + bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0); rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0); Index: branches/eam_branches/20090522/psphot/doc/notes.20090523.txt =================================================================== --- branches/eam_branches/20090522/psphot/doc/notes.20090523.txt (revision 24529) +++ branches/eam_branches/20090522/psphot/doc/notes.20090523.txt (revision 24529) @@ -0,0 +1,41 @@ + +20090523 : extended source analysis -- further development + + issues to address: + + * quality of moments used for initial psf selection + - what is the behavior of moments with Gaussian windowing as function of S/N & size? + - can we choose a sigma value that will give reliable moments for an expected range of seeing? + + * moments for galaxy model guesses + - trade-off between Gaussian windowing, tophat windowing, and S/N windowing in footprint area + + initial dev work: + + * add Gaussian windowing, S/N windowing, (tophat windowing) and options to select + + * check visualization of moments + +20090525 : some clarity of issues: + + I've made some progress understanding why the moments analysis was + giving poor results. I've added Gaussian (pseudo-gaussian) + windowing and the ability to set the per-pixel S/N limit. The first + important conclusion was: + + * the per-pixel S/N limit of 1.0 was driving the Mxx,Myy values low + for faint sources. This should be completely disabled. + + The next point is about the Gaussian vs tophat windowing. One + interesting realization was that the Gaussian windowing (or tophat + for that matter) have a stronger impact on the measured moments for + sources with profiles that are non-Gaussian. + + * Gauss with seeing (FWHM) = 1.0 arcsec = 4pix -> sigma = 1.702 + + Gauss window | measured Mxx + 2.0 1.75 + 4.0 2.58 + 6.0 2.85 + + * QGAUSS with see Index: branches/eam_branches/20090522/psphot/doc/stack.txt =================================================================== --- branches/eam_branches/20090522/psphot/doc/stack.txt (revision 24529) +++ branches/eam_branches/20090522/psphot/doc/stack.txt (revision 24529) @@ -0,0 +1,38 @@ + +20090606 : design notes on the multi-image photometry analysis + + assumptions: + + * each input image represents the same sky pixels : they are warped to a common frame. + + * each image has been previously processed, with the background subtracted and the psf model determined + + basic outline: + + * load images, masks, variance : for now, we should load all images. in the future, we might be able to split this into segments. + + * option: model and subtract background for each image (should not be needed) + + * option: smooth each image with a psf or fraction (not needed if stack smooths too much?) + + * generate the significance (chisq) image: X_i = sum(f_i^2 / var_i) + + * perform the peak detection on X_i + + * for each of the input images: + + * generate the significance image S_i = f_i^2 / var_i + + * perform the peak detection on S_i + + * merge the list of peaks + + * generate the footprints + + * linear fit to the peaks with the set of images (psf for each image) + + + design issues: + + * pmSource represents the analysis for an object on a single image + * extended pmSourceSet to group connected pmSources on multiple images? Index: branches/eam_branches/20090522/psphot/src/Makefile.am =================================================================== --- branches/eam_branches/20090522/psphot/src/Makefile.am (revision 24238) +++ branches/eam_branches/20090522/psphot/src/Makefile.am (revision 24529) @@ -16,14 +16,14 @@ # Force recompilation of psphotVersion.c, since it gets the version information -psphotVersion.c: psphotVersionDefinitions.h -psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE - -$(RM) psphotVersionDefinitions.h - $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h -FORCE: ; +# psphotVersion.c: psphotVersionDefinitions.h +# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE +# -$(RM) psphotVersionDefinitions.h +# $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h +# FORCE: ; libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotTest +bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -34,4 +34,8 @@ psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotTest_LDADD = libpsphot.la + +psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotMomentsStudy_LDADD = libpsphot.la psphot_SOURCES = \ @@ -54,4 +58,7 @@ psphotTest_SOURCES = \ psphotTest.c + +psphotMomentsStudy_SOURCES = \ + psphotMomentsStudy.c libpsphot_la_SOURCES = \ Index: branches/eam_branches/20090522/psphot/src/psphotArguments.c =================================================================== --- branches/eam_branches/20090522/psphot/src/psphotArguments.c (revision 24238) +++ branches/eam_branches/20090522/psphot/src/psphotArguments.c (revision 24529) @@ -105,5 +105,5 @@ pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE); if (config == NULL) { - psError(PSPHOT_ERR_CONFIG, false, "Can't read site configuration"); + psErrorStackPrint(stderr, "Can't read site configuration"); exit(PS_EXIT_CONFIG_ERROR); } @@ -121,6 +121,5 @@ if ((N = psArgumentGet (argc, argv, "-modeltest"))) { if (argc<=N+2) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 2 more arguments; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 2 more arguments; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } @@ -137,6 +136,5 @@ if ((N = psArgumentGet (argc, argv, "-model"))) { if (argc<=N+1) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } @@ -149,6 +147,5 @@ if ((N = psArgumentGet (argc, argv, "-fitmode"))) { if (argc<=N+1) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } @@ -159,6 +156,5 @@ if ((N = psArgumentGet (argc, argv, "-fitset"))) { if (argc<=N+1) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } @@ -172,6 +168,5 @@ if ((N = psArgumentGet (argc, argv, "-photcode"))) { if (argc<=N+1) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } @@ -190,6 +185,5 @@ if ((N = psArgumentGet (argc, argv, "-break"))) { if (argc<=N+1) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } @@ -202,6 +196,5 @@ if ((N = psArgumentGet (argc, argv, "-region"))) { if (argc<=N+1) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } @@ -214,6 +207,5 @@ if ((N = psArgumentGet (argc, argv, "-chip"))) { if (argc<=N+1) { - psError(PSPHOT_ERR_ARGUMENTS, true, - "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); exit(PS_EXIT_CONFIG_ERROR); } Index: branches/eam_branches/20090522/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- branches/eam_branches/20090522/psphot/src/psphotFitSourcesLinearStack.c (revision 24529) +++ branches/eam_branches/20090522/psphot/src/psphotFitSourcesLinearStack.c (revision 24529) @@ -0,0 +1,252 @@ +# include "psphotInternal.h" + +// fit flux (and optionally sky model) to all reasonable sources +// with the linear fitting process. sources must have an associated +// model with selected pixels, and the fit radius must be defined + +// given the set of sources, each of which points to the pixels in the +// science image, we construct a set of simulated sources with their own pixels. +// these are used to determine the simultaneous linear fit of fluxes. +// the analysis is performed wrt the simulated pixel values + +static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal); + +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { + + bool status; + float x; + float y; + float f; + // float r; + + psTimerStart ("psphot.linear"); + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // source analysis is done in spatial order + sources = psArraySort (sources, pmSourceSortByY); + + // storage array for fitSources + psArray *fitSources = psArrayAllocEmpty (sources->n); + + bool CONSTANT_PHOTOMETRIC_WEIGHTS = + psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + if (!status) { + psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); + } + + // select the sources which will be used for the fitting analysis + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // turn this bit off and turn it on again if we pass this test + source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // do not include CRs in the full ensemble fit + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + + if (final) { + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; + } else { + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + } + + // generate model for sources without, or skip if we can't + if (!source->modelFlux) { + if (!pmSourceCacheModel (source, maskVal)) continue; + } + + source->mode |= PM_SOURCE_MODE_LINEAR_FIT; + psArrayAdd (fitSources, 100, source); + } + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n); + + if (fitSources->n == 0) { + return true; + } + + // vectors to store stats for each object + // psVector *variance = psVectorAlloc (fitSources->n, PS_TYPE_F32); + psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32); + + // create the border matrix (includes the sparse matrix) + // for just sky: 1 row; for x,y terms: 3 rows + psSparse *sparse = psSparseAlloc (fitSources->n, 100); + psSparseBorder *border = psSparseBorderAlloc (sparse, 1); + + // fill out the sparse matrix elements and border elements (B) + // SRCi is the current source of interest + // SRCj is a possibly overlapping source + for (int i = 0; i < fitSources->n; i++) { + pmSource *SRCi = fitSources->data[i]; + + // diagonal elements of the sparse matrix (auto-cross-product) + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseMatrixElement (sparse, i, i, f); + + // the formal error depends on the weighting scheme + if (CONSTANT_PHOTOMETRIC_WEIGHTS) { + float var = pmSourceModelDotModel (SRCi, SRCi, false); + errors->data.F32[i] = 1.0 / sqrt(var); + } else { + errors->data.F32[i] = 1.0 / sqrt(f); + } + + + // find the image x model value + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseVectorElement (sparse, i, f); + + f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 0, f); + + // loop over all other stars following this one + for (int j = i + 1; j < fitSources->n; j++) { + pmSource *SRCj = fitSources->data[j]; + + // skip over disjoint source images, break after last possible overlap + if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) break; + if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue; + if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) continue; + if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue; + + // got an overlap; calculate cross-product and add to output array + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseMatrixElement (sparse, j, i, f); + } + } + + psSparseResort (sparse); + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); + + // set the sky, sky_x, sky_y components of border matrix + SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal); + + psSparseConstraint constraint; + constraint.paramMin = 0.0; + constraint.paramMax = 1e8; + constraint.paramDelta = 1e8; + + // solve for normalization terms (need include local sky?) + psVector *norm = psSparseSolve (NULL, constraint, sparse, 5); + + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); + + // XXXX **** philosophical question: + // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit: + // should retain the chisq and errors from the intermediate non-linear fit? + // the non-linear fit provides better values for the position errors, and for + // extended sources, the shape errors + + // adjust I0 for fitSources and subtract + for (int i = 0; i < fitSources->n; i++) { + pmSource *source = fitSources->data[i]; + pmModel *model = pmSourceGetModel (NULL, source); + + // assign linearly-fitted normalization + if (isnan(norm->data.F32[i])) { + psAbort("linear fitted source is nan"); + } + + model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; + model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; + // XXX is the value of 'errors' modified by the sky fit? + + // subtract object + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + } + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); + + // measure chisq for each source + for (int i = 0; final && (i < fitSources->n); i++) { + pmSource *source = fitSources->data[i]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; + pmModel *model = pmSourceGetModel (NULL, source); + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal); + } + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); + + // psFree (index); + psFree (sparse); + psFree (fitSources); + psFree (norm); + psFree (errors); + psFree (border); + + psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear")); + + psphotVisualShowResidualImage (readout); + psphotVisualShowFlags (sources); + + return true; +} + +// Calculate the weight terms for the sky fit component of the matrix. This function operates +// on the pixels which correspond to all of the sources of interest. These elements fill in +// the border matrix components in the sparse matrix equation. +static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal) { + + // generate the image-wide weight terms + // turn on MARK for all image pixels + psRegion fullArray = psRegionSet (0, 0, 0, 0); + fullArray = psRegionForImage (readout->mask, fullArray); + psImageMaskRegion (readout->mask, fullArray, "OR", markVal); + + // turn off MARK for all object pixels + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmModel *model = pmSourceGetModel (NULL, source); + if (model == NULL) continue; + float x = model->params->data.F32[PM_PAR_XPOS]; + float y = model->params->data.F32[PM_PAR_YPOS]; + psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_IMAGE_MASK(markVal)); + } + + // accumulate the image statistics from the masked regions + psF32 **image = readout->image->data.F32; + psF32 **variance = readout->variance->data.F32; + psImageMaskType **mask = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; + + double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy; + w = x = y = x2 = xy = y2 = fo = fx = fy = 0; + + int col0 = readout->image->col0; + int row0 = readout->image->row0; + + for (int j = 0; j < readout->image->numRows; j++) { + for (int i = 0; i < readout->image->numCols; i++) { + if (mask[j][i]) continue; + if (constant_weights) { + wt = 1.0; + } else { + wt = variance[j][i]; + } + f = image[j][i]; + w += 1/wt; + fo += f/wt; + } + } + + // turn off MARK for all image pixels + psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // set the Border T elements + psSparseBorderElementG (border, 0, fo); + psSparseBorderElementT (border, 0, 0, w); + + return true; +} Index: branches/eam_branches/20090522/psphot/src/psphotModelTest.c =================================================================== --- branches/eam_branches/20090522/psphot/src/psphotModelTest.c (revision 24238) +++ branches/eam_branches/20090522/psphot/src/psphotModelTest.c (revision 24529) @@ -139,5 +139,5 @@ // get the source moments - status = pmSourceMoments (source, mRADIUS); + status = pmSourceMoments (source, mRADIUS, 0.0, 1.0); if (!status) psAbort("psSourceMoments error"); source->peak->value = source->moments->Peak; Index: branches/eam_branches/20090522/psphot/src/psphotMomentsStudy.c =================================================================== --- branches/eam_branches/20090522/psphot/src/psphotMomentsStudy.c (revision 24529) +++ branches/eam_branches/20090522/psphot/src/psphotMomentsStudy.c (revision 24529) @@ -0,0 +1,171 @@ +# include "psphotInternal.h" + +int SetUpKapa (Graphdata *graphdata); +int PlotVectors (int kapa, Graphdata *graphdata, psVector *inVector, psVector *outVector); +int PlotModelSet (int kapa, Graphdata *graphdata, pmReadout *readout, char *name, float par7); + +FILE *output = NULL; + +int main (int argc, char **argv) { + + Graphdata graphdata; + int kapa; + + pmErrorRegister(); // register psModule's error codes/messages + pmModelClassInit(); + + // loop over sigma_window + // loop over sigma_input + // generate a fake source using the desired model + // add noise? + // measure the moments for the object + // accumulate Mxx, Myy + + // plot sigma_input vs sigma_output == sqrt(Mxx) + + // pmSourceMoments requires a source with a peak and pixels + + if (argc != 2) { + fprintf (stderr, "USAGE: psphotMomentsStudy (output)\n"); + exit (2); + } + + output = fopen (argv[1], "w"); + if (!output) exit (1); + fprintf (output, "# %4s %5s %6s %6s %6s\n", "type", "par7", "sigWin", "sigIn", "sigOut"); + + kapa = SetUpKapa(&graphdata); + + // create a containing image & associated readout + pmReadout *readout = pmReadoutAlloc(NULL); + readout->image = psImageAlloc(128, 128, PS_TYPE_F32); + + // create a dummy variance, but don't populate -- it is not used by pmSourceMoments if the sigma parameters is set to 0.0 + readout->variance = psImageAlloc(128, 128, PS_TYPE_F32); + + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_GAUSS", 0.0); + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PGAUSS", 0.0); + + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PS1_V1", 0.3); + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PS1_V1", 1.0); + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PS1_V1", 3.0); + + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_SERSIC", 0.125); + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_SERSIC", 0.250); + PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_SERSIC", 0.500); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + fclose (output); + + // psphotSaveImage(NULL, readout->image, "test.fits"); + exit (0); +} + +int SetUpKapa (Graphdata *graphdata) { + + int kapa = KapaOpenNamedSocket ("kapa", "psphot"); + + KapaClearPlots (kapa); + KapaInitGraph (graphdata); + KapaSetFont (kapa, "courier", 14); + + graphdata->color = KapaColorByName ("black"); + graphdata->xmin = -0.1; + graphdata->xmax = 1.6; + graphdata->ymin = -0.1; + graphdata->ymax = 1.6; + // graphdata->xmax = 10.1; + // graphdata->ymax = 10.1; + KapaSetLimits (kapa, graphdata); + + KapaBox (kapa, graphdata); + // KapaSendLabel (kapa, "&ss&h_in| (pixels)", KAPA_LABEL_XM); + // KapaSendLabel (kapa, "&ss&h_out| (pixels)", KAPA_LABEL_YM); + KapaSendLabel (kapa, "&ss&h_out| / &ss&h_window|", KAPA_LABEL_XM); + KapaSendLabel (kapa, "&ss&h_out| / &ss&h_in|", KAPA_LABEL_YM); + + return (kapa); +} + +int PlotVectors (int kapa, Graphdata *graphdata, psVector *inVector, psVector *outVector) { + + graphdata->color = KapaColorByName ("black"); + graphdata->ptype = 2; + graphdata->size = 1.0; + graphdata->style = 2; + KapaPrepPlot (kapa, inVector->n, graphdata); + KapaPlotVector (kapa, inVector->n, inVector->data.F32, "x"); + KapaPlotVector (kapa, inVector->n, outVector->data.F32, "y"); + + return 0; +} + + +int PlotModelSet (int kapa, Graphdata *graphdata, pmReadout *readout, char *name, float par7) { + + // create a model & associated source + // pmModelType type = pmModelClassGetType("PS_MODEL_GAUSS"); + + pmModelType type = pmModelClassGetType(name); + pmModel *model = pmModelAlloc(type); + + // set the model parameters + model->params->data.F32[PM_PAR_SKY] = 0.0; + model->params->data.F32[PM_PAR_I0] = 100.0; + model->params->data.F32[PM_PAR_XPOS] = 64.0; + model->params->data.F32[PM_PAR_YPOS] = 64.0; + + model->params->data.F32[PM_PAR_SXX] = 1.0 * M_SQRT2; + model->params->data.F32[PM_PAR_SYY] = 1.0 * M_SQRT2; + model->params->data.F32[PM_PAR_SXY] = 0.0; + + if (model->params->n > 7) { + model->params->data.F32[PM_PAR_7] = par7; + } + + pmSource *source = pmSourceFromModel(model, readout, 32.0, PM_SOURCE_TYPE_STAR); + + psVector *inVector = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *outVector = psVectorAllocEmpty(100, PS_TYPE_F32); + + for (float sigWindow = 1.0; sigWindow < 10.0; sigWindow += 1.0) { + + for (float sigIn = 0.5; sigIn < 10.0; sigIn += 0.1) { + + model->params->data.F32[PM_PAR_SXX] = sigIn * M_SQRT2; + model->params->data.F32[PM_PAR_SYY] = sigIn * M_SQRT2; + + // generate the modelFlux + pmSourceCacheModel(source, 0); + + // instantiate the source + pmSourceAdd(source, PM_MODEL_OP_FUNC, 0); + + pmSourceMoments (source, 32.0, sigWindow, 0.0); + // fprintf (stderr, "sigOut : %f\n", sqrt(source->moments->Mxx)); + + psVectorAppend(inVector, sqrt(source->moments->Mxx) / sigWindow); + psVectorAppend(outVector, sqrt(source->moments->Mxx) / sigIn); + + psImageInit (readout->image, 0.0); + + fprintf (output, " %4d %5.2f %6.3f %6.3f %6.3f\n", type, par7, sigWindow, sigIn, sqrt(source->moments->Mxx)); + } + } + + PlotVectors (kapa, graphdata, inVector, outVector); + inVector->n = 0; + outVector->n = 0; + + psFree (source); + psFree (inVector); + psFree (outVector); + + return 1; +} Index: branches/eam_branches/20090522/psphot/src/psphotSourceFits.c =================================================================== --- branches/eam_branches/20090522/psphot/src/psphotSourceFits.c (revision 24238) +++ branches/eam_branches/20090522/psphot/src/psphotSourceFits.c (revision 24529) @@ -224,5 +224,6 @@ // recalculate the source moments using the larger extended-source moments radius - if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false; + // at this stage, skip Gaussian windowing, and allow pixels with S/N > 0.5 + if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.5)) return false; psTrace ("psphot", 5, "trying blob...\n"); Index: branches/eam_branches/20090522/psphot/src/psphotSourceStats.c =================================================================== --- branches/eam_branches/20090522/psphot/src/psphotSourceStats.c (revision 24238) +++ branches/eam_branches/20090522/psphot/src/psphotSourceStats.c (revision 24529) @@ -89,25 +89,14 @@ if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS"); psFree (job); return NULL; } psFree(job); - -# if (0) - int nfail = 0; - int nmoments = 0; - if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - Nfail += nfail; - Nmoments += nmoments; -# endif } // wait for the threads to finish and manage results if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS"); return NULL; } @@ -149,9 +138,13 @@ psMetadata *recipe = job->args->data[1]; - float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); - if (!status) return false; - float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); - if (!status) return false; - float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); + if (!status) return false; + float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); + if (!status) return false; + float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + if (!status) return false; + float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN"); + if (!status) return false; + float SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA"); if (!status) return false; @@ -202,5 +195,5 @@ // measure basic source moments - status = pmSourceMoments (source, RADIUS); + status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN); if (status) { Nmoments ++; @@ -212,5 +205,5 @@ BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS); + status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN); if (status) { source->mode |= PM_SOURCE_MODE_BIG_RADIUS; @@ -292,5 +285,5 @@ // measure basic source moments - status = pmSourceMoments (source, RADIUS); + status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN); if (status) { Nmoments ++; @@ -302,5 +295,5 @@ BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS); + status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN); if (status) { Nmoments ++; Index: branches/eam_branches/20090522/psphot/src/psphotVisual.c =================================================================== --- branches/eam_branches/20090522/psphot/src/psphotVisual.c (revision 24238) +++ branches/eam_branches/20090522/psphot/src/psphotVisual.c (revision 24529) @@ -471,6 +471,7 @@ bool status; Graphdata graphdata; - - if (!pmVisualIsVisual()) return true; + KapaSection section; + + // if (!pmVisualIsVisual()) return true; if (kapa3 == -1) { @@ -485,27 +486,43 @@ KapaClearPlots (kapa3); KapaInitGraph (&graphdata); - - // there are N regions: use the first (guaranteed to exist) to get the overall limits - psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000"); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); - float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); - - // examine sources to set data range - graphdata.xmin = -0.05; - graphdata.ymin = -0.05; - graphdata.xmax = 2.0*psfX; - graphdata.ymax = 2.0*psfY; - KapaSetLimits (kapa3, &graphdata); - - KapaSetFont (kapa3, "helvetica", 14); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM); - KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM); - + KapaSetFont (kapa3, "courier", 14); + + float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); + + // select the max psfX,Y values for the plot limits + float Xmin = 0.0, Xmax = 0.0; + float Ymin = 0.0, Ymax = 0.0; + { + int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + + float X0 = psfX - 4.0*psfdX; + float X1 = psfX + 4.0*psfdX; + float Y0 = psfY - 4.0*psfdY; + float Y1 = psfY + 4.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } + } + } + + // storage vectors for data to be plotted psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32); psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32); psVector *xFaint = psVectorAllocEmpty (sources->n, PS_TYPE_F32); psVector *yFaint = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mFaint = psVectorAllocEmpty (sources->n, PS_TYPE_F32); // construct the vectors @@ -519,19 +536,45 @@ xFaint->data.F32[nF] = source->moments->Mxx; yFaint->data.F32[nF] = source->moments->Myy; + mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum); nF++; // XXX make this a user-defined cutoff - if (source->moments->SN < 50) + if (source->moments->SN < SN_LIM) continue; xBright->data.F32[nB] = source->moments->Mxx; yBright->data.F32[nB] = source->moments->Myy; + mBright->data.F32[nB] = -2.5*log10(source->moments->Sum); nB++; } xFaint->n = nF; yFaint->n = nF; + mFaint->n = nF; xBright->n = nB; yBright->n = nB; + mBright->n = nB; + + // three sections: MxxMyy, MagMxx, MagMyy + + // first section: MxxMyy + section.dx = 0.75; + section.dy = 0.75; + section.x = 0.00; + section.y = 0.00; + section.name = psStringCopy ("MxxMyy"); + KapaSetSection (kapa3, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (kapa3, &graphdata); + + KapaBox (kapa3, &graphdata); + KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("black"); @@ -551,6 +594,82 @@ KapaPlotVector (kapa3, nB, yBright->data.F32, "y"); - // XXX draw N circles to outline the clumps + // second section: MagMyy + section.dx = 0.75; + section.dy = 0.25; + section.x = 0.00; + section.y = 0.80; + section.name = psStringCopy ("MagMyy"); + KapaSetSection (kapa3, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -17.1; + graphdata.xmax = -7.9; + graphdata.ymin = Ymin; + graphdata.ymax = Ymax; + KapaSetLimits (kapa3, &graphdata); + + strcpy (graphdata.labels, "0210"); + KapaBox (kapa3, &graphdata); + KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP); + KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapa3, nF, &graphdata); + KapaPlotVector (kapa3, nF, mFaint->data.F32, "x"); + KapaPlotVector (kapa3, nF, yFaint->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapa3, nB, &graphdata); + KapaPlotVector (kapa3, nB, mBright->data.F32, "x"); + KapaPlotVector (kapa3, nB, yBright->data.F32, "y"); + + // third section: MagMxx + section.dx = 0.25; + section.dy = 0.75; + section.x = 0.80; + section.y = 0.00; + section.name = psStringCopy ("MagMxx"); + KapaSetSection (kapa3, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.xmax = Xmax; + graphdata.ymin = -7.9; + graphdata.ymax = -17.1; + KapaSetLimits (kapa3, &graphdata); + + strcpy (graphdata.labels, "2001"); + KapaBox (kapa3, &graphdata); + KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapa3, nF, &graphdata); + KapaPlotVector (kapa3, nF, xFaint->data.F32, "x"); + KapaPlotVector (kapa3, nF, mFaint->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapa3, nB, &graphdata); + KapaPlotVector (kapa3, nB, xBright->data.F32, "x"); + KapaPlotVector (kapa3, nB, mBright->data.F32, "y"); + + // draw N circles to outline the clumps { + KapaSelectSection (kapa3, "MxxMyy"); + // draw a circle centered on psfX,Y with size of the psf limit psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); @@ -562,4 +681,10 @@ graphdata.color = KapaColorByName ("blue"); graphdata.style = 0; + + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (kapa3, &graphdata); for (int n = 0; n < nRegions; n++) { @@ -626,6 +751,8 @@ psFree (xBright); psFree (yBright); + psFree (mBright); psFree (xFaint); psFree (yFaint); + psFree (mFaint); // pause and wait for user input: