Index: /trunk/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- /trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 31450) +++ /trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 31451) @@ -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: /trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 31450) +++ /trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 31451) @@ -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: /trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 31450) +++ /trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 31451) @@ -1197,8 +1197,8 @@ // fprintf (stderr, "%f,%f : %f %f : %f %f\n", source->peak->xf, source->peak->yf, - // source->psfMag, source->apMag, source->psfMag - source->apMag, source->errMag); + // source->psfMag, source->apMag, source->psfMag - source->apMag, source->psfMagErr); // XXX this is somewhat arbitrary... - if (source->errMag > 0.05) continue; + if (source->psfMagErr > 0.05) continue; if (fabs(source->psfMag - source->apMag) > 0.5) continue; Index: /trunk/psModules/src/objects/Makefile.am =================================================================== --- /trunk/psModules/src/objects/Makefile.am (revision 31450) +++ /trunk/psModules/src/objects/Makefile.am (revision 31451) @@ -114,4 +114,5 @@ pmTrend2D.h \ pmGrowthCurve.h \ + pmGrowthCurveGenerate.h \ pmSourceMatch.h \ pmDetEff.h \ Index: /trunk/psModules/src/objects/models/pmModel_DEV.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_DEV.c (revision 31450) +++ /trunk/psModules/src/objects/models/pmModel_DEV.c (revision 31451) @@ -268,5 +268,5 @@ float f0 = 1.0; float f1, f2; - for (z = DZ; z < 50; z += DZ) { + for (z = DZ; z < 150; z += DZ) { f1 = exp(-pow(z,ALPHA)); z += DZ; Index: /trunk/psModules/src/objects/models/pmModel_EXP.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_EXP.c (revision 31450) +++ /trunk/psModules/src/objects/models/pmModel_EXP.c (revision 31451) @@ -255,5 +255,5 @@ float f0 = 1.0; float f1, f2; - for (z = DZ; z < 50; z += DZ) { + for (z = DZ; z < 150; z += DZ) { f1 = exp(-sqrt(z)); z += DZ; Index: /trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 31450) +++ /trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 31451) @@ -239,5 +239,5 @@ float f0 = 1.0; float f1, f2; - for (z = DZ; z < 50; z += DZ) { + for (z = DZ; z < 150; z += DZ) { f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0); z += DZ; Index: /trunk/psModules/src/objects/models/pmModel_PS1_V1.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 31450) +++ /trunk/psModules/src/objects/models/pmModel_PS1_V1.c (revision 31451) @@ -261,5 +261,5 @@ float f0 = 1.0; float f1, f2; - for (z = DZ; z < 50; z += DZ) { + for (z = DZ; z < 150; z += DZ) { f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); z += DZ; Index: /trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 31450) +++ /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 31451) @@ -262,5 +262,5 @@ float f0 = 1.0; float f1, f2; - for (z = DZ; z < 50; z += DZ) { + for (z = DZ; z < 150; z += DZ) { f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); z += DZ; Index: /trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 31450) +++ /trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 31451) @@ -251,5 +251,5 @@ float f0 = 1.0; float f1, f2; - for (z = DZ; z < 50; z += DZ) { + for (z = DZ; z < 150; z += DZ) { f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7])); z += DZ; Index: /trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 31450) +++ /trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 31451) @@ -316,5 +316,5 @@ float f0 = 1.0; float f1, f2; - for (z = DZ; z < 50; z += DZ) { + for (z = DZ; z < 150; z += DZ) { // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25)); f1 = exp(-pow(z,PAR[PM_PAR_7])); Index: /trunk/psModules/src/objects/pmGrowthCurve.c =================================================================== --- /trunk/psModules/src/objects/pmGrowthCurve.c (revision 31450) +++ /trunk/psModules/src/objects/pmGrowthCurve.c (revision 31451) @@ -73,5 +73,5 @@ // Fractional pixel radii are not well defined; use integer pixel radii. Use 1 pixel steps // until the scaling factor steps in intervals larger than 1 pixel - float Rlin = 1.0 / (fR - 1.0); + // float Rlin = 1.0 / (fR - 1.0); growth->radius = psVectorAllocEmpty (NPTS, PS_DATA_F32); @@ -79,12 +79,13 @@ // there will be NPTS radii + a few extras float radius = minRadius; - while (radius < Rlin) { + while (radius < refRadius) { // fprintf (stderr, "r: %f\n", radius); - psVectorAppend (growth->radius, radius); + psVectorAppend (growth->radius, radius - 0.001); radius += 1.0; } + growth->refBin = growth->radius->n - 1; 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); @@ -96,6 +97,9 @@ growth->refRadius = refRadius; growth->maxRadius = maxRadius; - growth->apLoss = 0.0; - growth->fitMag = 0.0; + growth->fitMag = NAN; + growth->refMag = NAN; + growth->apRef = NAN; + growth->apLoss = NAN; + return growth; } Index: /trunk/psModules/src/objects/pmGrowthCurve.h =================================================================== --- /trunk/psModules/src/objects/pmGrowthCurve.h (revision 31450) +++ /trunk/psModules/src/objects/pmGrowthCurve.h (revision 31451) @@ -22,6 +22,8 @@ psF32 maxRadius; psF32 fitMag; + psF32 refMag; psF32 apRef; // apMag[refRadius] psF32 apLoss; // fitMag - apRef + int refBin; } pmGrowthCurve; Index: /trunk/psModules/src/objects/pmGrowthCurveGenerate.c =================================================================== --- /trunk/psModules/src/objects/pmGrowthCurveGenerate.c (revision 31450) +++ /trunk/psModules/src/objects/pmGrowthCurveGenerate.c (revision 31451) @@ -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" /*****************************************************************************/ @@ -197,5 +196,5 @@ // mask the given aperture and measure the apMag psImageKeepCircle (mask, xc, yc, radius, "OR", markVal); - if (!pmSourcePhotometryAper (&apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) { + if (!pmSourcePhotometryAper (NULL, &apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) { psFree (growth); psFree (view); @@ -226,2 +225,136 @@ return growth; } + +# define DEBUG 0 +# if (DEBUG) +static FILE *fgr = NULL; +# endif + +// 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); + +# if (DEBUG) + fgr = fopen ("growth.mags.dat", "w"); +# endif + + 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 (DEBUG) + fclose (fgr); +# 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->refMag); + } + 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; + } + psAssert(growth->refBin >= 0, "invalid growth reference bin"); + psAssert(growth->refBin < growth->apMag->n, "invalid growth reference bin"); + growth->refMag = growth->apMag->data.F32[growth->refBin]; + + // Loop over the range of radii +# if (DEBUG) + for (int i = 0; i < growth->radius->n; i++) { + fprintf (fgr, "%f %f %f %f %f %f\n", xc, yc, growth->radius->data.F32[i], growth->apMag->data.F32[i], growth->fitMag, growth->refMag); + } +# endif + + return growth; +} Index: /trunk/psModules/src/objects/pmGrowthCurveGenerate.h =================================================================== --- /trunk/psModules/src/objects/pmGrowthCurveGenerate.h (revision 31451) +++ /trunk/psModules/src/objects/pmGrowthCurveGenerate.h (revision 31451) @@ -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: /trunk/psModules/src/objects/pmPCMdata.c =================================================================== --- /trunk/psModules/src/objects/pmPCMdata.c (revision 31450) +++ /trunk/psModules/src/objects/pmPCMdata.c (revision 31451) @@ -263,8 +263,11 @@ bool pmPCMupdate(pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model) { - bool newWindow = (source->pixels->numRows != pcm->modelFlux->numRows) || (source->pixels->numCols != pcm->modelFlux->numCols); + bool sameWindow = (source->pixels->numRows == pcm->modelFlux->numRows); + sameWindow &= (source->pixels->numCols == pcm->modelFlux->numCols); + sameWindow &= (source->pixels->col0 == pcm->modelFlux->col0); + sameWindow &= (source->pixels->row0 == pcm->modelFlux->row0); // re-count the number of unmasked pixels: - if (newWindow) { + if (!sameWindow) { for (psS32 i = 0; i < source->pixels->numRows; i++) { for (psS32 j = 0; j < source->pixels->numCols; j++) { @@ -346,5 +349,5 @@ // has the source pixel window changed? - if (newWindow) { + if (!sameWindow) { // adjust all supporting images: Index: /trunk/psModules/src/objects/pmPSF.c =================================================================== --- /trunk/psModules/src/objects/pmPSF.c (revision 31450) +++ /trunk/psModules/src/objects/pmPSF.c (revision 31451) @@ -73,7 +73,5 @@ options->type = 0; - options->stats = NULL; - options->fitOptions = NULL; // XXX this has to be set before calling pmPSF fit functions options->psfTrendMode = PM_TREND_NONE; @@ -90,5 +88,8 @@ options->chiFluxTrend = true; - + options->fitOptions = NULL; // XXX this has to be set before calling pmPSF fit functions + + options->fitRadius = NAN; + options->apRadius = NAN; return options; } @@ -141,4 +142,5 @@ psf->type = options->type; + psf->chisq = 0.0; psf->ApResid = 0.0; Index: /trunk/psModules/src/objects/pmPSF.h =================================================================== --- /trunk/psModules/src/objects/pmPSF.h (revision 31450) +++ /trunk/psModules/src/objects/pmPSF.h (revision 31451) @@ -31,17 +31,28 @@ struct pmPSF { pmModelType type; ///< PSF Model in use - psArray *params; ///< Model parameters (psPolynomial2D) - psStats *psfTrendStats; ///< psf parameter trend clipping stats - pmTrend2DMode psfTrendMode; - psPolynomial1D *ChiTrend; ///< Chisq vs flux fit (correction for systematic errors) - pmTrend2D *ApTrend; ///< ApResid vs (x,y) - pmTrend2D *FluxScale; ///< Flux for PSF at (x,y) for normalization = 1.0 + + float chisq; ///< PSF goodness statistic (unused??) float ApResid; ///< apMag - psfMag (for PSF stars) float dApResid; ///< scatter of ApResid float skyBias; ///< implied residual sky offset from ApResid fit float skySat; ///< roll-over of ApResid fit - float chisq; ///< PSF goodness statistic (unused??) int nPSFstars; ///< number of stars used to measure PSF int nApResid; ///< number of stars used to measure ApResid + + bool poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting + bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting + bool poissonErrorsParams; ///< use poission errors for model parameter fitting + + pmTrend2D *ApTrend; ///< ApResid vs (x,y) + pmTrend2D *FluxScale; ///< Flux for PSF at (x,y) for normalization = 1.0 + psPolynomial1D *ChiTrend; ///< Chisq vs flux fit (correction for systematic errors) + + pmGrowthCurve *growth; ///< apMag vs Radius + pmResiduals *residuals; ///< normalized residual image (no spatial variation) + + psArray *params; ///< Model parameters (psPolynomial2D) + psStats *psfTrendStats; ///< psf parameter trend clipping stats + + pmTrend2DMode psfTrendMode; int trendNx; int trendNy; @@ -50,9 +61,4 @@ int fieldXo; int fieldYo; - bool poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting - bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting - bool poissonErrorsParams; ///< use poission errors for model parameter fitting - pmGrowthCurve *growth; ///< apMag vs Radius - pmResiduals *residuals; ///< normalized residual image (no spatial variation) }; @@ -60,4 +66,5 @@ pmModelType type; psStats *stats; // psfTrend clipping stats + pmTrend2DMode psfTrendMode; int psfTrendNx; @@ -67,11 +74,14 @@ int psfFieldXo; int psfFieldYo; + bool poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting bool poissonErrorsParams; ///< use poission errors for model parameter fitting + + bool chiFluxTrend; // Fit a trend in Chi2 as a function of flux? + pmSourceFitOptions *fitOptions; + float fitRadius; float apRadius; - bool chiFluxTrend; // Fit a trend in Chi2 as a function of flux? - pmSourceFitOptions *fitOptions; } pmPSFOptions; @@ -94,5 +104,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: /trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- /trunk/psModules/src/objects/pmPSF_IO.c (revision 31450) +++ /trunk/psModules/src/objects/pmPSF_IO.c (revision 31451) @@ -63,4 +63,9 @@ bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header); +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) @@ -197,4 +202,80 @@ psError(psErrorCodeLast(), false, "Failed to write PSF for chip"); return false; + } + return true; +} + +// XXX we save the model term identifiers (item) as S32, but they probably should be more flexible +bool pmTrend2DtoTable (psArray *table, pmTrend2D *trend, char *label, int item) { + + if (trend == NULL) return true; + + if (trend->mode == PM_TREND_MAP) { + // write the image components into a table: this is needed because they may each be a different size + psImageMap *map = trend->map; + for (int ix = 0; ix < map->map->numCols; ix++) { + for (int iy = 0; iy < map->map->numRows; iy++) { + psMetadata *row = psMetadataAlloc (); + psMetadataAddS32 (row, PS_LIST_TAIL, label, 0, "", item); + psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); + psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); + psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", map->map->data.F32[iy][ix]); + psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", map->error->data.F32[iy][ix]); + psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", 0); // no cells are masked + + psArrayAdd (table, 100, row); + psFree (row); + } + } + } else { + // write the polynomial components into a table + psPolynomial2D *poly = trend->poly; + for (int ix = 0; ix <= poly->nX; ix++) { + for (int iy = 0; iy <= poly->nY; iy++) { + psMetadata *row = psMetadataAlloc (); + psMetadataAddS32 (row, PS_LIST_TAIL, label, 0, "", item); + psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); + psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); + psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); + psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); + psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->coeffMask[ix][iy]); + + psArrayAdd (table, 100, row); + psFree (row); + } + } + } + return true; +} + +// extra trend2D elements from a row +bool pmTrend2DfromTableRow (pmTrend2D *trend, psMetadata *row) { + + bool status = false; + + int xPow = psMetadataLookupS32 (&status, row, "X_POWER"); + int yPow = psMetadataLookupS32 (&status, row, "Y_POWER"); + + if (trend->mode == PM_TREND_MAP) { + psImageMap *map = trend->map; + assert (map); + assert (map->map); + assert (map->error); + assert (xPow >= 0); + assert (yPow >= 0); + assert (xPow < map->map->numCols); + assert (yPow < map->map->numRows); + map->map->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "VALUE"); + map->error->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "ERROR"); + } else { + psPolynomial2D *poly = trend->poly; + assert (poly); + assert (xPow >= 0); + assert (yPow >= 0); + assert (xPow <= poly->nX); + assert (yPow <= poly->nY); + poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); + poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); + poly->coeffMask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); } return true; @@ -403,9 +484,9 @@ if (trend->mode == PM_TREND_MAP) { - nX = trend->map->map->numCols; - nY = trend->map->map->numRows; + nX = trend->map->map->numCols; + nY = trend->map->map->numRows; } else { - nX = trend->poly->nX; - nY = trend->poly->nY; + nX = trend->poly->nX; + nY = trend->poly->nY; } snprintf (name, 9, "PAR%02d_NX", i); @@ -428,45 +509,18 @@ psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); + float PSF_APERTURE = psMetadataLookupF32(&status, roAnalysis, "PSF_APERTURE"); + if (status) { + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE); + } + float PSF_FIT_RADIUS = psMetadataLookupF32(&status, roAnalysis, "PSF_FIT_RADIUS"); + if (status) { + psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS); + } + // build a FITS table of the PSF parameters psArray *psfTable = psArrayAllocEmpty (100); for (int i = 0; i < nPar; i++) { pmTrend2D *trend = psf->params->data[i]; - if (trend == NULL) continue; // skip unset parameters (eg, XPOS) - - if (trend->mode == PM_TREND_MAP) { - // write the image components into a table: this is needed because they may each be a different size - psImageMap *map = trend->map; - for (int ix = 0; ix < map->map->numCols; ix++) { - for (int iy = 0; iy < map->map->numRows; iy++) { - psMetadata *row = psMetadataAlloc (); - psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); - psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); - psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); - psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", map->map->data.F32[iy][ix]); - psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", map->error->data.F32[iy][ix]); - psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", 0); // no cells are masked - - psArrayAdd (psfTable, 100, row); - psFree (row); - } - } - } else { - // write the polynomial components into a table - psPolynomial2D *poly = trend->poly; - for (int ix = 0; ix <= poly->nX; ix++) { - for (int iy = 0; iy <= poly->nY; iy++) { - psMetadata *row = psMetadataAlloc (); - psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); - psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); - psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); - psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); - psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); - psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->coeffMask[ix][iy]); - - psArrayAdd (psfTable, 100, row); - psFree (row); - } - } - } + pmTrend2DtoTable (psfTable, trend, "MODEL_TERM", i); } @@ -547,20 +601,17 @@ } + 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; + } + // write a representation of the psf model { psMetadata *header = psMetadataAlloc (); - - if (0) { - // set some header keywords to make it clear there are no residuals? - if (!psFitsWriteBlank (file->fits, header, residName)) { - psError(psErrorCodeLast(), false, "Unable to write blank PSF residuals."); - psFree(residName); - psFree(header); - return false; - } - psFree (residName); - psFree (header); - return true; - } int DX = 65; @@ -651,4 +702,6 @@ } + + // if this file needs to have a PHU written out, write one bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config) @@ -934,4 +987,17 @@ psf->skyBias = psMetadataLookupF32 (&status, header, "SKY_BIAS"); + if (roAnalysis) { + float PSF_APERTURE = psMetadataLookupF32(&status, header, "PSF_APERTURE"); + if (status) { + psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE); + } + float PSF_FIT_RADIUS = psMetadataLookupF32(&status, header, "PSF_FIT_RADIUS"); + if (status) { + psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS); + } + } else { + psWarning ("unable to read PSF_APERTURE or PSF_FIT_RADIUS"); + } + // read the raw table data psArray *table = psFitsReadTable (file->fits); @@ -945,7 +1011,6 @@ for (int i = 0; i < table->n; i++) { psMetadata *row = table->data[i]; + int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM"); - int xPow = psMetadataLookupS32 (&status, row, "X_POWER"); - int yPow = psMetadataLookupS32 (&status, row, "Y_POWER"); pmTrend2D *trend = psf->params->data[iPar]; @@ -955,26 +1020,5 @@ } - if (trend->mode == PM_TREND_MAP) { - psImageMap *map = trend->map; - assert (map); - assert (map->map); - assert (map->error); - assert (xPow >= 0); - assert (yPow >= 0); - assert (xPow < map->map->numCols); - assert (yPow < map->map->numRows); - map->map->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "VALUE"); - map->error->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "ERROR"); - } else { - psPolynomial2D *poly = trend->poly; - assert (poly); - assert (xPow >= 0); - assert (yPow >= 0); - assert (xPow <= poly->nX); - assert (yPow <= poly->nY); - poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); - poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); - poly->coeffMask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); - } + pmTrend2DfromTableRow(trend, row); } psFree (header); @@ -1063,4 +1107,14 @@ } + if (!pmPSFmodelRead_ApTrend (psf, file)) { + psError(psErrorCodeLast(), false, "Unable to read PSF ApTrend data."); + return false; + } + + 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); @@ -1070,4 +1124,217 @@ psFree (header); + return true; +} + +// write aperture trend to a FITS table +bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf) { + + pmTrend2D *trend = psf->ApTrend; + if (trend == NULL) { + psWarning ("no PSF ApTrend to write out, skipping"); + return true; + } + + // we need to write a header for the table, + psMetadata *header = psMetadataAlloc(); + + int nX = 0, nY = 0; + if (trend->mode == PM_TREND_MAP) { + nX = trend->map->map->numCols; + nY = trend->map->map->numRows; + } else { + nX = trend->poly->nX; + nY = trend->poly->nY; + } + psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NX", 0, "", nX); + psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NY", 0, "", nY); + char *modeName = pmTrend2DModeToString (trend->mode); + psMetadataAddStr (header, PS_LIST_TAIL, "TREND_MD", 0, "", modeName); + psFree (modeName); + + // build a FITS table of the ApTrend (only 1) + psArray *table = psArrayAllocEmpty (100); + pmTrend2DtoTable (table, trend, "APTREND", 0); + + // 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 ApTrend data %s\n", "AP_TREND"); + if (!psFitsWriteTable(file->fits, header, table, "AP_TREND")) { + psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "AP_TREND"); + psFree(table); + psFree(header); + return false; + } + + psFree (table); + psFree (header); + return true; +} + +// read aperture trend to a FITS table +bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file) { + + bool status; + + // move fits pointer to AP_TREND section + // advance to the table data extension + if (!psFitsMoveExtNameClean (file->fits, "AP_TREND")) { + psWarning ("no Aperture Trend data in PSF file, skipping"); + return true; + } + + psMetadata *header = psFitsReadHeader (NULL, file->fits); + if (!header) { + psError(psErrorCodeLast(), false, "Unable to read AP_TREND header."); + return false; + } + + // read the raw table data + psArray *table = psFitsReadTable (file->fits); + if (!table) { + psError(psErrorCodeLast(), false, "Unable to read AP_TREND table."); + psFree(header); + return false; + } + + // XXX allow user to set this optionally? + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = psf->fieldNx; + binning->nYfine = psf->fieldNy; + binning->nXruff = psMetadataLookupS32 (&status, header, "TREND_NX"); + binning->nYruff = psMetadataLookupS32 (&status, header, "TREND_NY"); + psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); + char *modeName = psMetadataLookupStr (&status, header, "TREND_MD"); + if (!status) { + psError(PM_ERR_PROG, true, "inconsistent PSF header: NX & NY defined for AP TREND, but not MD"); + psFree (header); + psFree (stats); + psFree (table); + return false; + } + pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName); + if (psfTrendMode == PM_TREND_NONE) { + psfTrendMode = PM_TREND_POLY_ORD; + } + + // measure Trend2D for the current spatial scale + pmTrend2D *apTrend = pmTrend2DNoImageAlloc (PM_TREND_MAP, binning, stats); + + // fill in the matching psf->params entries + for (int i = 0; i < table->n; i++) { + psMetadata *row = table->data[i]; + pmTrend2DfromTableRow(apTrend, row); + } + psf->ApTrend = apTrend; + + psFree (binning); + psFree (header); + psFree (stats); + psFree (table); + return true; +} + +// 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; } Index: /trunk/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 31450) +++ /trunk/psModules/src/objects/pmPSFtryFitPSF.c (revision 31451) @@ -122,5 +122,5 @@ psfTry->fitMag->data.F32[i] = source->psfMag; psfTry->metric->data.F32[i] = source->apMag - source->psfMag; - psfTry->metricErr->data.F32[i] = source->errMag; + psfTry->metricErr->data.F32[i] = source->psfMagErr; // XXX this did not work: modifies shape of psf too much @@ -131,5 +131,5 @@ source->peak->xf, source->peak->yf, source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], - source->psfMag, source->apMag, source->errMag, + source->psfMag, source->apMag, source->psfMagErr, source->modelPSF->params->data.F32[PM_PAR_I0], source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], Index: /trunk/psModules/src/objects/pmPeaks.c =================================================================== --- /trunk/psModules/src/objects/pmPeaks.c (revision 31450) +++ /trunk/psModules/src/objects/pmPeaks.c (revision 31451) @@ -150,4 +150,8 @@ tmp->x = x; tmp->y = y; + tmp->xf = x; + tmp->yf = y; + tmp->dx = NAN; + tmp->dy = NAN; tmp->detValue = value; tmp->rawFlux = value; // set this by default: it is up to the user to supply a better value @@ -155,7 +159,4 @@ tmp->smoothFlux = value; // set this by default: it is up to the user to supply a better value tmp->smoothFluxStdev = NAN; - // tmp->SN = 0; - tmp->xf = x; - tmp->yf = y; tmp->assigned = false; tmp->type = type; @@ -166,4 +167,24 @@ psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__); return(tmp); +} + +// copy to an already allocated peak +bool pmPeakCopy(pmPeak *out, pmPeak *in) +{ + out->x = in->x; + out->y = in->y; + out->xf = in->xf; + out->yf = in->yf; + out->dx = in->dx; + out->dy = in->dy; + out->detValue = in->detValue; + out->rawFlux = in->rawFlux; + out->rawFluxStdev = in->rawFluxStdev; + out->smoothFlux = in->smoothFlux; + out->smoothFluxStdev = in->smoothFluxStdev; + out->assigned = in->assigned; + out->type = in->type; + + return true; } Index: /trunk/psModules/src/objects/pmPeaks.h =================================================================== --- /trunk/psModules/src/objects/pmPeaks.h (revision 31450) +++ /trunk/psModules/src/objects/pmPeaks.h (revision 31451) @@ -69,5 +69,4 @@ float smoothFlux; ///< peak flux in smoothed signal image float smoothFluxStdev; ///< peak stdev in smoothed signal image - // float SNestimated; ///< S/N estimated from the detection image bool assigned; ///< is peak assigned to a source? pmPeakType type; ///< Description of peak. @@ -89,4 +88,6 @@ bool psMemCheckPeak(psPtr ptr); + +bool pmPeakCopy(pmPeak *out, pmPeak *in); /** pmPeaksInVector() Index: /trunk/psModules/src/objects/pmSource.c =================================================================== --- /trunk/psModules/src/objects/pmSource.c (revision 31450) +++ /trunk/psModules/src/objects/pmSource.c (revision 31451) @@ -116,5 +116,4 @@ source->psfImage = NULL; source->moments = NULL; - source->blends = NULL; source->modelPSF = NULL; source->modelEXT = NULL; @@ -124,18 +123,11 @@ source->mode2 = PM_SOURCE_MODE_DEFAULT; source->tmpFlags = 0; - source->extpars = NULL; - source->diffStats = NULL; - source->radialAper = NULL; - source->parent = NULL; - - source->region = psRegionSet(NAN, NAN, NAN, NAN); - psMemSetDeallocator(source, (psFreeFunc) sourceFree); // default values are NAN source->psfMag = NAN; + source->psfMagErr = NAN; source->psfFlux = NAN; source->psfFluxErr = NAN; source->extMag = NAN; - source->errMag = NAN; source->apMag = NAN; source->apMagRaw = NAN; @@ -143,6 +135,5 @@ source->apFlux = NAN; source->apFluxErr = NAN; - source->sky = NAN; - source->skyErr = NAN; + source->pixWeightNotBad = NAN; source->pixWeightNotPoor = NAN; @@ -151,4 +142,16 @@ source->crNsigma = NAN; source->extNsigma = NAN; + source->sky = NAN; + source->skyErr = NAN; + + source->region = psRegionSet(NAN, NAN, NAN, NAN); + source->blends = NULL; + source->extpars = NULL; + source->diffStats = NULL; + source->radialAper = NULL; + source->parent = NULL; + source->imageID = -1; + + psMemSetDeallocator(source, (psFreeFunc) sourceFree); psTrace("psModules.objects", 10, "---- end ----\n"); @@ -172,11 +175,5 @@ if (in->peak != NULL) { source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->detValue, in->peak->type); - source->peak->xf = in->peak->xf; - source->peak->yf = in->peak->yf; - source->peak->rawFlux = in->peak->rawFlux; - source->peak->rawFluxStdev = in->peak->rawFluxStdev; - source->peak->smoothFlux = in->peak->smoothFlux; - source->peak->smoothFluxStdev = in->peak->smoothFluxStdev; - // source->peak->SN = in->peak->SN; + pmPeakCopy(source->peak, in->peak); } @@ -195,9 +192,31 @@ // the maskObj is a unique mask array; create a new mask image - source->maskObj = in->maskObj ? psImageCopy (NULL, in->maskObj, PS_TYPE_IMAGE_MASK) : NULL; - - source->type = in->type; - source->mode = in->mode; - source->imageID = in->imageID; + source->maskObj = in->maskObj ? psImageCopy (NULL, in->maskObj, PS_TYPE_IMAGE_MASK) : NULL; + + // NOTE : because of the const id element, we cannot just assign *source = *in + + source->type = in->type; + source->mode = in->mode; + source->mode2 = in->mode2; + source->tmpFlags = in->tmpFlags; + source->psfMag = in->psfMag; + source->psfMagErr = in->psfMagErr; + source->psfFlux = in->psfFlux; + source->psfFluxErr = in->psfFluxErr; + source->extMag = in->extMag; + source->apMag = in->apMag; + source->apMagRaw = in->apMagRaw; + source->apRadius = in->apRadius; + source->apFlux = in->apFlux; + source->apFluxErr = in->apFluxErr; + source->pixWeightNotBad = in->pixWeightNotBad; + source->pixWeightNotPoor = in->pixWeightNotPoor; + source->psfChisq = in->psfChisq; + source->crNsigma = in->crNsigma; + source->extNsigma = in->extNsigma; + source->sky = in->sky; + source->skyErr = in->skyErr; + + source->region = in->region; return(source); @@ -266,4 +285,16 @@ extend |= (mySource->maskObj == NULL); extend |= (mySource->maskView == NULL); + + // if ((fabs(x - 2020) < 5) && (fabs(y - 366) < 5)) { + // if (extend) { + // fprintf (stderr, "extend T, %f, %f : %f, %f vs %f, %f : %f, %f\n", + // newRegion.x0, newRegion.y0, newRegion.x1, newRegion.y1, + // mySource->region.x0, mySource->region.y0, mySource->region.x1, mySource->region.y1); + // } else { + // fprintf (stderr, "extend F, %f, %f : %f, %f vs %f, %f : %f, %f\n", + // newRegion.x0, newRegion.y0, newRegion.x1, newRegion.y1, + // mySource->region.x0, mySource->region.y0, mySource->region.x1, mySource->region.y1); + // } + // } if (extend) { @@ -486,5 +517,5 @@ // create vectors with Sx, Sy values in window // clip sources based on S/N - for (psS32 i = 0 ; i < sources->n ; i++) + for (psS32 i = 0; i < sources->n; i++) { pmSource *tmpSrc = (pmSource *) sources->data[i]; @@ -1175,5 +1206,5 @@ if (!source->moments) return false; // can't if there are no moments if (!source->moments->nPixels) return false; // can't if the moments were not measured - if (source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE) return false; // can't if the moments failed... + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false; // can't if the moments failed... if (source->mode & PM_SOURCE_MODE_SATSTAR) return true; // moments are best for SATSTARs Index: /trunk/psModules/src/objects/pmSource.h =================================================================== --- /trunk/psModules/src/objects/pmSource.h (revision 31450) +++ /trunk/psModules/src/objects/pmSource.h (revision 31451) @@ -80,10 +80,10 @@ pmSourceMode2 mode2; ///< analysis flags set for object. pmSourceTmpF tmpFlags; ///< internal-only flags - 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 @@ -98,6 +98,9 @@ float crNsigma; ///< Nsigma deviation from PSF to CR float extNsigma; ///< Nsigma deviation from PSF to EXT - float sky, skyErr; ///< The sky and its error at the center of the object + float sky; ///< The sky at the center of the object + float skyErr; ///< The sky error at the center of the object + psRegion region; ///< area on image covered by selected pixels + psArray *blends; ///< collection of sources thought to be confused with object pmSourceExtendedPars *extpars; ///< extended source parameters pmSourceDiffStats *diffStats; ///< extra parameters for difference detections @@ -254,4 +257,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 +277,6 @@ float xGuess, float yGuess); +float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM); + pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source); Index: /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 31451) @@ -125,5 +125,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 magnitude", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); @@ -256,10 +256,10 @@ 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); Index: /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 31451) @@ -124,5 +124,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 magnitude", source->psfFlux); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); @@ -274,10 +274,10 @@ 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); Index: /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 31451) @@ -125,5 +125,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); @@ -273,10 +273,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); @@ -720,6 +720,6 @@ sprintf (keyword1, "RMIN_%02d", i); sprintf (keyword2, "RMAX_%02d", i); - psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); + psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); + psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]); } Index: /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 31451) @@ -118,5 +118,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, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); @@ -241,10 +241,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); Index: /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 31451) @@ -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, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); @@ -247,10 +247,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); Index: /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c (revision 31451) @@ -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: /trunk/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 31451) @@ -137,5 +137,5 @@ axes = pmPSF_ModelToAxes (PAR, 20.0); - float errMag = isfinite(source->errMag) ? source->errMag : 999; + float psfMagErr = isfinite(source->psfMagErr) ? source->psfMagErr : 999; psLineInit (line); @@ -143,5 +143,5 @@ psLineAdd (line, "%6.1f ", PAR[PM_PAR_YPOS]); psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->psfMag + ZERO_POINT)); - psLineAdd (line, "%03d ", PS_MIN (999, (int)errMag)); + psLineAdd (line, "%03d ", PS_MIN (999, (int)psfMagErr)); psLineAdd (line, "%2d ", type); psLineAdd (line, "%3.1f ", lsky); @@ -280,5 +280,5 @@ source->psfMag = atof (array->data[2]); source->extMag = atof (array->data[6]); - source->errMag = atof (array->data[3]) / 1000.0; + source->psfMagErr = atof (array->data[3]) / 1000.0; source->apMag = atof (array->data[7]); axes.major = atof (array->data[8]); Index: /trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 31451) @@ -71,5 +71,5 @@ psF32 xPos, yPos; psF32 xErr, yErr; - psF32 errMag, chisq; + psF32 psfMagErr, chisq; pmChip *chip = readout->parent->parent; @@ -121,6 +121,6 @@ chisq = model->chisq; - // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag. - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + // need to determine the PSF photometry error: source->psfMagErr is the error on the 'best' model mag. + psfMagErr = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; } else { xPos = source->peak->xf; @@ -132,5 +132,5 @@ axes.theta = NAN; chisq = NAN; - errMag = NAN; + psfMagErr = NAN; } @@ -167,7 +167,7 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "Calibrated Magnitude from PSF Fit", calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "Calibrated Magnitude Error", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "Calibrated Magnitude Error", psfMagErr); 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", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", psfMagErr); psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); @@ -284,7 +284,7 @@ // 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"); 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); Index: /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c (revision 31451) @@ -113,5 +113,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", PS_MIN (99.0, source->psfMag)); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", PS_MIN (99.0, source->errMag)); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", PS_MIN (99.0, source->psfMagErr)); psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", PS_MIN (99.0, peakMag)); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); @@ -212,5 +212,5 @@ // 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"); pmPSF_AxesToModel (PAR, axes); Index: /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 31451) @@ -70,5 +70,5 @@ psF32 xPos, yPos; psF32 xErr, yErr; - psF32 errMag, chisq; + psF32 psfMagErr, chisq; // let's write these out in S/N order @@ -103,6 +103,6 @@ chisq = model->chisq; - // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag. - errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + // need to determine the PSF photometry error: source->psfMagErr is the error on the 'best' model mag. + psfMagErr = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; } else { xPos = source->peak->xf; @@ -114,5 +114,5 @@ axes.theta = NAN; chisq = NAN; - errMag = NAN; + psfMagErr = NAN; } @@ -128,5 +128,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); 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", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", psfMagErr); psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); @@ -253,7 +253,7 @@ // 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"); 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); Index: /trunk/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 31451) @@ -110,5 +110,5 @@ fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], source->sky, - source->psfMag, source->errMag, dPos); + source->psfMag, source->psfMagErr, dPos); for (j = 4; j < model->params->n; j++) { @@ -174,5 +174,5 @@ fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], source->sky, - source->extMag, source->errMag, dPos); + source->extMag, source->psfMagErr, dPos); for (j = 4; j < model->params->n; j++) { Index: /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c (revision 31451) @@ -106,5 +106,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", yPos); psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", PS_MIN (99.0, source->psfMag + ZERO_POINT)); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", PS_MIN (999, 1000*source->errMag)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", PS_MIN (999, 1000*source->psfMagErr)); psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", PS_MIN (99.0, source->extMag + ZERO_POINT)); psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", PS_MIN (99.0, source->apMag + ZERO_POINT)); @@ -194,5 +194,5 @@ source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; source->extMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT; - source->errMag = psMetadataLookupF32 (&status, row, "MAG_ERR") * 0.001; + source->psfMagErr = psMetadataLookupF32 (&status, row, "MAG_ERR") * 0.001; source->apMag = psMetadataLookupF32 (&status, row, "MAG_AP") - ZERO_POINT; Index: /trunk/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_SX.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceIO_SX.c (revision 31451) @@ -88,5 +88,5 @@ psLineAdd (line, "%11.3f", PAR[PM_PAR_YPOS]); psLineAdd (line, "%9.4f", source->psfMag); - psLineAdd (line, "%9.4f", source->errMag); + psLineAdd (line, "%9.4f", source->psfMagErr); psLineAdd (line, "%13.4f", source->sky); psLineAdd (line, "%9.2f", axes.major); Index: /trunk/psModules/src/objects/pmSourceMasks.h =================================================================== --- /trunk/psModules/src/objects/pmSourceMasks.h (revision 31450) +++ /trunk/psModules/src/objects/pmSourceMasks.h (revision 31451) @@ -14,5 +14,5 @@ PM_SOURCE_MODE_PSFSTAR = 0x00000040, ///< Source used to define PSF model PM_SOURCE_MODE_SATSTAR = 0x00000080, ///< Source model peak is above saturation - PM_SOURCE_MODE_BLEND = 0x00000100, ///< Source is a blend with other sourcers + PM_SOURCE_MODE_BLEND = 0x00000100, ///< Source is a blend with other sources PM_SOURCE_MODE_EXTERNAL = 0x00000200, ///< Source based on supplied input position PM_SOURCE_MODE_BADPSF = 0x00000400, ///< Failed to get good estimate of object's PSF @@ -51,4 +51,6 @@ PM_SOURCE_MODE2_ON_BURNTOOL = 0x00000020, ///< > 25% of (PSF-weighted) pixels land on burntool PM_SOURCE_MODE2_ON_CONVPOOR = 0x00000040, ///< > 25% of (PSF-weighted) pixels land on convpoor + + PM_SOURCE_MODE2_PASS1_SRC = 0x00000080, ///< source detected in first pass analysis } pmSourceMode2; Index: /trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceMatch.c (revision 31451) @@ -88,5 +88,5 @@ pmSource *source = sources->data[i]; // Source of interest if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) || - !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) { + !isfinite(source->psfMagErr) || source->psfMag > SOURCE_FAINTEST) { continue; } @@ -102,5 +102,5 @@ (*y)->data.F32[num] = ySrc; (*mag)->data.F32[num] = source->psfMag; - (*magErr)->data.F32[num] = source->errMag; + (*magErr)->data.F32[num] = source->psfMagErr; (*indices)->data.S32[num] = i; num++; Index: /trunk/psModules/src/objects/pmSourceMoments.c =================================================================== --- /trunk/psModules/src/objects/pmSourceMoments.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceMoments.c (revision 31451) @@ -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,35 +109,48 @@ // Xn = SUM (x - xc)^n * (z - sky) - 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 + + // float dX = source->moments->Mx - source->peak->xf; + // float dY = source->moments->My - source->peak->yf; + // float dR = hypot(dX, dY); + // float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf; + // float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf; + float Xo = source->moments->Mx; + float Yo = source->moments->My; // 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 = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage - psF32 yCM = source->moments->My - 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]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -151,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 @@ -171,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; @@ -182,30 +195,36 @@ // 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 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; RH += rh; RS += rs; + + RFW += rfw; + RHW += rhw; XX += xx; @@ -244,8 +263,12 @@ source->moments->Myyyy = YYYY/Sum; - // 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; @@ -259,10 +282,12 @@ 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]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -276,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 @@ -290,5 +315,5 @@ if (PS_SQR(pDiff) < minSN2*wDiff) continue; - psF32 r = sqrt(r2); + float r = sqrt(r2); if (r < radKron) { Sum += pDiff; @@ -335,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) @@ -342,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 @@ -358,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"); @@ -370,10 +395,12 @@ 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]; + // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row]; for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { @@ -387,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 @@ -408,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; @@ -419,6 +446,6 @@ Sum += pDiff; - psF32 xWght = xDiff * pDiff; - psF32 yWght = yDiff * pDiff; + float xWght = xDiff * pDiff; + float yWght = yDiff * pDiff; X1 += xWght; @@ -477,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: /trunk/psModules/src/objects/pmSourceOutputs.c =================================================================== --- /trunk/psModules/src/objects/pmSourceOutputs.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourceOutputs.c (revision 31451) @@ -98,12 +98,13 @@ outputs->xPos = PAR[PM_PAR_XPOS]; outputs->yPos = PAR[PM_PAR_YPOS]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + if ((source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) && !(source->mode & PM_SOURCE_MODE_EXTMODEL)) { + // we only do non-linear PSF fits for non-extended objects outputs->xErr = dPAR[PM_PAR_XPOS]; outputs->yErr = dPAR[PM_PAR_YPOS]; } else { - outputs->xErr = fwhmMajor * source->errMag / 2.35; - outputs->yErr = fwhmMinor * source->errMag / 2.35; + outputs->xErr = fwhmMajor * source->psfMagErr / 2.35; + outputs->yErr = fwhmMinor * source->psfMagErr / 2.35; } - if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) { + if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY])) { axes = pmPSF_ModelToAxes (PAR, 20.0); outputs->psfMajor = axes.major; @@ -125,6 +126,6 @@ outputs->xPos = source->moments->Mx; outputs->yPos = source->moments->My; - outputs->xErr = fwhmMajor * source->errMag / 2.35; - outputs->yErr = fwhmMinor * source->errMag / 2.35; + outputs->xErr = fwhmMajor * source->psfMagErr / 2.35; + outputs->yErr = fwhmMinor * source->psfMagErr / 2.35; } else { outputs->xPos = source->peak->xf; Index: /trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /trunk/psModules/src/objects/pmSourcePhotometry.c (revision 31450) +++ /trunk/psModules/src/objects/pmSourcePhotometry.c (revision 31451) @@ -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,13 +97,10 @@ int status = false; - bool isPSF; float x, y; - float rflux; float SN; - pmModel *model; source->psfMag = NAN; source->extMag = NAN; - source->errMag = NAN; + source->psfMagErr = NAN; source->apMag = NAN; source->apMagRaw = NAN; @@ -114,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) { @@ -153,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) { @@ -173,15 +172,16 @@ } - // 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); source->psfMag += apTrend; - } - - // measure the contribution of included pixels + source->psfFlux *= pow(10.0, -0.4*apTrend); + source->psfFluxErr *= pow(10.0, -0.4*apTrend); + } + + // 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); } @@ -217,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(); @@ -232,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"); @@ -241,16 +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? - } - if (mode & PM_SOURCE_PHOT_APCORR) { - // XXX this should be removed -- we no longer fit for the 'sky bias' - // XXX is this happening??? - rflux = pow (10.0, 0.4*source->psfMag); - psAssert (psf->skyBias == 0.0, "sky bias not 0"); - psAssert (psf->skySat == 0.0, "sky sat not 0"); - source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux; + 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); } } @@ -309,11 +304,14 @@ bool status; - status = pmSourcePhotometryAper(&source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal); - + int nPix = 0; + status = pmSourcePhotometryAper(&nPix, &source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal); + if (status) { + source->mode |= PM_SOURCE_MODE_AP_MAGS; + } return status; } // return source aperture magnitude -bool pmSourcePhotometryAper (float *apMag, float *apFluxOut, float *apFluxErr, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal) +bool pmSourcePhotometryAper (int *nPixOut, float *apMag, float *apFluxOut, float *apFluxErr, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(apMag, false); @@ -328,4 +326,5 @@ float apFlux = 0; float apFluxVar = 0; + int nPix = 0; if (DO_SKY) { @@ -345,8 +344,11 @@ apFlux += imData[iy][ix] - sky; apFluxVar += varData[iy][ix]; - } - } + nPix ++; + } + } + if (apFluxOut) *apFluxOut = apFlux; if (apFluxErr) *apFluxErr = sqrt(fabs(apFluxVar)); + if (nPixOut) *nPixOut = nPix; if (apFlux <= 0) { @@ -412,4 +414,6 @@ psImageMaskType maskBad = maskVal; maskBad &= ~maskSuspect; + + psImageMaskType maskPoor = maskVal | maskSuspect; // measure modelSum and validSum. this function is applied to a sources' subimage. the @@ -446,5 +450,5 @@ // count pixels which are masked with an mask bit (bad or poor) - if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) { + if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskPoor)) { notPoorSum += value; } @@ -487,5 +491,5 @@ if (isfinite(source->pixWeightNotBad) && isfinite(source->pixWeightNotPoor)) { - psAssert (source->pixWeightNotBad <= source->pixWeightNotPoor, "error: all bad pixels should also be poor"); + psAssert (source->pixWeightNotPoor <= source->pixWeightNotBad, "error: all bad pixels should also be poor"); } Index: /trunk/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /trunk/psModules/src/objects/pmSourcePhotometry.h (revision 31450) +++ /trunk/psModules/src/objects/pmSourcePhotometry.h (revision 31451) @@ -45,4 +45,5 @@ bool pmSourcePhotometryAper( + int *nPixOut, float *apMag, float *apFluxOut, Index: /trunk/psModules/src/objects/pmSourceUtils.h =================================================================== --- /trunk/psModules/src/objects/pmSourceUtils.h (revision 31450) +++ /trunk/psModules/src/objects/pmSourceUtils.h (revision 31451) @@ -38,3 +38,3 @@ /// @} -# endif /* PM_MODEL_CLASS_H */ +# endif /* PM_SOURCE_UTILS_H */ Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 31450) +++ /trunk/psModules/src/psmodules.h (revision 31451) @@ -148,4 +148,5 @@ #include #include +#include #include #include Index: /trunk/psModules/test/objects/tap_pmSource.c =================================================================== --- /trunk/psModules/test/objects/tap_pmSource.c (revision 31450) +++ /trunk/psModules/test/objects/tap_pmSource.c (revision 31451) @@ -91,5 +91,5 @@ ok(isnan(src->psfMag), "pmSourceAlloc() pmSource->psfMag correctly"); ok(isnan(src->extMag), "pmSourceAlloc() pmSource->extMag correctly"); - ok(isnan(src->errMag), "pmSourceAlloc() pmSource->errMag correctly"); + ok(isnan(src->psfMagErr), "pmSourceAlloc() pmSource->psfMagErr correctly"); ok(isnan(src->apMag), "pmSourceAlloc() pmSource->apMag correctly"); ok(isnan(src->sky), "pmSourceAlloc() pmSource->sky correctly"); Index: /trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_0.c =================================================================== --- /trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_0.c (revision 31450) +++ /trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_0.c (revision 31451) @@ -167,5 +167,5 @@ model->dparams->data.F32[PM_PAR_YPOS] = TEST_BASE_Y_ERR + (float) i; src->psfMag = TEST_BASE_PSF_MAG + (float) i; - src->errMag = TEST_BASE_ERR_MAG + (float) i; + src->psfMagErr = TEST_BASE_ERR_MAG + (float) i; src->peak->flux = TEST_BASE_PEAK_FLUX + (float) i; src->sky = TEST_BASE_SKY + (float) i; @@ -215,6 +215,6 @@ ok(src->psfMag == (TEST_BASE_PSF_MAG + (float) i), "pmSourcesRead_PS1_DEV_0() set src->psfMag correctly (is %.2f, should be %.2f)", src->psfMag, (TEST_BASE_PSF_MAG + (float) i)); - ok(src->errMag == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_0() set src->errMag correctly (is %.2f, should be %.2f)", - src->errMag, (TEST_BASE_ERR_MAG + (float) i)); + ok(src->psfMagErr == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_0() set src->psfMagErr correctly (is %.2f, should be %.2f)", + src->psfMagErr, (TEST_BASE_ERR_MAG + (float) i)); // XXX: Source code always sets src->modelPSF. Is that right? Index: /trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c =================================================================== --- /trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c (revision 31450) +++ /trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c (revision 31451) @@ -196,5 +196,5 @@ model->dparams->data.F32[PM_PAR_YPOS] = TEST_BASE_Y_ERR + (float) i; src->psfMag = TEST_BASE_PSF_MAG + (float) i; - src->errMag = TEST_BASE_ERR_MAG + (float) i; + src->psfMagErr = TEST_BASE_ERR_MAG + (float) i; src->peak->flux = TEST_BASE_PEAK_FLUX + (float) i; src->sky = TEST_BASE_SKY + (float) i; @@ -203,5 +203,5 @@ sources->data[i] = (psPtr *) src; src->psfMag = TEST_BASE_PSF_MAG + (float) i; - src->errMag = TEST_BASE_ERR_MAG + (float) i; + src->psfMagErr = TEST_BASE_ERR_MAG + (float) i; src->sky = TEST_BASE_SKY + (float) i; src->skyErr = TEST_BASE_SKY_ERR + (float) i; @@ -280,6 +280,6 @@ ok(src->psfMag == (TEST_BASE_PSF_MAG + (float) i), "pmSourcesRead_PS1_DEV_1() set src->psfMag correctly (is %.2f, should be %.2f)", src->psfMag, (TEST_BASE_PSF_MAG + (float) i)); - ok(src->errMag == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_1() set src->errMag correctly (is %.2f, should be %.2f)", - src->errMag, (TEST_BASE_ERR_MAG + (float) i)); + ok(src->psfMagErr == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_1() set src->psfMagErr correctly (is %.2f, should be %.2f)", + src->psfMagErr, (TEST_BASE_ERR_MAG + (float) i)); // XXX: Source code always sets src->modelPSF. Is that right? Index: /trunk/psModules/test/objects/tap_pmSourceIO_SMPDATA.c =================================================================== --- /trunk/psModules/test/objects/tap_pmSourceIO_SMPDATA.c (revision 31450) +++ /trunk/psModules/test/objects/tap_pmSourceIO_SMPDATA.c (revision 31451) @@ -167,5 +167,5 @@ model->dparams->data.F32[PM_PAR_YPOS] = TEST_BASE_Y_ERR + (float) i; src->psfMag = TEST_BASE_PSF_MAG + (float) i; - src->errMag = TEST_BASE_ERR_MAG + (float) i; + src->psfMagErr = TEST_BASE_ERR_MAG + (float) i; src->peak->flux = TEST_BASE_PEAK_FLUX + (float) i; src->sky = TEST_BASE_SKY + (float) i; @@ -226,6 +226,6 @@ src->psfMag, (TEST_BASE_PSF_MAG + (float) i)); float tmpF = 0.001 * PS_MIN(999, (1000 * (TEST_BASE_ERR_MAG + (float) i))); - ok(src->errMag == tmpF, "pmSourcesRead_SMPDATA() set src->errMag correctly (is %.2f, should be %.2f)", - src->errMag, tmpF); + ok(src->psfMagErr == tmpF, "pmSourcesRead_SMPDATA() set src->psfMagErr correctly (is %.2f, should be %.2f)", + src->psfMagErr, tmpF); tmpF = PS_MIN(99.0, (TEST_BASE_EXT_MAG + ZERO_POINT)) - ZERO_POINT; ok(src->extMag == tmpF, "pmSourcesRead_SMPDATA() set src->extMag correctly (is %.2f, should be %.2f)", Index: /trunk/psModules/test/objects/tap_pmSourcePhotometry.c =================================================================== --- /trunk/psModules/test/objects/tap_pmSourcePhotometry.c (revision 31450) +++ /trunk/psModules/test/objects/tap_pmSourcePhotometry.c (revision 31451) @@ -103,5 +103,5 @@ ok_float_tol(source->psfMag, fitMag, 0.0002, "source fitMag is %f", source->psfMag); ok_float_tol(source->apMag, apMag, 0.0002, "source apMag is %f", source->apMag); - ok_float(source->errMag, 0.001, "source errMag is %f", source->errMag); + ok_float(source->psfMagErr, 0.001, "source psfMagErr is %f", source->psfMagErr); float refMag = source->apMag;