Changeset 29060 for branches/sc_branches/trunkTest/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

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

branches/sc_branches/trunkTest/psModules/src/objects/pmSourcePhotometry.c

@@ -22 +22 @@
 #include "pmFPA.h"
 #include "pmFPAMaskWeight.h"
+
+#include "pmTrend2D.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
 #include "pmSpan.h"
+#include "pmFootprintSpans.h"
 #include "pmFootprint.h"
 #include "pmPeaks.h"
 #include "pmMoments.h"
-#include "pmGrowthCurve.h"
-#include "pmResiduals.h"
-#include "pmTrend2D.h"
+#include "pmModelFuncs.h"
+#include "pmModel.h"
+#include "pmModelUtils.h"
+#include "pmModelClass.h"
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
+#include "pmSource.h"
+#include "pmSourceFitModel.h"
 #include "pmPSF.h"
-#include "pmModel.h"
-#include "pmSource.h"
-#include "pmModelClass.h"
+#include "pmPSFtry.h"
+
 #include "pmSourcePhotometry.h"
 
@@ -66 +76 @@
 
 // XXX masked region should be (optionally) elliptical
-bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal)
+bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -122 +132 @@
         for (int i = 0; i < source->modelFits->n; i++) {
             pmModel *model = source->modelFits->data[i];
+            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
             status = pmSourcePhotometryModel (&model->mag, NULL, model);
             if (model == source->modelEXT) foundEXT = true;
@@ -145 +156 @@
     // measure the contribution of included pixels
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
+        pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal);
     }
 
     // measure the contribution of included pixels
     if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
-        pmSourceMeasureDiffStats (source, maskVal);
+        pmSourceMeasureDiffStats (source, maskVal, markVal);
     }
 
@@ -191 +202 @@
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAper  (&source->apMag, model, flux, mask, maskVal);
+    status = pmSourcePhotometryAper  (&source->apMagRaw, model, flux, mask, maskVal);
     if (!status) {
         psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
@@ -199 +210 @@
     // if the aper mag is NAN, the flux < 0.  this can happen for sources near the
     // detection limits (esp near bright neighbors)
+    source->apMag = source->apMagRaw;
     if (isfinite (source->apMag) && isPSF && psf) {
         if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
-            source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
+            source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius);
         }
         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;
         }
@@ -257 +272 @@
     PS_ASSERT_PTR_NON_NULL(image, false);
     PS_ASSERT_PTR_NON_NULL(mask, false);
-    PS_ASSERT_PTR_NON_NULL(model, false);
+
+    if (DO_SKY) {
+        PS_ASSERT_PTR_NON_NULL(model, false);
+    }
 
     float apSum = 0;
@@ -271 +289 @@
     psF32 **imData = image->data.F32;
     psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA;
+    int nAperPix = 0;
 
     // measure apMag
-    for (int ix = 0; ix < image->numCols; ix++) {
-        for (int iy = 0; iy < image->numRows; iy++) {
+    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);
         }
     }
@@ -290 +311 @@
 
 // return source aperture magnitude
-bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *mask, psImageMaskType maskVal)
-{
-    PS_ASSERT_PTR_NON_NULL(pixWeight, false);
+bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal)
+{
+    PS_ASSERT_PTR_NON_NULL(pixWeightNotBad, false);
+    PS_ASSERT_PTR_NON_NULL(pixWeightNotPoor, false);
     PS_ASSERT_PTR_NON_NULL(mask, false);
     PS_ASSERT_PTR_NON_NULL(model, false);
 
     float modelSum = 0;
-    float validSum = 0;
+    float notBadSum = 0;
+    float notPoorSum = 0;
     float sky = 0;
     float value;
@@ -305 +328 @@
     int dY, DY, NY;
 
-    *pixWeight = 0.0;
+    *pixWeightNotBad = 0.0;
+    *pixWeightNotPoor = 0.0;
 
     // we only care about the value of the object model, not the local sky
@@ -345 +369 @@
 
             // for the full model, add all points
-            value = model->modelFunc (NULL, params, coord) - sky;
+            value = fabs(model->modelFunc (NULL, params, coord) - sky);
             modelSum += value;
 
             // include count only the unmasked pixels within the image area
-            if (mx < 0)
-                continue;
-            if (my < 0)
-                continue;
-            if (mx >= NX)
-                continue;
-            if (my >= NY)
-                continue;
-            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)
-                continue;
-
-            validSum += value;
+            if (mx < 0) continue;
+            if (my < 0) continue;
+            if (mx >= NX) continue;
+            if (my >= NY) continue;
+
+            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) {
+                notBadSum += value;
+            }
+            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & ~markVal)) {
+                notPoorSum += value;
+            }
         }
     }
     psFree (coord);
 
-    if (validSum <= 0)
-        return false;
-
-    *pixWeight = validSum / modelSum;
+    *pixWeightNotBad  = notBadSum  / modelSum;
+    *pixWeightNotPoor = notPoorSum / modelSum;
+
     return (true);
 }
@@ -374 +396 @@
 # define FLUX_LIMIT 3.0
 
-// return source aperture magnitude
-bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
+// measure stats that may be using in difference images for distinguishing real sources from bad residuals
+bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -399 +421 @@
     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)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & markVal) {
+                continue;
+            }
             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
                 nMask ++;