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

Timestamp:

May 9, 2025, 11:21:40 AM (14 months ago)

Author:

eugene

Message:

add psphotFindFeatures; enable forced-photometry between cameras; do NOT skip negative-flux sources in full-force photometry; make PSF residual component generation more robust: IRLS vs median; exclude pixels with poor values (scaled value > 1.5, relative to normalized PSF); exclude pixels with too few inputs (<5); allow PSF clump to use clipped stats or robust stats

Location:

trunk/psphot

Files:

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

trunk/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");
@@ -144 +145 @@
 
     if (setWindow) {
+        // identify sources on lines, drop from MomentsWindow analysis
+        psphotFindFeatures (readout, sources);
+
         if (!psphotSetMomentsWindow(recipe, readout->analysis, sources, maskVal)) {
             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
@@ -503 +507 @@
     }
 
+    // use clipped stats or robust stats to define the clump
+    bool PSF_CLUMP_USE_CLIPPED_STATS = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_USE_CLIPPED_STATS");
+    if (!status) { PSF_CLUMP_USE_CLIPPED_STATS = true; }
+
     // 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?)
@@ -539 +547 @@
     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 +585 @@
           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 +601 @@
 
         // 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, PSF_CLUMP_USE_CLIPPED_STATS);
         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 +620 @@
         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 +628 @@
         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)