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

Timestamp:

Oct 2, 2009, 3:12:47 PM (17 years ago)

Author:

eugene

Message:

check in changes from gene@development branch : extensive changes to moments calculation, psf model generation, aperture residuals; improvements to extended source analysis

Location:

trunk/psphot/src

Files:

• . (modified) (1 prop)
• psphotSourceStats.c (modified) (7 diffs)

trunk/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy

trunk/psphot/src/psphotSourceStats.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
+bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources);
+
+psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) {
 
     bool status = false;
@@ -21 +23 @@
     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 +65 @@
         psphotVisualShowMoments (sources);
         return sources;
+    }
+
+    if (setWindow) {
+        if (!psphotSetMomentsWindow(recipe, sources)) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
+            return NULL;
+        }
     }
 
@@ -144 +156 @@
     float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     if (!status) return false;
-    float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN");
-    if (!status) return false;
     float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
     if (!status) return false;
@@ -194 +204 @@
         }
 
-        // measure basic source moments
-        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
+        // measure basic source moments (no S/N clipping on input pixels)
+        status = pmSourceMoments (source, RADIUS, SIGMA, 0.0);
         if (status) {
             Nmoments ++;
@@ -205 +215 @@
         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, MIN_PIXEL_SN);
+        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, 0.0);
         if (status) {
             source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
@@ -231 +241 @@
 }
 
-# if (0)
-bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) {
-
-    bool status = false;
-    float BIG_RADIUS;
-
-    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    if (!status) return false;
-    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
-    if (!status) return false;
-    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
-    if (!status) return false;
-
-    // bit-masks to test for good/bad pixels
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
-
-    // bit-mask to mark pixels not used in analysis
-    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
-
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
-
-    // threaded measurement of the sources moments
-    int Nfail = 0;
-    int Nmoments = 0;
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-
-        // skip faint sources for moments measurement
-        if (source->peak->SN < MIN_SN) {
-            continue;
-        }
-
-        // measure a local sky value
-        // the local sky is now ignored; kept here for reference only
-        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
-        if (!status) {
-            psErrorClear(); // XXX re-consider the errors raised here
-            Nfail ++;
-            continue;
-        }
-
-        // measure the local sky variance (needed if noise is not sqrt(signal))
-        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
-        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
-        if (!status) {
-            Nfail ++;
-            psErrorClear();
-            continue;
-        }
-
-        // measure basic source moments
-        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
-        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, MIN_PIXEL_SN);
-        if (status) {
-            Nmoments ++;
-            continue;
-        }
-
-        Nfail ++;
-        psErrorClear();
-        continue;
-    }
-
-    // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    *nmoments = Nmoments;
-    *nfail = Nfail;
-
+// this function attempts to iteratively determine the best value for sigma of the moments weighting Gaussian
+bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources) {
+
+    bool status;
+
+    float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
+    if (!status) return false;
+
+    // XXX move this to a config file?
+    float sigma[4] = {1.0, 2.0, 3.0, 4.5};
+    float Sout[4];
+
+    // this sorts by peak->SN
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // loop over radii:
+    for (int i = 0; i < 4; i++) {
+
+        // XXX move max source number to config
+        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 (no S/N clipping on input pixels)
+            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);
     return true;
 }
-# endif

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy