Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 29004) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 29546) @@ -62,18 +62,19 @@ /** - this function is used to calculate the three defined source magnitudes: - - 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 + this function is used to calculate the three defined source magnitudes: + - 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 **/ // XXX masked region should be (optionally) elliptical +// if mode is PM_SOURCE_PHOT_PSFONLY, we skip all other magnitudes bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal) { @@ -88,8 +89,11 @@ pmModel *model; - source->psfMag = NAN; - source->extMag = NAN; - source->errMag = NAN; - source->apMag = NAN; + source->psfMag = NAN; + source->extMag = NAN; + source->errMag = NAN; + source->apMag = NAN; + source->apMagRaw = NAN; + source->apFlux = NAN; + source->apFluxErr = NAN; // we must have a valid model @@ -114,4 +118,5 @@ // 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 @@ -127,4 +132,8 @@ } + if (mode == PM_SOURCE_PHOT_PSFONLY) { + return true; + } + // if we have a collection of model fits, check if one of them is a pointer to modelEXT if (source->modelFits) { @@ -148,4 +157,5 @@ // 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) { // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center @@ -181,5 +191,7 @@ // if we are measuring aperture photometry and applying the growth correction, // we need to shift the flux in the selected pixels (but not the mask) - psImage *flux = NULL, *mask = NULL; // Star flux and mask images, to photometer + psImage *flux = NULL; + psImage *variance = NULL; + psImage *mask = NULL; // Star flux and mask images, to photometer if (mode & PM_SOURCE_PHOT_INTERP) { float dx = 0.5 - x + (int)x; @@ -198,9 +210,10 @@ } else { flux = source->pixels; + variance = source->variance; mask = source->maskObj; } // measure object aperture photometry - status = pmSourcePhotometryAper (&source->apMagRaw, model, flux, mask, maskVal); + status = pmSourcePhotometryAperSource (source, model, flux, variance, mask, maskVal); if (!status) { psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag"); @@ -214,4 +227,5 @@ 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) { @@ -233,7 +247,7 @@ /* -aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) -(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) -(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) + aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) + (aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) + (aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) */ @@ -267,5 +281,22 @@ // return source aperture magnitude -bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask, psImageMaskType maskVal) +bool pmSourcePhotometryAperSource (pmSource *source, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(mask, false); + + if (DO_SKY) { + PS_ASSERT_PTR_NON_NULL(model, false); + } + + bool status; + status = pmSourcePhotometryAper(&source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal); + + return status; +} + +// return source aperture magnitude +bool pmSourcePhotometryAper (float *apMag, float *apFluxOut, float *apFluxErr, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(apMag, false); @@ -277,34 +308,34 @@ } - float apSum = 0; float sky = 0; - *apMag = NAN; + float apFlux = 0; + float apFluxVar = 0; if (DO_SKY) { sky = model->params->data.F32[PM_PAR_SKY]; - } else { - sky = 0; } psF32 **imData = image->data.F32; psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA; - int nAperPix = 0; - - // measure apMag + psF32 **varData = (variance) ? variance->data.F32 : image->data.F32; // if variance is not supplied, assume gain of 1.0, no read noise + + // measure apFlux and apFluxVar, save apMag if not NAN + // XXX note that these fluxes/mags are uncorrected for masked pixels + // XXX raise a bit if the aperture has a masked pixel (not marked)? for (int iy = 0; iy < image->numRows; iy++) { for (int ix = 0; ix < image->numCols; ix++) { - if (mkData[iy][ix] & maskVal) - continue; - apSum += imData[iy][ix] - sky; - nAperPix ++; - // fprintf (stderr, "aper: %d %d %f %f %f\n", ix, iy, sky, imData[iy][ix], apSum); - } - } - if (apSum <= 0) { + if (mkData[iy][ix] & maskVal) continue; + apFlux += imData[iy][ix] - sky; + apFluxVar += varData[iy][ix]; + } + } + if (apFluxOut) *apFluxOut = apFlux; + if (apFluxErr) *apFluxErr = sqrt(apFluxVar); + + if (apFlux <= 0) { *apMag = NAN; - return true; - } - - *apMag = -2.5*log10(apSum); + } else { + *apMag = -2.5*log10(apFlux); + } return true; } @@ -419,7 +450,13 @@ } + // NOTE: until 2010.10.01, these measurements included a 3sigma-per-pixel significance + // this followed what we understood as the definition given to us + // by Armin, but it always seemed a poor idea -- a faint source is unlikely to have any 3sigma pixels. + // changed to remove the per-pixel filter. + for (int iy = 0; iy < flux->numRows; iy++) { for (int ix = 0; ix < flux->numCols; ix++) { // only count up the stats in the unmarked region (ie, the aperture) + // skip the marked pixels; these are not relevant if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & markVal) { continue; @@ -430,14 +467,12 @@ } - float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - - if (SN > +FLUX_LIMIT) { + float value = flux->data.F32[iy][ix]; + + if (value > 0.0) { nGood ++; - fGood += fabs(flux->data.F32[iy][ix]); - } - - if (SN < -FLUX_LIMIT) { + fGood += fabs(value); + } else { nBad ++; - fBad += fabs(flux->data.F32[iy][ix]); + fBad += fabs(value); } } @@ -613,14 +648,14 @@ switch (term) { - case 0: + case 0: factor = 1; break; - case 1: + case 1: factor = xi + Pi->col0; break; - case 2: + case 2: factor = yi + Pi->row0; break; - default: + default: psAbort("invalid term for pmSourceWeight"); } @@ -691,14 +726,14 @@ switch (term) { - case 0: + case 0: factor = 1; break; - case 1: + case 1: factor = xi + Pi->col0; break; - case 2: + case 2: factor = yi + Pi->row0; break; - default: + default: psAbort("invalid term for pmSourceWeight"); }