Changeset 25358 for branches/eam_branches/20090715

Timestamp:

Sep 13, 2009, 4:55:31 PM (17 years ago)

Author:

eugene

Message:

iteratively determine optimal Gaussian window sigma value; adjust clump parameters based on window sigma

File:

• branches/eam_branches/20090715/psphot/src/psphotSourceStats.c (modified) (2 diffs)

branches/eam_branches/20090715/psphot/src/psphotSourceStats.c

@@ -21 +21 @@
     float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     if (!status) return NULL;
+
+    OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
+
     char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
     if (!status) return NULL;
@@ -60 +63 @@
         psphotVisualShowMoments (sources);
         return sources;
+    }
+
+    // XXX *** this section is an attempt to iteratively determine the best value for sigma of the moments weighting Gaussian
+    {
+        float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
+        if (!status) return false;
+
+        float sigma[4] = {1.0, 2.0, 3.0, 4.5};
+        float Sout[4];
+
+        // XXX check that this sorts by peak->SN
+        sources = psArraySort (sources, pmSourceSortBySN);
+
+        // loop over radii:
+        for (int i = 0; i < 4; i++) {
+
+            for (int j = 0; (j < sources->n) && (j < 400); j++) {
+ 
+                pmSource *source = sources->data[j];
+                psAssert (source->moments, "force moments to exist");
+                source->moments->nPixels = 0;
+
+                // skip faint sources for moments measurement
+                if (source->peak->SN < MIN_SN) {
+                    continue;
+                }
+
+                // measure basic source moments
+                status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
+            }
+
+            // choose a grid scale that is a fixed fraction of the psf sigma^2
+            psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
+            psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
+            psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
+
+            // determine the PSF parameters from the source moment values
+            pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
+            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]);
+
+            psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
+            if (!regionMD) {
+                regionMD = psMetadataAlloc();
+                psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
+                psFree (regionMD);
+            }
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
+            psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
+            
+            psphotVisualPlotMoments (recipe, sources);
+
+            Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
+        }
+
+        // we are looking for sigma for which Sout = 0.65 (or some other value)
+
+        float Sigma = NAN;
+        float minS = Sout[0];
+        float maxS = Sout[0];
+        for (int i = 0; i < 4; i++) {
+            minS = PS_MIN(Sout[i], minS);
+            maxS = PS_MAX(Sout[i], maxS);
+        }
+        if (minS > 0.65) Sigma = sigma[3];
+        if (maxS < 0.65) Sigma = sigma[0];
+
+        for (int i = 0; i < 3; i++) {
+            if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
+            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]);
+        }
+        psAssert (isfinite(Sigma), "did we miss a case?");
+        
+        // choose a grid scale that is a fixed fraction of the psf sigma^2
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
+        psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
+
+        psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
     }