Changeset 6556 for branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Mar 8, 2006, 5:14:23 PM (20 years ago)

Author:

magnier

Message:

major rework of objects code

File:

• branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c (modified) (9 diffs)

branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c

@@ -1 @@
-/** @file  pmSourcePhotometry.C
+/** @file  pmSourcePhotometry.c
  *
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-07 06:33:35 $
+ *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-03-09 03:14:23 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -10 +10 @@
  */
 
-// XXX EAM : the apMag should only be calculated for the brighter sources?
-// XXX EAM : SN limit set by user?
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include "pslib.h"
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmGrowthCurve.h"
+#include "pmModel.h"
+#include "pmPSF.h"
+#include "pmSource.h"
+#include "pmModelGroup.h"
+#include "pmSourcePhotometry.h"
+
+static float AP_MIN_SN = 0.0;
+
+bool pmSourceMagnitudesInit (psMetadata *config)
+{
+
+    bool status;
+
+    float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN");
+    if (status) {
+        AP_MIN_SN = limit;
+    }
+    return true;
+}
+
+/**
+    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 (how??)
+    - psfMag : all sources with non-NULL modelPSF
+             : is optionally corrected for aperture residual if:
+        - the source is a STAR (PSF)
+        - the option is selected (how??)
+       
+    - extMag : all sources with non-NULL modelEXT 
+**/
+
 // XXX EAM : masked region should be (optionally) elliptical
+// XXX curve of growth is corrected to
 pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius)
 {
@@ -21 +63 @@
     float rflux;
     float radius;
+    float SN;
     pmModel *model;
 
     switch (source->type) {
-    case PM_SOURCE_STAR:
+    case PM_SOURCE_TYPE_STAR:
         model = source->modelPSF;
         if (model == NULL)
@@ -32 +75 @@
         break;
 
-    case PM_SOURCE_EXTENDED:
+    case PM_SOURCE_TYPE_EXTENDED:
         model = source->modelEXT;
         if (model == NULL)
@@ -44 +87 @@
     }
 
+    if (model->dparams->data.F32[1] > 0) {
+        SN = model->params->data.F32[1] / model->dparams->data.F32[1];
+        source->errMag = 1.0 / SN;
+    } else {
+        SN = 0.0;
+        source->errMag = 0.0;
+    }
     x = model->params->data.F32[2];
     y = model->params->data.F32[3];
 
     // replace source flux
+    // XXX test to see if source has been subtracted?
     pmModelAdd (source->pixels, source->mask, model, false, false);
 
-    // set aperture mask circle of PSF_FIT_RADIUS
-    psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED);
-
-    // measure object photometry
-    status = pmSourcePhotometry (&source->fitMag, &source->apMag, model, source->pixels, source->mask);
-
-    // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS
+    // set aperture mask circle to scaled radius
+    psImageKeepCircle (source->mask, x, y, model->radius, "OR", PM_SOURCE_MASK_MARKED);
+
+    // measure object model photometry
+    status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
+    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+
+    // measure object aperture photometry
+    if (SN > AP_MIN_SN) {
+        status = pmSourcePhotometryAper  (&source->apMag,  model, source->pixels, source->mask);
+    } else {
+        source->apMag = 99.0;
+    }
+
+    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
     if (isPSF && (psf != NULL)) {
-        if (psf->growth != NULL) {
-            source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius);
-        }
-
-        rflux   = pow (10.0, 0.4*source->fitMag);
-        source->apMag  -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux;
-        source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0);
+        if (SN > AP_MIN_SN) {
+            if (psf->growth != NULL) {
+                source->apMag += pmGrowthCurveCorrect (psf->growth, radius);
+            }
+            rflux   = pow (10.0, 0.4*source->psfMag);
+            source->apMag  -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux;
+        }
+        source->psfMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0);
     }
 
     // unmask aperture
-    psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED);
+    psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PM_SOURCE_MASK_MARKED);
 
     // subtract object, leave local sky
@@ -85 +145 @@
 */
 
-// XXX split into ap, psf, ext mags?
-bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask)
+// return source model magnitude
+bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
+{
+
+    float fitSum = 0;
+    *fitMag = 99.0;
+
+    if (model == NULL) {
+        return false;
+    }
+
+    // measure fitMag
+    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
+    fitSum = modelFluxFunc (model->params);
+    if (fitSum <= 0)
+        return false;
+    if (!isfinite(fitSum))
+        return false;
+    *fitMag = -2.5*log10(fitSum);
+
+    return (true);
+}
+
+// return source aperture magnitude
+bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask)
+{
+
+    float apSum = 0;
+    *apMag = 99.0;
+
+    if (model == NULL) {
+        return false;
+    }
+
+    float sky = model->params->data.F32[0];
+
+    // measure apMag
+    for (int ix = 0; ix < image->numCols; ix++) {
+        for (int iy = 0; iy < image->numRows; iy++) {
+            if (mask->data.U8[iy][ix])
+                continue;
+            apSum += image->data.F32[iy][ix] - sky;
+        }
+    }
+    if (apSum <= 0)
+        return false;
+    *apMag = -2.5*log10(apSum);
+
+    return (true);
+}
+
+# if (0)
+    // XXX split into ap, psf, ext mags?
+    bool pmSourcePhotometry (float *fitMag, float *apMag, pmModel *model, psImage *image, psImage *mask)
 {
 
@@ -96 +208 @@
     *obsMag = 99.0;
 
+    // measure fitMag
     pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
     fitSum = modelFluxFunc (model->params);
@@ -104 +217 @@
     *fitMag = -2.5*log10(fitSum);
 
+    // measure apMag
     for (int ix = 0; ix < image->numCols; ix++) {
         for (int iy = 0; iy < image->numRows; iy++) {
@@ -117 +231 @@
     return (true);
 }
+# endif
 
 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj)