Changeset 42843 for trunk/psModules

Timestamp:

May 9, 2025, 11:28:49 AM (12 months ago)

Author:

eugene

Message:

optionally allow PSF star positions to be fitted during PSF model construction (default false); option to allow PSF clump analysis to use clipped stats (default) or robust stats; use select circular region of 2nd moment space, not square region; for bright (S/N > 25) externally-supplied sources, use the centroid (Mx,My) value instead of the peak

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/pmPSF.h (modified) (1 diff)
• src/objects/pmPSFtryFitEXT.c (modified) (1 diff)
• src/objects/pmSource.c (modified) (5 diffs)
• src/objects/pmSource.h (modified) (1 diff)
• src/objects/pmSourceFitModel.c (modified) (2 diffs)
• src/objects/pmSourceFitModel.h (modified) (1 diff)
• src/objects/pmSourceMasks.h (modified) (1 diff)
• src/objects/pmSourceMoments.c (modified) (5 diffs)

trunk/psModules/src/objects/pmPSF.h

@@ -79 +79 @@
     bool          poissonErrorsParams; ///< use poission errors for model parameter fitting
 
+    bool          fitPSFstarCoords; // Allow the PSF star positions to fit in the full fit? (pmPSFtryFitEXT)
     bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
     pmSourceFitOptions *fitOptions;

trunk/psModules/src/objects/pmPSFtryFitEXT.c

@@ -73 +73 @@
 
     // in this segment, we are fitting the full PSF model class (shape unconstrained)
-    options->fitOptions->mode = PM_SOURCE_FIT_EXT;
-
+    if (options->fitPSFstarCoords) {
+      options->fitOptions->mode = PM_SOURCE_FIT_FULL;
+    } else {
+      options->fitOptions->mode = PM_SOURCE_FIT_EXT;
+    }
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;

trunk/psModules/src/objects/pmSource.c

@@ -407 +407 @@
 *****************************************************************************/
 
-pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 SX_MIN, psF32 SY_MIN, psF32 AR_MAX)
+pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 SX_MIN, psF32 SY_MIN, psF32 AR_MAX, bool useClippedMean)
 {
     psTrace("psModules.objects", 10, "---- begin ----\n");
@@ -455 +455 @@
                 continue;
             }
+
+            // skip sources associated with possible detector features 
+            if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
 
             float Mxx = source->moments->Mxx, Myy = source->moments->Myy; // Second moments
@@ -598 +601 @@
         }
 
-        // measures stats of Sx, Sy
-        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-
-        if (!psVectorStats (stats, tmpSx, NULL, NULL, 0)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats");
-            return (emptyClump);
-        }
-        psfClump.X  = stats->clippedMean;
-        psfClump.dX = hypot(stats->clippedStdev, PSF_CLUMP_GRID_SCALE);
-
-        if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats");
-            return (emptyClump);
-        }
-        psfClump.Y  = stats->clippedMean;
-        psfClump.dY = hypot(stats->clippedStdev, PSF_CLUMP_GRID_SCALE);
+        // measures stats of Sx, Sy
+        if (useClippedMean) {
+          stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        } else {
+          stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+        }
+
+        if (!psVectorStats (stats, tmpSx, NULL, NULL, 0)) {
+          psError(PS_ERR_UNKNOWN, false, "failed to measure Sx stats");
+          return (emptyClump);
+        }
+        psfClump.X  = psStatsGetValue (stats, psStatsMeanOption (stats->options)); 
+        psfClump.dX = psStatsGetValue (stats, psStatsStdevOption (stats->options));
+        psfClump.dX = hypot(psfClump.dX, PSF_CLUMP_GRID_SCALE);
+
+        if (!psVectorStats (stats, tmpSy, NULL, NULL, 0)) {
+          psError(PS_ERR_UNKNOWN, false, "failed to measure Sy stats");
+          return (emptyClump);
+        }
+        psfClump.Y  = psStatsGetValue (stats, psStatsMeanOption (stats->options)); 
+        psfClump.dY = psStatsGetValue (stats, psStatsStdevOption (stats->options));
+        psfClump.dY = hypot(psfClump.dX, PSF_CLUMP_GRID_SCALE);              
 
         psTrace ("psModules.objects", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
@@ -746 +755 @@
             }
 
