Changeset 32202

Timestamp:

Aug 29, 2011, 2:54:38 AM (15 years ago)

Author:

eugene

Message:

window for centroid, robustifying kron

File:

• branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceMoments.c (modified) (4 diffs)

branches/eam_branches/ipp-20110710/psModules/src/objects/pmSourceMoments.c

@@ -80 +80 @@
     }
 
+    if (source->moments->nPixels != 0) {
+        fprintf (stderr, "remeasure moments: %f,%f\n", source->peak->xf, source->peak->yf);
+    }
+
     float Sum = 0.0;
     float Var = 0.0;
@@ -96 +100 @@
     // (int) so they can be used in the image index below.
 
-    // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to 
-    // get an unbiased (but probably noisy) centroid
-    if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
-        return false;
-    }
-    // second pass applies the Gaussian window and uses the centroid from the first pass
-    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
+    // XXX // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to 
+    // XXX // get an unbiased (but probably noisy) centroid
+    // XXX if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
+    // XXX      return false;
+    // XXX }
+    // XXX // second pass applies the Gaussian window and uses the centroid from the first pass
+    // XXX if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
+    // XXX      return false;
+    // XXX }
+
+    // If we use a large radius for the centroid, it will be biased by any neighbors.  The flux
+    // of any object drops pretty quickly outside 1-2 sigmas.  (The exception is bright
+    // saturated stars, for which we need to use a very large radius here)
+    if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
         return false;
     }
@@ -318 +329 @@
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
 
+# define WEIGHTED 1
+# if (WEIGHTED)
+            float z = r2 * rsigma2 / 4.0;
+            assert (z >= 0.0);
+            float weight  = exp(-z);
+# endif
+
             float r  = sqrt(r2);
             if (r < radKron) {
+# if (WEIGHTED)
+                Sum += pDiff*weight;
+                Var += wDiff*weight;
+# else
                 Sum += pDiff;
                 Var += wDiff;
+# endif
                 nKronPix ++;
                 // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
@@ -350 +373 @@
     source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
     source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
+
+    // XXX not sure I should save this here...
+    source->moments->KronFluxPSF    = source->moments->KronFlux;
+    source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+    source->moments->KronRadiusPSF  = source->moments->Mrf;
 
     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",