Changeset 29230 for branches/eam_branches/ipp-20100823/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Sep 24, 2010, 2:49:53 PM (16 years ago)

Author:

eugene

Message:

added apFlux and apFluxErr to pmSource; added AP_FLUX, AP_FLUX_SIG to output in CMF_DV2

File:

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

branches/eam_branches/ipp-20100823/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
 **/
 
@@ -88 +88 @@
     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->apRadius  = NAN;
+    source->apFlux    = NAN;
+    source->apFluxErr = NAN;
 
     // we must have a valid model
@@ -182 +186 @@
     // 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;
@@ -199 +205 @@
     } 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");
@@ -215 +222 @@
         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) {
@@ -234 +242 @@
 
 /*
-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)
 
 */
@@ -268 +276 @@
 
 // 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);
@@ -278 +303 @@
     }
 
-    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
-    // XXX save the apFlux and apFluxErr
+    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;
 }
@@ -617 +639 @@
 
             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");
             }
@@ -695 +717 @@
 
             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");
             }