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);