Changeset 42710

Timestamp:

Jul 5, 2024, 12:29:41 PM (2 years ago)

Author:

eugene

Message:

add option in psphot to fit the PSF star positions as well as shapes (default is false); add option to use robustMedian vs clippedMean for moments-based source sizes and PSF star selection; reduce chance of choosing galaxies for PSF stars

Location:

branches/eam_branches/ipp-20230313/psModules/src/objects

Files:

• pmPSF.h (modified) (1 diff)
• pmPSFtryFitEXT.c (modified) (1 diff)
• pmSource.c (modified) (4 diffs)
• pmSource.h (modified) (1 diff)
• pmSourceFitModel.c (modified) (2 diffs)
• pmSourceFitModel.h (modified) (1 diff)

branches/eam_branches/ipp-20230313/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;

branches/eam_branches/ipp-20230313/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;

branches/eam_branches/ipp-20230313/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");
@@ -601 +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);
@@ -749 +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 ++;
@@ -768 +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;

branches/eam_branches/ipp-20230313/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
 );
 

branches/eam_branches/ipp-20230313/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

branches/eam_branches/ipp-20230313/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,