+            psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
+
+            // XXX this is crude cut to avoid using extended sources in the PSF model
+            // psphot will update this analysis based on the PSF - Kron mags
             // likely unsaturated extended source (too large to be stellar)
-            if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) {
+            // XXX old test: if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) {
+
+            if (radius > 2.0) {
                 source->type = PM_SOURCE_TYPE_EXTENDED;
                 Next ++;
                 continue;
             }
+
+            // skip sources associated with possible detector features 
+            if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
 
             // the rest are probable stellar objects
@@ -762 +780 @@
 
             // PSF star (within 1.5 sigma of clump center, S/N > limit)
-            psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
             if ((source->moments->SN > PSF_SN_LIM) && (radius < PSF_CLUMP_NSIGMA)) {
                 source->type = PM_SOURCE_TYPE_STAR;

trunk/psModules/src/objects/pmSource.h

@@ -245 +245 @@
     psF32 SX_MIN, 
     psF32 SY_MIN, 
-    psF32 AR_MAX
+    psF32 AR_MAX,
+    bool useClippedMean                 ///< use clipped mean & stdev to define clump, else robust
 );
 

trunk/psModules/src/objects/pmSourceFitModel.c

@@ -189 +189 @@
         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
+        break;
+      case PM_SOURCE_FIT_FULL:
+        // PSFstar Full model fits all shape params and Xo, Yo, Io (not sky)
+        nParams = params->n - 1;
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
         break;
       case PM_SOURCE_FIT_EXT:
@@ -270 +276 @@
     model->nPar = nParams;
 
-    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
+    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f, status: %d", myMin->iter, myMin->value, fitStatus);
 
     // save the resulting chisq, nDOF, nIter

trunk/psModules/src/objects/pmSourceFitModel.h

@@ -18 +18 @@
     PM_SOURCE_FIT_PSF,
     PM_SOURCE_FIT_EXT,
+    PM_SOURCE_FIT_FULL,
     PM_SOURCE_FIT_PSF_AND_SKY,
     PM_SOURCE_FIT_EXT_AND_SKY,

trunk/psModules/src/objects/pmSourceMasks.h

@@ -47 +47 @@
     PM_SOURCE_MODE2_MATCHED               = 0x00000004, ///< source generated based on another image
     PM_SOURCE_MODE2_ON_SPIKE              = 0x00000008, ///< > 25% of (PSF-weighted) pixels land on diffraction spike
+
     PM_SOURCE_MODE2_ON_STARCORE           = 0x00000010, ///< > 25% of (PSF-weighted) pixels land on starcore
     PM_SOURCE_MODE2_ON_BURNTOOL           = 0x00000020, ///< > 25% of (PSF-weighted) pixels land on burntool
     PM_SOURCE_MODE2_ON_CONVPOOR           = 0x00000040, ///< > 25% of (PSF-weighted) pixels land on convpoor
     PM_SOURCE_MODE2_PASS1_SRC             = 0x00000080, ///< source detected in first pass analysis
+
     PM_SOURCE_MODE2_HAS_BRIGHTER_NEIGHBOR = 0x00000100, ///< peak is not the brightest in its footprint
     PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_1     = 0x00000200, ///< flux_n / (r^2 flux_p) > 1
     PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_10    = 0x00000400, ///< flux_n / (r^2 flux_p) > 10
     PM_SOURCE_MODE2_DIFF_SELF_MATCH       = 0x00000800, ///< positive detection match is probably this source 
+
     PM_SOURCE_MODE2_SATSTAR_PROFILE       = 0x00001000, ///< saturated source is modeled with a radial profile
-
     PM_SOURCE_MODE2_ECONTOUR_FEW_PTS      = 0x00002000, ///< too few points to measure the elliptical contour
     PM_SOURCE_MODE2_RADBIN_NAN_CENTER     = 0x00004000, ///< radial bins failed with too many NaN center bin
     PM_SOURCE_MODE2_PETRO_NAN_CENTER      = 0x00008000, ///< petrosian radial bins failed with too many NaN center bin
+
     PM_SOURCE_MODE2_PETRO_NO_PROFILE      = 0x00010000, ///< petrosian not build because radial bins missing
-
     PM_SOURCE_MODE2_PETRO_INSIG_RATIO     = 0x00020000, ///< insignificant measurement of petrosian ratio
     PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN   = 0x00040000, ///< petrosian ratio in the 0th bin (likely bad)
-    
     PM_SOURCE_MODE2_EXT_FITS_RUN          = 0x00080000, ///< we attempted to run extended fits on this source
+
     PM_SOURCE_MODE2_EXT_FITS_FAIL         = 0x00100000, ///< at least one of the model fits failed
     PM_SOURCE_MODE2_EXT_FITS_RETRY        = 0x00200000, ///< one of the model fits was re-tried with new window
     PM_SOURCE_MODE2_EXT_FITS_NONE         = 0x00400000, ///< ALL of the model fits failed
-    
     PM_SOURCE_MODE2_ON_GHOST              = 0x00800000, ///< > 25% of (PSF-weighted) pixels land on ghost
+
     PM_SOURCE_MODE2_ON_CROSSTALK          = 0x01000000, ///< peaks land on electronic crostalk
     PM_SOURCE_MODE2_ON_CTE                = 0x02000000, ///< peaks land on CTE region
+    PM_SOURCE_MODE2_ON_LINE               = 0x04000000, ///< peaks land on Linear feature
 
     

trunk/psModules/src/objects/pmSourceMoments.c

@@ -231 +231 @@
     }
 
