Changeset 42435

Timestamp:

Mar 24, 2023, 9:45:49 AM (3 years ago)

Author:

eugene

Message:

use a selected subset of sources to calculate source stats; new feature detection code to remove streaks from the sources used to set the PSF stars

Location:

branches/eam_branches/ipp-20230313/psphot/src

Files:

• Makefile.am (modified) (1 diff)
• psphot.h (modified) (1 diff)
• psphotFindFeatures.c (added)
• psphotSourceStats.c (modified) (7 diffs)

branches/eam_branches/ipp-20230313/psphot/src/Makefile.am

@@ -247 +247 @@
         psphotSetNFrames.c             \
         psphotSourceMemory.c           \
+        psphotFindFeatures.c           \
         psphotChipParams.c              \
         psphotGalaxyParams.c

branches/eam_branches/ipp-20230313/psphot/src/psphot.h

@@ -592 +592 @@
 extern bool psphotINpsphotStack;
 
+bool psphotFindFeatures (psArray *sources);
+
 
 #endif

branches/eam_branches/ipp-20230313/psphot/src/psphotSourceStats.c

@@ -120 +120 @@
         source->moments = pmMomentsAlloc();
 
+        // XXX can this ever be set at this point? (we just allocated the source and moments)
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
             fprintf (stderr, "moment failure\n");
@@ -142 +143 @@
         return true;
     }
+
+    // XXX I want to identify sources which are on lines.  these can be excluded from
+    // the Moments and PSF analysis
+
+    psphotFindFeatures (sources);
 
     if (setWindow) {
@@ -539 +545 @@
     sources = psArraySort (sources, pmSourceSortByFlux);
 
+    // generate an array of the sources which we want to use for this analysis
+    // XXX move max source number to config
+    psArray *sourceSubset = psArrayAllocEmpty (100);
+    for (int i = 0; (i < sources->n) && (i < 400); i++) {
+
+        pmSource *source = sources->data[i];
+
+        // skip faint sources for moments measurement
+        if (sqrt(source->peak->detValue) < MIN_SN) continue;
+
+        // skip sources on lines
+        if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
+        psArrayAdd (sourceSubset, 100, source);
+    }
+
     // loop over radii:
     for (int i = 0; i < nsigma; i++) {
-
-        // XXX move max source number to config
-        for (int j = 0; (j < sources->n) && (j < 400); j++) {
-
-            pmSource *source = sources->data[j];
+        for (int j = 0; j < sourceSubset->n; j++) {
+            pmSource *source = sourceSubset->data[j];
             psAssert (source->moments, "force moments to exist");
             source->moments->nPixels = 0;
@@ -565 +583 @@
           sprintf (name, "moments.v%d.dat", i);
           FILE *fout = fopen (name, "w");
-          for (int j = 0; j < sources->n; j++) {
-            pmSource *source = sources->data[j];
+          for (int j = 0; j < sourceSubset->n; j++) {
+            pmSource *source = sourceSubset->data[j];
             psAssert (source->moments, "force moments to exist");
             source->moments->nPixels = 0;
@@ -581 +599 @@
 
         // 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_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
+        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sourceSubset, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
         psLogMsg ("psphot", 3, "sigma guess (pix) %.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);
+        Rmin[i] = pmSourceMinKronRadius(sourceSubset, PSF_SN_LIM);
 
 #if 0
@@ -600 +618 @@
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
         if (pmVisualTestLevel("psphot.moments.full", 2)) {
-            psphotVisualPlotMoments (recipe, analysis, sources);
+            psphotVisualPlotMoments (recipe, analysis, sourceSubset);
         }
 #endif
@@ -608 +626 @@
         Sout[i] = (Nout[i] == 0) ? NAN : sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     }
+
+    psFree (sourceSubset);
 
     // we are looking for sigma for which Sout = 0.65 (or some other value)