Index: /branches/eam_branches/20100225/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/Makefile.am (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/Makefile.am (revision 27494) @@ -22,4 +22,5 @@ pmSourceMasks.c \ pmSourceMoments.c \ + pmSourceDiffStats.c \ pmSourceExtendedPars.c \ pmSourceUtils.c \ @@ -40,4 +41,5 @@ pmSourceIO_CMF_PS1_V1.c \ pmSourceIO_CMF_PS1_V2.c \ + pmSourceIO_CMF_PS1_DV1.c \ pmSourceIO_MatchedRefs.c \ pmSourcePlots.c \ @@ -79,4 +81,5 @@ pmSource.h \ pmSourceMasks.h \ + pmSourceDiffStats.h \ pmSourceExtendedPars.h \ pmSourceUtils.h \ Index: /branches/eam_branches/20100225/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/pmGrowthCurveGenerate.c (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/pmGrowthCurveGenerate.c (revision 27494) @@ -157,5 +157,5 @@ // measure the fitMag for this model - pmSourcePhotometryModel (&fitMag, model); + pmSourcePhotometryModel (&fitMag, NULL, model); growth->fitMag = fitMag; Index: /branches/eam_branches/20100225/psModules/src/objects/pmSource.c =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/pmSource.c (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/pmSource.c (revision 27494) @@ -46,5 +46,5 @@ psFree(tmp->maskView); psFree(tmp->modelFlux); - psFree(tmp->psfFlux); + psFree(tmp->psfImage); psFree(tmp->moments); psFree(tmp->modelPSF); @@ -52,5 +52,6 @@ psFree(tmp->modelFits); psFree(tmp->extpars); - psFree(tmp->blends); + psFree(tmp->moments); + psFree(tmp->diffStats); psTrace("psModules.objects", 10, "---- end ----\n"); } @@ -68,5 +69,5 @@ psFree (source->maskView); psFree (source->modelFlux); - psFree (source->psfFlux); + psFree (source->psfImage); source->pixels = NULL; @@ -75,5 +76,5 @@ source->maskView = NULL; source->modelFlux = NULL; - source->psfFlux = NULL; + source->psfImage = NULL; return; } @@ -103,5 +104,5 @@ source->maskView = NULL; source->modelFlux = NULL; - source->psfFlux = NULL; + source->psfImage = NULL; source->moments = NULL; source->blends = NULL; @@ -113,9 +114,13 @@ source->tmpFlags = 0; source->extpars = NULL; + source->diffStats = NULL; + source->region = psRegionSet(NAN, NAN, NAN, NAN); psMemSetDeallocator(source, (psFreeFunc) sourceFree); // default values are NAN - source->psfMag = NAN; + source->psfMag = NAN; + source->psfFlux = NAN; + source->psfFluxErr = NAN; source->extMag = NAN; source->errMag = NAN; @@ -259,7 +264,7 @@ mySource->modelFlux = NULL; - // drop the old psfFlux pixels and force the user to re-create - psFree (mySource->psfFlux); - mySource->psfFlux = NULL; + // drop the old psfImage pixels and force the user to re-create + psFree (mySource->psfImage); + mySource->psfImage = NULL; } return extend; @@ -873,11 +878,11 @@ // if we already have a cached image, re-use that memory - source->psfFlux = psImageCopy (source->psfFlux, source->pixels, PS_TYPE_F32); - psImageInit (source->psfFlux, 0.0); + source->psfImage = psImageCopy (source->psfImage, source->pixels, PS_TYPE_F32); + psImageInit (source->psfImage, 0.0); // in some places (psphotEnsemble), we need a normalized version // in others, we just want the model. which is more commonly used? - // psfFlux always has unity normalization (I0 = 1.0) - pmModelAdd (source->psfFlux, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); + // psfImage always has unity normalization (I0 = 1.0) + pmModelAdd (source->psfImage, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); return true; } Index: /branches/eam_branches/20100225/psModules/src/objects/pmSource.h =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/pmSource.h (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/pmSource.h (revision 27494) @@ -16,4 +16,5 @@ #include "pmMoments.h" #include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" /// @addtogroup Objects Object Detection / Analysis Functions @@ -65,5 +66,5 @@ psImage *maskView; ///< view into global image mask for this object region psImage *modelFlux; ///< cached copy of the best model for this source - psImage *psfFlux; ///< cached copy of the psf model for this source + psImage *psfImage; ///< cached copy of the psf model for this source pmMoments *moments; ///< Basic moments measured for the object. pmModel *modelPSF; ///< PSF Model fit (parameters and type) @@ -75,4 +76,6 @@ psArray *blends; ///< collection of sources thought to be confused with object float psfMag; ///< calculated from flux in modelPSF + 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) @@ -86,4 +89,5 @@ psRegion region; ///< area on image covered by selected pixels pmSourceExtendedPars *extpars; ///< extended source parameters + pmSourceDiffStats *diffStats; ///< extra parameters for difference detections }; Index: /branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.c =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.c (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.c (revision 27494) @@ -539,4 +539,7 @@ status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname); } + if (!strcmp (exttype, "PS1_DV1")) { + status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname); + } if (xsrcname) { @@ -553,4 +556,7 @@ status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe); } + if (!strcmp (exttype, "PS1_DV1")) { + status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, sources, xsrcname, recipe); + } } if (xfitname) { @@ -567,4 +573,7 @@ status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname); } + if (!strcmp (exttype, "PS1_DV1")) { + status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, sources, xfitname); + } } psFree (outhead); @@ -1015,4 +1024,7 @@ sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header); } + if (!strcmp (exttype, "PS1_DV1")) { + sources = pmSourcesRead_CMF_PS1_DV1 (file->fits, hdu->header); + } if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) { Index: /branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.h =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.h (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/pmSourceIO.h (revision 27494) @@ -43,4 +43,8 @@ bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); +bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); +bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname); + bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config); @@ -53,4 +57,5 @@ psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header); bool pmSourcesWritePSFs (psArray *sources, char *filename); Index: /branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.c (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.c (revision 27494) @@ -109,7 +109,10 @@ 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->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]); + 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->modelPSF); + 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]); } @@ -119,5 +122,5 @@ for (int i = 0; i < source->modelFits->n; i++) { pmModel *model = source->modelFits->data[i]; - status = pmSourcePhotometryModel (&model->mag, model); + status = pmSourcePhotometryModel (&model->mag, NULL, model); if (model == source->modelEXT) foundEXT = true; } @@ -125,9 +128,9 @@ source->extMag = source->modelEXT->mag; } else { - status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); } } else { if (source->modelEXT) { - status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT); } } @@ -143,4 +146,9 @@ if (mode & PM_SOURCE_PHOT_WEIGHT) { pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal); + } + + // measure the contribution of included pixels + if (mode & PM_SOURCE_PHOT_DIFFSTATS) { + pmSourceMeasureDiffStats (source, maskVal); } @@ -217,21 +225,26 @@ // return source model magnitude -bool pmSourcePhotometryModel (float *fitMag, pmModel *model) -{ - PS_ASSERT_PTR_NON_NULL(fitMag, false); - if (model == NULL) { - return false; - } - - float fitSum = 0; - *fitMag = NAN; +bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model) +{ + psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)"); + if (model == NULL) return false; + + float mag = NAN; + float flux = NAN; // measure fitMag - fitSum = model->modelFlux (model->params); - if (fitSum <= 0) - return false; - if (!isfinite(fitSum)) - return false; - *fitMag = -2.5*log10(fitSum); + flux = model->modelFlux (model->params); + if (flux > 0) { + mag = -2.5*log10(flux); + } + if (fitMag) { + *fitMag = mag; + } + if (fitFlux) { + *fitFlux = flux; + } + + if (flux <= 0) return false; + if (!isfinite(flux)) return false; return (true); @@ -356,4 +369,55 @@ *pixWeight = validSum / modelSum; + return (true); +} + +# define FLUX_LIMIT 3.0 + +// return source aperture magnitude +bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + + if (source->diffStats == NULL) { + source->diffStats = pmSourceDiffStatsAlloc(); + } + + float fGood = 0.0; + float fBad = 0.0; + int nGood = 0; + int nMask = 0; + int nBad = 0; + + psImage *flux = source->pixels; + psImage *variance = source->variance; + psImage *mask = source->maskObj; + + for (int iy = 0; iy < flux->numRows; iy++) { + for (int ix = 0; ix < flux->numCols; ix++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) { + nMask ++; + continue; + } + + float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + + if (SN > +FLUX_LIMIT) { + nGood ++; + fGood += flux->data.F32[iy][ix]; + } + + if (SN < -FLUX_LIMIT) { + nBad ++; + fBad += flux->data.F32[iy][ix]; + } + } + } + + source->diffStats->fRatio = fGood / (fGood + fBad); + source->diffStats->nRatioBad = nGood / (nGood + nBad); + source->diffStats->nRatioMask = nGood / (nGood + nMask); + source->diffStats->nRatioAll = nGood / (nGood + nMask + nBad); + source->diffStats->nGood = nGood; + return (true); } Index: /branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.h (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/objects/pmSourcePhotometry.h (revision 27494) @@ -29,13 +29,15 @@ typedef enum { - PM_SOURCE_PHOT_NONE = 0x0000, - PM_SOURCE_PHOT_GROWTH = 0x0001, - PM_SOURCE_PHOT_APCORR = 0x0002, - PM_SOURCE_PHOT_WEIGHT = 0x0004, - PM_SOURCE_PHOT_INTERP = 0x0008, + PM_SOURCE_PHOT_NONE = 0x0000, + PM_SOURCE_PHOT_GROWTH = 0x0001, + PM_SOURCE_PHOT_APCORR = 0x0002, + PM_SOURCE_PHOT_WEIGHT = 0x0004, + PM_SOURCE_PHOT_INTERP = 0x0008, + PM_SOURCE_PHOT_DIFFSTATS = 0x0010, } pmSourcePhotometryMode; bool pmSourcePhotometryModel( float *fitMag, ///< integrated fit magnitude + float *fitFlux, ///< integrated fit magnitude pmModel *model ///< model used for photometry ); @@ -54,4 +56,5 @@ bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor); +bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal); double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor); Index: /branches/eam_branches/20100225/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/20100225/psModules/src/psmodules.h (revision 27493) +++ /branches/eam_branches/20100225/psModules/src/psmodules.h (revision 27494) @@ -115,4 +115,6 @@ #include #include +#include +#include #include #include