Changeset 31361 for branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Apr 24, 2011, 10:09:38 PM (15 years ago)

Author:

eugene

Message:

determine the min kron radius = radius for bright PSF stars; mask interpolation used index instead of pixel coordinates in psImageInterpolate; add minKronRadius to pmSource; change errMag to psfMagErr; in growth curve, avoid perfect int radii (they can be inconsistent on + and - due to float precision; add pmGrowthCurveForSources; psfMagErr is always calculated from psfModel; apply ApTrend to all source psf mags and fluxes (not just psf sources); always use the psfModel for PSF_QF,_PERFECT; use BILINEAR to interpolate for aperture mags (BIQUAD implementation is wrong; apply growth curve to apMag & apFlux for all sources; save GrowthCurve with PSF model and read from PSF model

File:

• branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourcePhotometry.c (modified) (9 diffs)

branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourcePhotometry.c

@@ -82 +82 @@
    - 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 +97 @@
 
     int status = false;
-    bool isPSF;
     float x, y;
     float SN;
-    pmModel *model;
 
     source->psfMag    = NAN;
     source->extMag    = NAN;
-    source->errMag    = NAN;
+    source->psfMagErr = NAN;
     source->apMag     = NAN;
     source->apMagRaw  = NAN;
@@ -113 +108 @@
     source->apFluxErr = NAN;
 
-    // we must have a valid model
-    // XXX allow aperture magnitudes for sources without a model
-    model = pmSourceGetModel (&isPSF, source);
-    if (model == NULL) {
-        psTrace ("psModules.objects", 3, "fail mag : no valid model");
+    // XXXXXX review:
+    // Select the 'best' model -- this is used for PSF_QF,_PERFECT & ???. isPSF is true if this
+    // object is a PSF (not extended).  We must have a valid model.  XXX NOTE: allow aperture
+    // magnitudes for sources without a model
+
+    // select the psf model
+    pmModel *modelPSF = source->modelPSF;
+    if (modelPSF == NULL) {
+        psTrace ("psModules.objects", 3, "fail mag : no valid PSF model");
         return false;
     }
 
-    // XXX handle negative flux, low-significance
-    if (model->dparams->data.F32[PM_PAR_I0] > 0) {
-        SN = fabs(model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0]);
-        source->errMag = 1.0 / SN;
+    // get the error on the PSF model magnitude
+    if (modelPSF->dparams->data.F32[PM_PAR_I0] > 0) {
+        SN = fabs(modelPSF->params->data.F32[PM_PAR_I0] / modelPSF->dparams->data.F32[PM_PAR_I0]);
+        source->psfMagErr = 1.0 / SN;
     } else {
         SN = NAN;
-        source->errMag = NAN;
-    }
-    x = model->params->data.F32[PM_PAR_XPOS];
-    y = model->params->data.F32[PM_PAR_YPOS];
+        source->psfMagErr = NAN;
+    }
+    // the source position is used to recenter the aperture for ap photometry
+    x = modelPSF->params->data.F32[PM_PAR_XPOS];
+    y = modelPSF->params->data.F32[PM_PAR_YPOS];
 
     // measure PSF model photometry
-    // XXX TEST: do not use flux scale
-    // XXX NOTE: turn this back on?
-    if (0 && psf->FluxScale) {
-        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
-        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
-        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
-        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
-        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
-        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
-        source->psfMag = -2.5*log10(source->psfFlux);
-    } else {
-        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
-        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
-    }
+    status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, modelPSF);
+    source->psfFluxErr = source->psfFlux * source->psfMagErr;
+
+# if (0)
+    // XXX NOTE: old code to use the flux scale.  test & turn this back on?  if so, need to save with psf model
+    // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
+    double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
+    psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
+    psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
+    source->psfFlux = fluxScale * modelPSF->params->data.F32[PM_PAR_I0];
+    source->psfFluxErr = fluxScale * modelPSF->dparams->data.F32[PM_PAR_I0];
+    source->psfMag = -2.5*log10(source->psfFlux);
+# endif
 
     if (mode == PM_SOURCE_PHOT_PSFONLY) {
@@ -152 +151 @@
     }
 
+    // get the EXT model photometry (all EXT models)
     // if we have a collection of model fits, check if one of them is a pointer to modelEXT
     if (source->modelFits) {
@@ -172 +172 @@
     }
 
-    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
-    // XXX add a flag for "ap_mag is corrected?"
-    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) {
+    // Correct psfMag to match aperture magnitude system (NOTE : Growth curve is already applied to ApTrend)
+    if ((mode & PM_SOURCE_PHOT_APCORR) && psf && psf->ApTrend) {
         // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
         double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
@@ -182 +181 @@
     }
 
-    // measure the contribution of included pixels
+    // measure the contribution of included pixels to the PSF model fit
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (source, model, source->maskObj, maskVal, radius);
+        pmSourcePixelWeight (source, modelPSF, source->maskObj, maskVal, radius);
     }
 
@@ -218 +217 @@
         y += dy;
 
-        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy,
-                              NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) {
+        // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_LANCZOS2)) {
+        // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) {
+        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BILINEAR)) {
             // Not much we can do about it
             psErrorClear();
@@ -233 +233 @@
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAperSource (source, model, flux, variance, mask, maskVal);
+    status = pmSourcePhotometryAperSource (source, modelPSF, flux, variance, mask, maskVal);
     if (!status) {
         psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
@@ -242 +242 @@
     // detection limits (esp near bright neighbors)
     source->apMag = source->apMagRaw;
-    if (isfinite (source->apMag) && isPSF && psf) {
+    if (isfinite (source->apMag) && psf) {
         if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
-            source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius);
-            // XXX correct the apFlux?
+            float apOffset = pmGrowthCurveCorrect (psf->growth, source->apRadius);
+            source->apMag = source->apMagRaw + apOffset;
+            source->apFlux *= pow(10.0, -0.4*apOffset);
+            source->apFluxErr *= pow(10.0, -0.4*apOffset);
         }
     }