Changeset 29546 for trunk/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Oct 25, 2010, 2:54:22 PM (16 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20100823

File:

• trunk/psModules/src/objects/pmSourcePhotometry.c (modified) (15 diffs)

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -62 +62 @@
 
 /**
-    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 +89 @@
     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 +118 @@
     // 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 +132 @@
     }
 
+    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 +157 @@
 
     // 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 +191 @@
     // 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 +210 @@
     } 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 +227 @@
         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 +247 @@
 
 /*
-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 +281 @@
 
 // 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 +308 @@
     }
 
-    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 +450 @@
     }
 
+    // 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 +467 @@
             }
 
-            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 +648 @@
 
             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 +726 @@
 
             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");
             }