Changeset 29060 for branches/sc_branches/trunkTest/psModules/src/objects/pmSourceMoments.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/pmSourceMoments.c (modified) (11 diffs)

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

@@ -22 +22 @@
 #include <strings.h>
 #include <pslib.h>
+
 #include "pmHDU.h"
 #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 "pmResiduals.h"
-#include "pmGrowthCurve.h"
-#include "pmTrend2D.h"
-#include "pmPSF.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"
 
@@ -54 +61 @@
 # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
 
+static bool beVerbose = false;
+void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
+
 bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
 {
@@ -61 +71 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    // use sky from moments if defined, 0.0 otherwise
+    // use sky from moments if defined, 0.0 otherwise 
+
+    // XXX this value comes from the sky model at the source center, and tends to over-estimate
+    // the sky in the vicinity of bright sources.  we are better off assuming the model worked
+    // well:
     psF32 sky = 0.0;
     if (source->moments == NULL) {
-        source->moments = pmMomentsAlloc();
-    } else {
-        sky = source->moments->Sky;
-    }
+      source->moments = pmMomentsAlloc();
+    }
+    // XXX if (source->moments == NULL) {
+    // XXX     source->moments = pmMomentsAlloc();
+    // XXX } else {
+    // XXX     sky = source->moments->Sky;
+    // XXX }
 
     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
@@ -131 +148 @@
             psF32 wDiff = *vWgt;
 
-            // skip pixels below specified significance level.  this is allowed, but should be
-            // avoided -- the over-weights the wings of bright stars compared to those of faint
-            // stars.
-            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-            // if (pDiff < 0) continue; // XXX : MWV says I should include < 0.0 valued points...
+            // skip pixels below specified significance level.  for a PSFs, this
+            // over-weights the wings of bright stars compared to those of faint stars.
+            // for the estimator used for extended source analysis (where the window
+            // function is allowed to be arbitrarily large), we need to clip to avoid
+            // negative second moments.
+            if (PS_SQR(pDiff) < minSN2*wDiff) continue; // 
+            if ((minSN > 0.0) && (pDiff < 0)) continue; // 
 
             // Apply a Gaussian window function.  Be careful with the window function.  S/N
@@ -197 +216 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
+    psF32 RF = 0.0;
+    psF32 RH = 0.0;
+    psF32 RS = 0.0;
     psF32 XX = 0.0;
     psF32 XY = 0.0;
@@ -244 +266 @@
             if (r > radius) continue;
 
-            psF32 pDiff = *vPix - sky;
+            psF32 fDiff = *vPix - sky;
+            psF32 pDiff = fDiff;
             psF32 wDiff = *vWgt;
 
@@ -257 +280 @@
             if (sigma > 0.0) {
                 // XXX a lot of extra flops; can we do pre-calculate?
-                psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
+                // XXX we were re-calculating r2 (maybe the compiler caught this?)
+                // psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
+                psF32 z  = r2 * rsigma2;
                 assert (z >= 0.0);
                 psF32 weight  = exp(-z);
@@ -266 +291 @@
 
             Sum += pDiff;
+
+# if (1)
+# if (0)
+            if (fDiff < 0) continue;
+# endif
+            psF32 rf = r * fDiff;
+            psF32 rh = sqrt(r) * fDiff;
+            psF32 rs = fDiff;
+# else
+            psF32 rf = r * pDiff;
+            psF32 rh = sqrt(r) * pDiff;
+            psF32 rs = pDiff;
+# endif
 
             psF32 x = xDiff * pDiff;
@@ -284 +322 @@
             psF32 xyyy = xDiff * yyy / r2;
             psF32 yyyy = yDiff * yyy / r2;
+
+            RF  += rf;
+            RH  += rh;
+            RS  += rs;
 
             XX  += xx;
@@ -302 +344 @@
     }
 
+    source->moments->Mrf = RF/RS;
+    source->moments->Mrh = RH/RS;
+
     source->moments->Mxx = XX/Sum;
     source->moments->Mxy = XY/Sum;
@@ -317 +362 @@
     source->moments->Myyyy = YYYY/Sum;
 
-    // if (source->moments->Mxx < 0) {
-    //  fprintf (stderr, "error: neg second moment??\n");
-    // }
-    // if (source->moments->Myy < 0) {
-    //  fprintf (stderr, "error: neg second moment??\n");
-    // }
-
-    psTrace ("psModules.objects", 4, "Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
+    // Calculate the Kron magnitude (make this block optional?)
+    // float radKron = 2.5*source->moments->Mrf;
+    float radKinner = 1.0*source->moments->Mrf;
+    float radKron   = 2.5*source->moments->Mrf;
+    float radKouter = 4.0*source->moments->Mrf;
+
+    int nKronPix = 0;
+    Sum = Var = 0.0;
+    float SumInner = 0.0;
+    float SumOuter = 0.0;
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+        psF32 yDiff = row - yCM;
+        if (fabs(yDiff) > radKron) continue;
+
+        psF32 *vPix = source->pixels->data.F32[row];
+        psF32 *vWgt = source->variance->data.F32[row];
+        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+
+        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
+            if (vMsk) {
+                if (*vMsk & maskVal) {
+                    vMsk++;
+                    continue;
+                }
+                vMsk++;
+            }
+            if (isnan(*vPix)) continue;
+
+            psF32 xDiff = col - xCM;
+            if (fabs(xDiff) > radKron) continue;
+
+            // radKron is just a function of (xDiff, yDiff)
+            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            psF32 r  = sqrt(r2);
+
+            psF32 pDiff = *vPix - sky;
+            psF32 wDiff = *vWgt;
+
+            // skip pixels below specified significance level.  this is allowed, but should be
+            // avoided -- the over-weights the wings of bright stars compared to those of faint
+            // stars.
+            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
+
+            if (r < radKron) {
+                Sum += pDiff;
+                Var += wDiff;
+                nKronPix ++;
+                // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
+            }
+
+            if ((r > radKinner) && (r < radKron)) {
+                SumInner += pDiff;
+            }
+            if ((r > radKron)  && (r < radKouter)) {
+                SumOuter += pDiff;
+            }
+        }
+    }
+    source->moments->KronFlux = Sum;
+    source->moments->KronFinner = SumInner;
+    source->moments->KronFouter = SumOuter;
+    source->moments->KronFluxErr = sqrt(Var);
+
+    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
+             source->moments->Mrf,   source->moments->KronFlux, 
              source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
              source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,