Changeset 31606 for branches/czw_branch/20110406/psphot/src/psphotSourceStats.c

Timestamp:

Jun 6, 2011, 1:56:22 PM (15 years ago)

Author:

watersc1

Message:

merge from trunk before working on new stacking code

Location:

branches/czw_branch/20110406

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotSourceStats.c (modified) (16 diffs)

branches/czw_branch/20110406/psphot/src/psphotSourceStats.c

@@ -98 +98 @@
 
     // generate the array of sources, define the associated pixel
+    bool firstPass = false;
     if (!detections->newSources) {
         detections->newSources = psArrayAllocEmpty (peaks->n);
+        firstPass = true;
     }
     sources = detections->newSources;
@@ -125 +127 @@
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
             fprintf (stderr, "moment failure\n");
+        }
+
+        if (firstPass) {
+            source->mode2 |= PM_SOURCE_MODE2_PASS1_SRC;
         }
 
@@ -159 +165 @@
         RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     }
+    float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON");
+    if (!status) {
+        MIN_KRON_RADIUS = 0.75*SIGMA;
+    }
 
     // threaded measurement of the source magnitudes
@@ -186 +196 @@
             PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32);
 
             PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
@@ -215 +226 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[7];
+                scalar = job->args->data[8];
                 Nmoments += scalar->data.S32;
-                scalar = job->args->data[8];
+                scalar = job->args->data[9];
                 Nfail += scalar->data.S32;
-                scalar = job->args->data[9];
+                scalar = job->args->data[10];
                 Nfaint += scalar->data.S32;
             }
@@ -231 +242 @@
 
     psphotVisualShowMoments (sources);
+
+    // clear the mark bits
+    // psImageMaskPixels (readout->mask, "AND", PS_NOT_IMAGE_MASK(markVal));
 
     if (detections->allSources) {
@@ -373 +387 @@
     float RADIUS                    = PS_SCALAR_VALUE(job->args->data[3],F32);
     float SIGMA                     = PS_SCALAR_VALUE(job->args->data[4],F32);
-
-    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
+    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[5],F32);
+
+    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
 
     // if no valid pixels, massive swing or very large Mrf, likely saturated source,
@@ -431 +446 @@
         if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) {
             // measure basic source moments (no S/N clipping on input pixels)
-            status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
+            status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
         } else { 
             // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but
@@ -447 +462 @@
 
             psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
-            status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal);
+            status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
             source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
         }
@@ -456 +471 @@
             continue;
         }
+
+        // XXX test of masking neighbors when measureing moments (does this fail?)
+        // psEllipseMoments moments = {source->moments->Mxx, source->moments->Myy, source->moments->Mxy};
+        // psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
+        // psImageMaskCircle (source->maskView, source->peak->x, source->peak->y, 3.0*axes.major, "OR", markVal);
+
         Nmoments ++;
 
@@ -492 +513 @@
     float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
     float Sout[NSIGMA];
+    float Rmin[NSIGMA];
     int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
 
@@ -513 +535 @@
 
             // measure basic source moments (no S/N clipping on input pixels)
-            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, maskVal);
+            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, 0.0, maskVal);
         }
 
@@ -524 +546 @@
         pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
         psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
+
+        Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);
 
 #if 0
@@ -549 +573 @@
     // we are looking for sigma for which Sout = 0.65 (or some other value)
     int Nstars = 0;
+    float minKronRadius = NAN;
     float Sigma = NAN;
     float minS = Sout[0];
@@ -565 +590 @@
         Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
         Nstars = 0.5*(Nout[i] + Nout[i+1]);
+        minKronRadius = Rmin[i] + (0.65 - Sout[i])*(Rmin[i+1] - Rmin[i])/(Sout[i+1] - Sout[i]);
     }
     psAssert (isfinite(Sigma), "did we miss a case?");
@@ -575 +601 @@
     psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
     psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars);
-
-    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_MIN_KRON", PS_META_REPLACE, "minimum Kron Radius", minKronRadius);
+
+    psLogMsg ("psphot", 3, "using window function with sigma = %f (Min Kron Radius = %f, Nstars = %d)\n", Sigma, minKronRadius, Nstars);
     return true;
 }