Index: branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c (revision 31384) @@ -96,6 +96,8 @@ growth->refRadius = refRadius; growth->maxRadius = maxRadius; - growth->apLoss = 0.0; - growth->fitMag = 0.0; + growth->fitMag = NAN; + growth->apRef = NAN; + growth->apLoss = NAN; + return growth; } Index: branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.c =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.c (revision 31384) @@ -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: branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.h =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.h (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.h (revision 31384) @@ -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; Index: branches/eam_branches/ipp-20110404/psModules/src/objects/pmPeaks.c =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmPeaks.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmPeaks.c (revision 31384) @@ -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: branches/eam_branches/ipp-20110404/psModules/src/objects/pmPeaks.h =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmPeaks.h (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmPeaks.h (revision 31384) @@ -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: branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.c =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.c (revision 31384) @@ -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->psfMagErr = 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,29 +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->mode2 = in->mode2; - source->tmpFlags = in->tmpFlags; - source->imageID = in->imageID; - - source->psfMag = in->psfMag; - source->psfFlux = in->psfFlux; - source->psfFluxErr = in->psfFluxErr; - source->extMag = in->extMag; - source->psfMagErr = in->psfMagErr; - 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->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); @@ -1195,5 +1194,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: branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.h =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.h (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.h (revision 31384) @@ -80,5 +80,5 @@ 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 @@ -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 Index: branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 31384) @@ -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: branches/eam_branches/ipp-20110404/psphot/src/psphotBlendFit.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotBlendFit.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotBlendFit.c (revision 31384) @@ -242,4 +242,7 @@ if (source->mode & PM_SOURCE_MODE_PAIR) continue; + // do not include MOMENTS_FAILURES in the fit + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + // limit selection to some SN limit if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue; Index: branches/eam_branches/ipp-20110404/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotFitSourcesLinear.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotFitSourcesLinear.c (revision 31384) @@ -137,4 +137,7 @@ if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + // do not include MOMENTS_FAILURES in the fit + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + // XXX count saturated stars if (source->mode & PM_SOURCE_MODE_SATSTAR) { @@ -166,13 +169,20 @@ // apply mask? float modelSum = 0.0; + float maskedSum = 0.0; for (int iy = 0; iy < source->modelFlux->numRows; iy++) { for (int ix = 0; ix < source->modelFlux->numCols; ix++) { modelSum += source->modelFlux->data.F32[iy][ix]; + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + maskedSum += source->modelFlux->data.F32[iy][ix]; } } if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit) if (modelSum < 0.8) { - fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n", - source->peak->xf, source->peak->yf, modelSum, source->peak->rawFlux); + fprintf (stderr, "low-sig model @ %f, %f (%f masked sum, %f sum, %f peak)\n", + source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux); + } + if (maskedSum < 0.5) { + fprintf (stderr, "worrying model @ %f, %f (%f masked sum, %f sum, %f peak)\n", + source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux); } Index: branches/eam_branches/ipp-20110404/psphot/src/psphotGuessModels.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotGuessModels.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotGuessModels.c (revision 31384) @@ -160,8 +160,4 @@ pmSource *source = sources->data[i]; - if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) { - fprintf (stderr, "moment failure\n"); - } - // this is used to mark sources for which the model is measured. We check later that // all are used. @@ -169,4 +165,5 @@ // skip non-astronomical objects (very likely defects) + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; @@ -196,4 +193,5 @@ } +# if (0) if (source->mode & PM_SOURCE_MODE_SATSTAR) { fprintf (stderr, "satstar: %f,%f vs %f,%f : %c\n", @@ -202,4 +200,5 @@ (useMoments ? 'T' : 'F')); } +# endif // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo) Index: branches/eam_branches/ipp-20110404/psphot/src/psphotMakeGrowthCurve.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotMakeGrowthCurve.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotMakeGrowthCurve.c (revision 31384) @@ -50,8 +50,8 @@ # endif + psLogMsg ("psphot", PS_LOG_MINUTIA, "built growth curve: %f sec\n", psTimerMark ("psphot.growth")); + float offset = pmGrowthCurveCorrect (psf->growth, PSF_APERTURE); - fprintf (stderr, "correction from %f to %f: %f mags\n", PSF_APERTURE, REF_RADIUS, offset); - - psLogMsg ("psphot", PS_LOG_MINUTIA, "built growth curve: %f sec\n", psTimerMark ("psphot.growth")); + psLogMsg ("psphot", PS_LOG_DETAIL, "correction from %f to %f: %f mags\n", PSF_APERTURE, REF_RADIUS, offset); return true; Index: branches/eam_branches/ipp-20110404/psphot/src/psphotRadialApertures.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotRadialApertures.c (revision 31381) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotRadialApertures.c (revision 31384) @@ -95,5 +95,7 @@ if (source->type == PM_SOURCE_TYPE_SATURATED) continue; if (source->mode & PM_SOURCE_MODE_DEFECT) continue; - if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + // XXX measure radial apertures even for saturated stars + // if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; // limit selection to some SN limit @@ -257,9 +259,12 @@ float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean)); + // XXX report the total flux or the mask-corrected flux? // flux->data.F32[i] = SBmean * Area; + // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix; + + fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]); fluxStd->data.F32[i] = SBstdv * Area; - fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix; - - // fill->data.F32[i] /= Area; + fill->data.F32[i] /= Area; + psTrace ("psphot", 5, "radial bins: %3d %5.1f : %8.1f +/- %7.2f : %8.1f +/- %8.1f : %4.2f %6.1f\n", i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], SBmean, SBstdv, fill->data.F32[i], Area);