-    // we only update the centroid if the position is not supplied from elsewhere
+    // For faint externally-supplied sources, the centroids are unreliable.
+    // In those cases, we should keep the peak position instead of the centroid position.
     bool skipCentroid = false;
     skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
@@ -245 +246 @@
     }
 
-    // only skip negative sum sources if the sources are not forced
+    // Skip negative-sum sources, but only if the sources are not forced (externally supplied)
     if (!skipCentroid && (Sum <= 0)) {
         psTrace ("psModules.objects", 3, "insufficient significant pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
@@ -259 +260 @@
     psTrace ("psModules.objects", 5, "id: %d, sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", source->id, sky, Mx, My, Sum, X1, Y1, numPixels);
 
+    // XXX EAM : I need to move this block so I can use SN to choose if we keep the centroids or not
+    source->moments->Sum = Sum;
+    source->moments->SN  = Sum / sqrt(Var);
+    source->moments->Peak = peakPixel;
+    source->moments->nPixels = numPixels;
+
     // add back offset of peak in primary image
     // also offset from pixel index to pixel coordinate
@@ -264 +271 @@
     // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
 
-    if (skipCentroid) {
+    // XXX for a test, do NOT skip the centroid for relatively strong sources, even if forced
+    // XXX this is a somewhat arbitrary limit.  use PSF_SN_LIM (see below) instead?
+    if (skipCentroid && (source->moments->SN < 25)) {
         source->moments->Mx = source->peak->xf;
         source->moments->My = source->peak->yf;
@@ -270 +279 @@
         if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
             psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
+
+            // XXX should we really be rejecting these sources, or just keeping the peak position?
             return (false);
+
+            source->moments->Mx = source->peak->xf;
+            source->moments->My = source->peak->yf;
+            return (true);
         }
         source->moments->Mx = Mx + xGuess;
         source->moments->My = My + yGuess;
     }
-
-    source->moments->Sum = Sum;
-    source->moments->SN  = Sum / sqrt(Var);
-    source->moments->Peak = peakPixel;
-    source->moments->nPixels = numPixels;
 
     return true;