Changeset 31606 for branches/czw_branch/20110406/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Jun 6, 2011, 1:56:22 PM (15 years ago)

Author:

watersc1

Message:

merge from trunk before working on new stacking code

Location:

branches/czw_branch/20110406

Files:

• . (modified) (1 prop)
• psModules/src/objects/pmSourcePhotometry.c (modified) (14 diffs)

branches/czw_branch/20110406/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 rflux;
     float SN;
-    pmModel *model;
 
     source->psfMag    = NAN;
     source->extMag    = NAN;
-    source->errMag    = NAN;
+    source->psfMagErr = NAN;
     source->apMag     = NAN;
     source->apMagRaw  = NAN;
@@ -114 +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) {
@@ -153 +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) {
@@ -173 +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);
         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 +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();
@@ -232 +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");
@@ -241 +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?
-        }
-        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 +304 @@
 
     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 +326 @@
     float apFlux = 0;
     float apFluxVar = 0;
+    int nPix = 0;
 
     if (DO_SKY) {
@@ -345 +344 @@
             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 +414 @@
     psImageMaskType maskBad = maskVal;
     maskBad &= ~maskSuspect;
+
+    psImageMaskType maskPoor = maskVal | maskSuspect;
 
     // measure modelSum and validSum.  this function is applied to a sources' subimage.  the
@@ -446 +450 @@
 
             // 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 +491 @@
 
     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");
     }