Index: /branches/eam_branches/ipp-20110404/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/camera/pmFPAMaskWeight.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/camera/pmFPAMaskWeight.c (revision 31361) @@ -587,5 +587,7 @@ double imageValue, varianceValue; // Image and variance value from interpolation psImageMaskType maskValue = 0; // Mask value from interpolation - psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x, y, interp); + + // interpolate to pixel center (index + 0.5) + psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x + 0.5, y + 0.5, interp); if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) { psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y); Index: /branches/eam_branches/ipp-20110404/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/imcombine/pmPSFEnvelope.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/imcombine/pmPSFEnvelope.c (revision 31361) @@ -380,5 +380,5 @@ // measure the source moments: tophat windowing, no pixel S/N cutoff - if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, maskVal)) { + if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, 0.0, maskVal)) { // Can't do anything about it; limp along as best we can psErrorClear(); Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/Makefile.am (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/Makefile.am (revision 31361) @@ -114,4 +114,5 @@ pmTrend2D.h \ pmGrowthCurve.h \ + pmGrowthCurveGenerate.h \ pmSourceMatch.h \ pmDetEff.h \ Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c (revision 31361) @@ -81,10 +81,10 @@ while (radius < Rlin) { // fprintf (stderr, "r: %f\n", radius); - psVectorAppend (growth->radius, radius); + psVectorAppend (growth->radius, radius - 0.001); radius += 1.0; } while (radius < maxRadius) { // fprintf (stderr, "r: %f\n", radius); - psVectorAppend (growth->radius, radius); + psVectorAppend (growth->radius, radius - 0.001); radius *= fR; radius = (int) (radius + 0.5); Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c (revision 31361) @@ -50,6 +50,5 @@ #include "pmSourcePhotometry.h" - -pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc); +#include "pmGrowthCurveGenerate.h" /*****************************************************************************/ @@ -226,2 +225,188 @@ return growth; } + +// we generate the growth curve for the center of the image with the specified psf model +bool pmGrowthCurveGenerateFromSources (pmReadout *readout, pmPSF *psf, psArray *sources, bool INTERPOLATE_AP, psImageMaskType maskVal, psImageMaskType markVal) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(readout->image, false); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + pmSourcePhotometryMode photMode = INTERPOLATE_AP ? PM_SOURCE_PHOT_INTERP : 0; + + // measure the growth curve for each PSF source and average them together + psArray *growths = psArrayAllocEmpty (100); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + pmGrowthCurve *growth = pmGrowthCurveForSource (source, psf, photMode, maskVal, markVal); + if (!growth) continue; + + psArrayAdd (growths, 100, growth); + psFree (growth); + } + psAssert (growths->n, "cannot build growth curve (no valid PSF stars?)"); + +# if (0) + FILE *f = fopen ("growth.dat", "w"); + assert (f); + + for (int j = 0; j < psf->growth->radius->n; j++) { + + fprintf (f, "%f : ", psf->growth->radius->data.F32[j]); + + // XXX dump the growth curves: + for (int i = 0; i < growths->n; i++) { + + pmGrowthCurve *growth = growths->data[i]; + if (!growth) continue; + + fprintf (f, "%8.4f ", growth->apMag->data.F32[j]); + } + fprintf (f, "\n"); + } + fclose (f); +# endif + + // just use a simple sample median to get the 'best' value from each growth curve... + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + psVector *values = psVectorAlloc (growths->n, PS_DATA_F32); + + // loop over a range of source fluxes + // no need to interpolate since we have forced the object center + // to 0.5, 0.5 above + for (int i = 0; i < psf->growth->radius->n; i++) { + + // median the values for each radial bin + values->n = 0; + for (int j = 0; j < growths->n; j++) { + pmGrowthCurve *growth = growths->data[j]; + if (!isfinite(growth->apMag->data.F32[i])) continue; + psVectorAppend (values, growth->apMag->data.F32[i] - growth->fitMag); + } + if (values->n == 0) { + psf->growth->apMag->data.F32[i] = NAN; + } else { + if (!psVectorStats (stats, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + psf->growth->apMag->data.F32[i] = stats->sampleMedian; + } + } + + psf->growth->fitMag = psf->growth->apMag->data.F32[psf->growth->radius->n-1]; + psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius); + psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef; + + psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f (fitMag - apMag @ ref : %f - %f)\n", psf->growth->apLoss, psf->growth->fitMag, psf->growth->apRef); + + psFree (growths); + psFree (stats); + psFree (values); + + return true; +} + +pmGrowthCurve *pmGrowthCurveForSource (pmSource *source, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal, psImageMaskType markVal) { + + float radius; + + assert (psf->growth); + + float minRadius = psf->growth->radius->data.F32[0]; + pmGrowthCurve *growth = pmGrowthCurveAlloc (minRadius, psf->growth->maxRadius, psf->growth->refRadius); + + // measure the fitMag for this source (for normalization) + // pmSourcePhotometryModel (&fitMag, NULL, source->psfModel); + growth->fitMag = source->psfMag; + + float xc = source->peak->xf; + float yc = source->peak->yf; + + // Loop over the range of radii + for (int i = 0; i < growth->radius->n; i++) { + + radius = growth->radius->data.F32[i]; + + // mask the given aperture and measure the apMag + psImageKeepCircle (source->maskObj, xc, yc, radius, "OR", markVal); + + if (!pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, radius)) { + psFree (growth); + return NULL; + } + + // if (!pmSourcePhotometryAper (NULL, &apMag, NULL, NULL, NULL, source->pixels, NULL, source->maskObj, maskVal)) { + // psFree (growth); + // return NULL; + // } + + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + + growth->apMag->data.F32[i] = source->apMag; + } + return growth; +} + +# if (0) + // solve for the best growth and the offsets per star to minimize the scatter + + // generate the per-star and per-radius sums + int Nradius = psf->growth->radius->n; + int Nstar = growths->n; + + psVector *Mstar = psVectorAlloc (Nstar, PS_DATA_F32); + psVector *Mradius = psVectorAlloc (Nradius, PS_DATA_F32); + + psVectorInit (Mstar, 0.0); + psVectorInit (Mradius, 0.0); + + for (int i = 0; i < Nstar; i++) { + pmGrowthCurve *growth = growths->data[i]; + for (int j = 0; j < Nradius; j++) { + Mstar->data.F32[i] += growth->apMag->data.F32[j]; + Mradius->data.F32[j] += growth->apMag->data.F32[j]; + } + } + + psImage *A = psImageAlloc(Nstar + Nradius, Nstar + Nradius, PS_DATA_F32); + psVector *B = psVectorAlloc(Nstar + Nradius, PS_DATA_F32); + + psImageInit (A, 0.0); + psVectorInit (B, 0.0); + + for (int i = 0; i < Nstar; i++) { + A->data.F32[i][i] = Nradius; + B->data.F32[i] = Mstar->data.F32[i]; + } + for (int i = 0; i < Nradius; i++) { + int j = i + Nstar; + A->data.F32[j][j] = Nstar; + B->data.F32[j] = Mradius->data.F32[i]; + } + + if (!psMatrixGJSolve(A, B)) { + psAbort("GJ failure"); + } + + for (int i = 0; i < Nradius; i++) { + psf->growth->apMag->data.F32[i] = B->data.F32[Nstar+i]; + } + + // XXX this analysis produces a biased growth curve: points with small radius are more + // likely to err on the side of being too faint (bacause of mis-aligned aperture), while on + // the bright side the tend to be too bright (because of contaminating neighbors). Bright + // stars are less likely to be biased, + + // XXX at this point, we could iterate and exclude some of the stars with low data quality + // XXX we could also weight the above by flux or something + +# endif Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.h =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.h (revision 31361) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.h (revision 31361) @@ -0,0 +1,22 @@ +/* @file pmGrowthCurveGenerate.h + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-24 00:11:02 $ + * Copyright 2007 IfA, University of Hawaii + */ + +# ifndef PM_GROWTH_CURVE_GENERATE_H +# define PM_GROWTH_CURVE_GENERATE_H + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +bool pmGrowthCurveGenerate (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType mark); +pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc); +bool pmGrowthCurveGenerateFromSources (pmReadout *readout, pmPSF *psf, psArray *sources, bool INTERPOLATE_AP, psImageMaskType maskVal, psImageMaskType markVal); +pmGrowthCurve *pmGrowthCurveForSource (pmSource *source, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal, psImageMaskType markVal); + +/// @} +# endif /* PM_GROWTH_CURVE_GENERATE_H */ Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.h =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.h (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.h (revision 31361) @@ -94,5 +94,4 @@ double pmPSF_SXYtoModel (psF32 *fittedPar); -bool pmGrowthCurveGenerate (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType mark); pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...); Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF_IO.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF_IO.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF_IO.c (revision 31361) @@ -65,4 +65,7 @@ bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file); bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf); + +bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file); +bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf); bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file) @@ -600,4 +603,9 @@ if (!pmPSFmodelWrite_ApTrend(file, psf)) { psError(psErrorCodeLast(), false, "Unable to write PSF ApTrend"); + return false; + } + + if (!pmPSFmodelWrite_GrowthCurve(file, psf)) { + psError(psErrorCodeLast(), false, "Unable to write PSF Growth Curve"); return false; } @@ -1104,4 +1112,9 @@ } + if (!pmPSFmodelRead_GrowthCurve(psf, file)) { + psError(psErrorCodeLast(), false, "Unable to read PSF Growth Curve"); + return false; + } + psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); psFree (psf); @@ -1230,4 +1243,101 @@ } +// write aperture trend to a FITS table +bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf) { + + pmGrowthCurve *growth = psf->growth; + if (growth == NULL) { + psWarning ("no PSF Growth Curve to write out, skipping"); + return true; + } + + // we need to write a header for the table, + psMetadata *header = psMetadataAlloc(); + + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MIN_RAD", 0, "", growth->radius->data.F32[0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MAX_RAD", 0, "", growth->maxRadius); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_REF_RAD", 0, "", growth->refRadius); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_LOSS", 0, "", growth->apLoss); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_REF", 0, "", growth->apRef); + psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_FIT_MAG", 0, "", growth->fitMag); + + // build a FITS table of the ApTrend (only 1) + psArray *table = psArrayAllocEmpty (100); + for (int i = 0; i < growth->apMag->n; i++) { + psMetadata *row = psMetadataAlloc (); + psMetadataAddF32 (row, PS_LIST_TAIL, "RADIUS", 0, "", growth->radius->data.F32[i]); + psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG", 0, "", growth->apMag->data.F32[i]); + psArrayAdd (table, 100, row); + psFree (row); + } + + // write an empty FITS segment if we have no PSF information + if (table->n == 0) { + psError(PM_ERR_PROG, true, "No PSF data to write."); + psFree(table); + psFree(header); + return false; + } + + psTrace ("pmFPAfile", 5, "writing psf Growth Curve data %s\n", "GROWTH_CURVE"); + if (!psFitsWriteTable(file->fits, header, table, "GROWTH_CURVE")) { + psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "GROWTH_CURVE"); + psFree(table); + psFree(header); + return false; + } + + psFree (table); + psFree (header); + return true; +} + +// read aperture trend to a FITS table +bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file) { + + bool status; + + // move fits pointer to AP_TREND section + // advance to the table data extension + if (!psFitsMoveExtNameClean (file->fits, "GROWTH_CURVE")) { + psWarning ("no Growth Curve data in PSF file, skipping"); + return true; + } + + psMetadata *header = psFitsReadHeader (NULL, file->fits); + if (!header) { + psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE header."); + return false; + } + + // read the raw table data + psArray *table = psFitsReadTable (file->fits); + if (!table) { + psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE table."); + psFree(header); + return false; + } + + float minRadius = psMetadataLookupF32 (&status, header, "GROWTH_MIN_RAD"); if (!status) return false; + float maxRadius = psMetadataLookupF32 (&status, header, "GROWTH_MAX_RAD"); if (!status) return false; + float refRadius = psMetadataLookupF32 (&status, header, "GROWTH_REF_RAD"); if (!status) return false; + + psf->growth = pmGrowthCurveAlloc(minRadius, maxRadius, refRadius); + + psf->growth->apLoss = psMetadataLookupF32 (&status, header, "GROWTH_AP_LOSS"); if (!status) return false; + psf->growth->apRef = psMetadataLookupF32 (&status, header, "GROWTH_AP_REF"); if (!status) return false; + psf->growth->fitMag = psMetadataLookupF32 (&status, header, "GROWTH_FIT_MAG"); if (!status) return false; + + // fill in the matching psf->params entries + for (int i = 0; i < table->n; i++) { + psMetadata *row = table->data[i]; + psf->growth->apMag->data.F32[i] = psMetadataLookupF32 (&status, row, "AP_MAG"); if (!status) return false; + } + + psFree (header); + psFree (table); + return true; +} + bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header) { Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.h =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.h (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.h (revision 31361) @@ -82,8 +82,8 @@ psArray *blends; ///< collection of sources thought to be confused with object float psfMag; ///< calculated from flux in modelPSF + float psfMagErr; ///< error in psfMag float psfFlux; ///< calculated from flux in modelPSF - float psfFluxErr; ///< calculated from flux in modelPSF - float extMag; ///< calculated from flux in modelEXT - float errMag; ///< error in psfMag OR extMag (depending on type) + float psfFluxErr; ///< error in psfFlux + float extMag; ///< calculated from flux in modelEXT -- NOTE this is not actually used float apMag; ///< apMag corresponding to psfMag or extMag (depending on type) float apMagRaw; ///< raw mag in given aperture @@ -254,4 +254,5 @@ float sigma, ///< size of Gaussian window function (<= 0.0 -> skip window) float minSN, ///< minimum pixel significance + float minKronRadius, ///< minimum pixel significance psImageMaskType maskVal ); @@ -273,4 +274,6 @@ float xGuess, float yGuess); +float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM); + pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source); Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31361) @@ -116,5 +116,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); @@ -161,6 +161,8 @@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); + + // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.KronCore); - // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR", PS_DATA_F32, "Kron Error (in 1.0 R1)", moments.KronCoreErr); + // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR", PS_DATA_F32, "Kron Error (in 1.0 R1)", moments.KronCoreErr); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); @@ -270,10 +272,10 @@ // XXX use these to determine PAR[PM_PAR_I0]? source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->errMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN; + dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; pmPSF_AxesToModel (PAR, axes); @@ -292,6 +294,6 @@ // we no longer sort by S/N, only flux - // if (isfinite (source->errMag) && (source->errMag > 0.0)) { - // source->peak->SN = 1.0 / source->errMag; + // if (isfinite (source->psfMagErr) && (source->psfMagErr > 0.0)) { + // source->peak->SN = 1.0 / source->psfMagErr; // } else { // source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceMoments.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceMoments.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceMoments.c (revision 31361) @@ -66,5 +66,5 @@ // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid) -bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) +bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -74,5 +74,5 @@ // this function assumes the sky has been well-subtracted for the image - psF32 sky = 0.0; + float sky = 0.0; if (source->moments == NULL) { @@ -80,11 +80,11 @@ } - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 SumCore = 0.0; - psF32 VarCore = 0.0; - psF32 R2 = PS_SQR(radius); - psF32 minSN2 = PS_SQR(minSN); - psF32 rsigma2 = 0.5 / PS_SQR(sigma); + float Sum = 0.0; + float Var = 0.0; + float SumCore = 0.0; + float VarCore = 0.0; + float R2 = PS_SQR(radius); + float minSN2 = PS_SQR(minSN); + float rsigma2 = 0.5 / PS_SQR(sigma); // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -109,22 +109,22 @@ // Xn = SUM (x - xc)^n * (z - sky) - psF32 RFW = 0.0; - psF32 RHW = 0.0; - - psF32 RF = 0.0; - psF32 RH = 0.0; - psF32 RS = 0.0; - psF32 XX = 0.0; - psF32 XY = 0.0; - psF32 YY = 0.0; - psF32 XXX = 0.0; - psF32 XXY = 0.0; - psF32 XYY = 0.0; - psF32 YYY = 0.0; - psF32 XXXX = 0.0; - psF32 XXXY = 0.0; - psF32 XXYY = 0.0; - psF32 XYYY = 0.0; - psF32 YYYY = 0.0; + float RFW = 0.0; + float RHW = 0.0; + + float RF = 0.0; + float RH = 0.0; + float RS = 0.0; + float XX = 0.0; + float XY = 0.0; + float YY = 0.0; + float XXX = 0.0; + float XXY = 0.0; + float XYY = 0.0; + float YYY = 0.0; + float XXXX = 0.0; + float XXXY = 0.0; + float XXYY = 0.0; + float XYYY = 0.0; + float YYYY = 0.0; Sum = 0.0; // the second pass may include slightly different pixels, re-determine Sum @@ -140,14 +140,14 @@ // center of mass in subimage. Note: the calculation below uses pixel index, so we correct // xCM, yCM from pixel coords to pixel index here. - psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage - psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage + float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage + float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage for (psS32 row = 0; row < source->pixels->numRows ; row++) { - psF32 yDiff = row - yCM; + float yDiff = row - yCM; if (fabs(yDiff) > radius) continue; - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; + float *vPix = source->pixels->data.F32[row]; + float *vWgt = source->variance->data.F32[row]; psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; @@ -164,14 +164,14 @@ if (isnan(*vPix)) continue; - psF32 xDiff = col - xCM; + float xDiff = col - xCM; if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); if (r2 > R2) continue; - psF32 fDiff = *vPix - sky; - psF32 pDiff = fDiff; - psF32 wDiff = *vWgt; + float fDiff = *vPix - sky; + float pDiff = fDiff; + float wDiff = *vWgt; // skip pixels below specified significance level. this is allowed, but should be @@ -184,7 +184,7 @@ // weighting over weights the sky for faint sources if (sigma > 0.0) { - psF32 z = r2 * rsigma2; + float z = r2 * rsigma2; assert (z >= 0.0); - psF32 weight = exp(-z); + float weight = exp(-z); wDiff *= weight; @@ -195,29 +195,29 @@ // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?) - psF32 r = sqrt(r2); - psF32 rf = r * fDiff; - psF32 rh = sqrt(r) * fDiff; - psF32 rs = fDiff; - - psF32 rfw = r * pDiff; - psF32 rhw = sqrt(r) * pDiff; - - psF32 x = xDiff * pDiff; - psF32 y = yDiff * pDiff; - - psF32 xx = xDiff * x; - psF32 xy = xDiff * y; - psF32 yy = yDiff * y; - - psF32 xxx = xDiff * xx / r; - psF32 xxy = xDiff * xy / r; - psF32 xyy = xDiff * yy / r; - psF32 yyy = yDiff * yy / r; - - psF32 xxxx = xDiff * xxx / r2; - psF32 xxxy = xDiff * xxy / r2; - psF32 xxyy = xDiff * xyy / r2; - psF32 xyyy = xDiff * yyy / r2; - psF32 yyyy = yDiff * yyy / r2; + float r = sqrt(r2); + float rf = r * fDiff; + float rh = sqrt(r) * fDiff; + float rs = fDiff; + + float rfw = r * pDiff; + float rhw = sqrt(r) * pDiff; + + float x = xDiff * pDiff; + float y = yDiff * pDiff; + + float xx = xDiff * x; + float xy = xDiff * y; + float yy = yDiff * y; + + float xxx = xDiff * xx / r; + float xxy = xDiff * xy / r; + float xyy = xDiff * yy / r; + float yyy = yDiff * yy / r; + + float xxxx = xDiff * xxx / r2; + float xxxy = xDiff * xxy / r2; + float xxyy = xDiff * xyy / r2; + float xyyy = xDiff * yyy / r2; + float yyyy = yDiff * yyy / r2; RF += rf; @@ -245,9 +245,6 @@ } - // if Mrf (first radial moment) is << sigma, we are getting into low-significance - // territory. saturate at 0.75*sigma. conversely, if Mrf is > radius, we are clearly - // making an error. saturate at radius. - source->moments->Mrf = MIN(radius, MAX(0.75*sigma, RF/RS)); - source->moments->Mrh = MIN(radius, MAX(0.75*sigma, RH/RS)); + source->moments->Mrf = RF/RS; + source->moments->Mrh = RH/RS; source->moments->Mxx = XX/Sum; @@ -266,12 +263,12 @@ source->moments->Myyyy = YYYY/Sum; - // XXX TEST: - // source->moments->KronFinner = RFW/Sum; - // source->moments->KronFouter = sigma; - - // Calculate the Kron magnitude (make this block optional?) - float radKinner = 1.0*source->moments->Mrf; - float radKron = 2.5*source->moments->Mrf; - float radKouter = 4.0*source->moments->Mrf; + // if Mrf (first radial moment) is very small, we are getting into low-significance + // territory. saturate at minKronRadius. conversely, if Mrf is > radius, we are clearly + // making an error. saturate at radius. + float kronRefRadius = MIN(radius, MAX(minKronRadius, source->moments->Mrf)); + + float radKinner = 1.0*kronRefRadius; + float radKron = 2.5*kronRefRadius; + float radKouter = 4.0*kronRefRadius; int nKronPix = 0; @@ -285,9 +282,9 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { - psF32 yDiff = row - yCM; + float yDiff = row - yCM; if (fabs(yDiff) > radKouter) continue; - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; + float *vPix = source->pixels->data.F32[row]; + float *vWgt = source->variance->data.F32[row]; psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; @@ -304,12 +301,12 @@ if (isnan(*vPix)) continue; - psF32 xDiff = col - xCM; + float xDiff = col - xCM; if (fabs(xDiff) > radKouter) continue; // radKron is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); - - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + + float pDiff = *vPix - sky; + float wDiff = *vWgt; // skip pixels below specified significance level. this is allowed, but should be @@ -318,5 +315,5 @@ if (PS_SQR(pDiff) < minSN2*wDiff) continue; - psF32 r = sqrt(r2); + float r = sqrt(r2); if (r < radKron) { Sum += pDiff; @@ -363,5 +360,5 @@ } -bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) { +bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { // First Pass: calculate the first moments (these are subtracted from the coordinates below) @@ -370,15 +367,15 @@ // .. etc - psF32 sky = 0.0; - - psF32 peakPixel = -PS_MAX_F32; + float sky = 0.0; + + float peakPixel = -PS_MAX_F32; psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 Var = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 R2 = PS_SQR(radius); - psF32 minSN2 = PS_SQR(minSN); - psF32 rsigma2 = 0.5 / PS_SQR(sigma); + float Sum = 0.0; + float Var = 0.0; + float X1 = 0.0; + float Y1 = 0.0; + float R2 = PS_SQR(radius); + float minSN2 = PS_SQR(minSN); + float rsigma2 = 0.5 / PS_SQR(sigma); float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage @@ -386,5 +383,5 @@ // we are guaranteed to have a valid pixel and variance at this location (right? right?) - // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); + // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]); // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel"); @@ -398,9 +395,9 @@ for (psS32 row = 0; row < source->pixels->numRows ; row++) { - psF32 yDiff = row + 0.5 - yPeak; + float yDiff = row + 0.5 - yPeak; if (fabs(yDiff) > radius) continue; - psF32 *vPix = source->pixels->data.F32[row]; - psF32 *vWgt = source->variance->data.F32[row]; + float *vPix = source->pixels->data.F32[row]; + float *vWgt = source->variance->data.F32[row]; psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row]; @@ -417,13 +414,13 @@ if (isnan(*vPix)) continue; - psF32 xDiff = col + 0.5 - xPeak; + float xDiff = col + 0.5 - xPeak; if (fabs(xDiff) > radius) continue; // radius is just a function of (xDiff, yDiff) - psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); if (r2 > R2) continue; - psF32 pDiff = *vPix - sky; - psF32 wDiff = *vWgt; + float pDiff = *vPix - sky; + float wDiff = *vWgt; // skip pixels below specified significance level. for a PSFs, this @@ -438,7 +435,7 @@ // weighting over weights the sky for faint sources if (sigma > 0.0) { - psF32 z = r2*rsigma2; + float z = r2*rsigma2; assert (z >= 0.0); - psF32 weight = exp(-z); + float weight = exp(-z); wDiff *= weight; @@ -449,6 +446,6 @@ Sum += pDiff; - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; + float xWght = xDiff * pDiff; + float yWght = yDiff * pDiff; X1 += xWght; @@ -507,2 +504,42 @@ return true; } + +float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) { + + psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32); + + for (int i = 0; i < sources->n; i++) { + pmSource *src = sources->data[i]; // Source of interest + if (!src || !src->moments) { + continue; + } + + if (src->mode & PM_SOURCE_MODE_BLEND) { + continue; + } + + if (!src->moments->nPixels) continue; + + if (src->moments->SN < PSF_SN_LIM) continue; + + // XXX put in Mxx,Myy cut based on clump location + + psVectorAppend(radii, src->moments->Mrf); + } + + // find the peak in this image + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + if (!psVectorStats (stats, radii, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n"); + psFree(stats); + return NAN; + } + + float minRadius = stats->sampleMedian; + + psFree(radii); + psFree(stats); + return minRadius; +} + Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourcePhotometry.c (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourcePhotometry.c (revision 31361) @@ -82,11 +82,8 @@ - apMag : only if S/N > AP_MIN_SN : is optionally corrected for curve-of-growth if: - - the source is a STAR (PSF) - the option is selected (mode & PM_SOURCE_PHOT_GROWTH) - psfMag : all sources with non-NULL modelPSF : is optionally corrected for aperture residual if: - - the source is a STAR (PSF) - the option is selected (mode & PM_SOURCE_PHOT_APCORR) - - extMag : all sources with non-NULL modelEXT **/ @@ -100,12 +97,10 @@ int status = false; - bool isPSF; float x, y; float SN; - pmModel *model; source->psfMag = NAN; source->extMag = NAN; - source->errMag = NAN; + source->psfMagErr = NAN; source->apMag = NAN; source->apMagRaw = NAN; @@ -113,38 +108,42 @@ source->apFluxErr = NAN; - // we must have a valid model - // XXX allow aperture magnitudes for sources without a model - model = pmSourceGetModel (&isPSF, source); - if (model == NULL) { - psTrace ("psModules.objects", 3, "fail mag : no valid model"); + // XXXXXX review: + // Select the 'best' model -- this is used for PSF_QF,_PERFECT & ???. isPSF is true if this + // object is a PSF (not extended). We must have a valid model. XXX NOTE: allow aperture + // magnitudes for sources without a model + + // select the psf model + pmModel *modelPSF = source->modelPSF; + if (modelPSF == NULL) { + psTrace ("psModules.objects", 3, "fail mag : no valid PSF model"); return false; } - // XXX handle negative flux, low-significance - if (model->dparams->data.F32[PM_PAR_I0] > 0) { - SN = fabs(model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0]); - source->errMag = 1.0 / SN; + // get the error on the PSF model magnitude + if (modelPSF->dparams->data.F32[PM_PAR_I0] > 0) { + SN = fabs(modelPSF->params->data.F32[PM_PAR_I0] / modelPSF->dparams->data.F32[PM_PAR_I0]); + source->psfMagErr = 1.0 / SN; } else { SN = NAN; - source->errMag = NAN; - } - x = model->params->data.F32[PM_PAR_XPOS]; - y = model->params->data.F32[PM_PAR_YPOS]; + source->psfMagErr = NAN; + } + // the source position is used to recenter the aperture for ap photometry + x = modelPSF->params->data.F32[PM_PAR_XPOS]; + y = modelPSF->params->data.F32[PM_PAR_YPOS]; // measure PSF model photometry - // XXX TEST: do not use flux scale - // XXX NOTE: turn this back on? - if (0 && psf->FluxScale) { - // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center - double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y); - psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y); - psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y); - source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]; - source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0]; - source->psfMag = -2.5*log10(source->psfFlux); - } else { - status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF); - source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]); - } + status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, modelPSF); + source->psfFluxErr = source->psfFlux * source->psfMagErr; + +# if (0) + // XXX NOTE: old code to use the flux scale. test & turn this back on? if so, need to save with psf model + // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center + double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y); + psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y); + psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y); + source->psfFlux = fluxScale * modelPSF->params->data.F32[PM_PAR_I0]; + source->psfFluxErr = fluxScale * modelPSF->dparams->data.F32[PM_PAR_I0]; + source->psfMag = -2.5*log10(source->psfFlux); +# endif if (mode == PM_SOURCE_PHOT_PSFONLY) { @@ -152,4 +151,5 @@ } + // get the EXT model photometry (all EXT models) // if we have a collection of model fits, check if one of them is a pointer to modelEXT if (source->modelFits) { @@ -172,7 +172,6 @@ } - // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS - // XXX add a flag for "ap_mag is corrected?" - if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) { + // Correct psfMag to match aperture magnitude system (NOTE : Growth curve is already applied to ApTrend) + if ((mode & PM_SOURCE_PHOT_APCORR) && psf && psf->ApTrend) { // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y); @@ -182,7 +181,7 @@ } - // measure the contribution of included pixels + // measure the contribution of included pixels to the PSF model fit if (mode & PM_SOURCE_PHOT_WEIGHT) { - pmSourcePixelWeight (source, model, source->maskObj, maskVal, radius); + pmSourcePixelWeight (source, modelPSF, source->maskObj, maskVal, radius); } @@ -218,6 +217,7 @@ y += dy; - if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, - NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) { + // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_LANCZOS2)) { + // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) { + if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BILINEAR)) { // Not much we can do about it psErrorClear(); @@ -233,5 +233,5 @@ // measure object aperture photometry - status = pmSourcePhotometryAperSource (source, model, flux, variance, mask, maskVal); + status = pmSourcePhotometryAperSource (source, modelPSF, flux, variance, mask, maskVal); if (!status) { psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag"); @@ -242,8 +242,10 @@ // detection limits (esp near bright neighbors) source->apMag = source->apMagRaw; - if (isfinite (source->apMag) && isPSF && psf) { + if (isfinite (source->apMag) && psf) { if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) { - source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius); - // XXX correct the apFlux? + float apOffset = pmGrowthCurveCorrect (psf->growth, source->apRadius); + source->apMag = source->apMagRaw + apOffset; + source->apFlux *= pow(10.0, -0.4*apOffset); + source->apFluxErr *= pow(10.0, -0.4*apOffset); } } Index: /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceUtils.h =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceUtils.h (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceUtils.h (revision 31361) @@ -38,3 +38,3 @@ /// @} -# endif /* PM_MODEL_CLASS_H */ +# endif /* PM_SOURCE_UTILS_H */ Index: /branches/eam_branches/ipp-20110404/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/ipp-20110404/psModules/src/psmodules.h (revision 31360) +++ /branches/eam_branches/ipp-20110404/psModules/src/psmodules.h (revision 31361) @@ -148,4 +148,5 @@ #include #include +#include #include #include