Changeset 31154 for trunk/psphot/src/psphotSourceStats.c

Timestamp:

Apr 4, 2011, 1:12:39 PM (15 years ago)

Author:

eugene

Message:

use the smoothed image for footprint culling; use a more stringent culling for saturated stars; more tweaks to get sat stars to be fitted; various updates to psphotStack; unify psphotImageLoop varients; be a bit careful about number of stars in a PSF clump

Location:

trunk/psphot

Files:

• . (modified) (2 props)
• src/psphotSourceStats.c (modified) (12 diffs)

trunk/psphot

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

trunk/psphot/src/psphotSourceStats.c

@@ -122 +122 @@
         // allocate space for moments
         source->moments = pmMomentsAlloc();
+
+        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
+            fprintf (stderr, "moment failure\n");
+        }
 
         // allocate image, weight, mask arrays for each peak (square of radius OUTER)
@@ -361 +365 @@
 
     bool status = false;
-    float BIG_RADIUS;
     psScalar *scalar = NULL;
 
@@ -373 +376 @@
     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);
+
+    // if no valid pixels, massive swing or very large Mrf, likely saturated source,
+    // try a much larger box
+    float BIG_RADIUS = 3.0*RADIUS;
+    float BIG_SIGMA  = 3.0*SIGMA;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -390 +398 @@
 
         // skip faint sources for moments measurement
-        if (source->peak->SN < MIN_SN) {
+        if (sqrt(source->peak->detValue) < MIN_SN) {
             source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN;
             Nfaint++;
@@ -416 +424 @@
         }
 
-        // measure basic source moments (no S/N clipping on input pixels)
-        status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
-        // XXX moments / aperture test:
-        if (0) {
-            // clear the mask bit and set the circular mask pixels
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, 4.0*source->moments->Mrf, "OR", markVal);
-            
-            float apMag = NAN;
-            pmSourcePhotometryAper (&apMag, NULL, NULL, NULL, source->pixels, source->variance, source->maskObj, maskVal);
-            fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux));
-
-            // clear the mask bit
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-        }           
-        if (status) {
-            Nmoments ++;
-            continue;
-        }
-
-        // if no valid pixels, or massive swing, likely saturated source,
-        // try a much larger box
-        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
-        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
-        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, 0.0, maskVal);
-        if (status) {
-            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
-            Nmoments ++;
-            continue;
-        }
-
-        source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
-        Nfail ++;
-        psErrorClear();
-        continue;
+        // XXX how can this be set: this is raised in pmSourceRoughClass, called after this function?
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
+            fprintf (stderr, "satstar: %f,%f\n", source->peak->xf, source->peak->yf);
+        }
+
+        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);
+        } else { 
+            // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but
+            // only slightly: pmSourceRedefinePixels uses the readout to pass the pointers to
+            // the parent image data.  I guess the API could be simplified: we could recover
+            // this from the source in the function
+
+            pmReadout tmpReadout;
+            tmpReadout.image    = (psImage *)source->pixels->parent;
+            tmpReadout.mask     = (psImage *)source->maskView->parent;
+            tmpReadout.variance = (psImage *)source->variance->parent;
+
+            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+            pmSourceRedefinePixels (source, &tmpReadout, source->peak->x, source->peak->y, BIG_RADIUS + 2);
+
+            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);
+            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
+        }
+        if (!status) {
+            source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
+            Nfail ++;
+            psErrorClear();
+            continue;
+        }
+        Nmoments ++;
+
+        // re-try big sources or not??
+        // if (source->moments->Mrf < 2.0*SIGMA)
     }
 
@@ -473 +482 @@
 
     float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
-    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
     psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
+
+    // when we set the window, we are not attempting to measure spatial variations; we can use a somewhat higher S/N limit
+    // since we are using all sources (true?)
+    float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
 
     // XXX move this to a config file?
@@ -480 +492 @@
     float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
     float Sout[NSIGMA];
-
-    // this sorts by peak->SN
-    sources = psArraySort (sources, pmSourceSortBySN);
+    int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
+
+    // this sorts by peak->rawFlux
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // loop over radii:
@@ -495 +508 @@
 
             // skip faint sources for moments measurement
-            if (source->peak->SN < MIN_SN) {
+            if (sqrt(source->peak->detValue) < MIN_SN) {
                 continue;
             }
@@ -510 +523 @@
         // determine the PSF parameters from the source moment values
         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, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
+        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]);
 
 #if 0
@@ -531 +544 @@
 
         Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
+        Nout[i] = psfClump.nStars;
     }
 
     // we are looking for sigma for which Sout = 0.65 (or some other value)
+    int Nstars = 0;
     float Sigma = NAN;
     float minS = Sout[0];
@@ -549 +564 @@
         if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
         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]);
     }
     psAssert (isfinite(Sigma), "did we miss a case?");
@@ -558 +574 @@
     psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
     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);

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM