From 87e3bf13a7c5a18f192c8ce613e3dc7d03486f7b Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Mon, 2 Mar 2026 20:09:21 -0800 Subject: [PATCH 1/8] Add single-stamp fallback when coadd detection fails When source detection on the coadd image fails for a guider, fall back to trying detection on individual stamps evenly spread across the sequence. Stamps with low peak pixel SNR are skipped as truly empty. Returns the detection with the highest SNR. This recovers stars that are lost in the coadd due to drift or artifacts. Adds nFallbackStamps and peakSnrThreshold to GuiderStarTrackerConfig. (DM-54263) --- python/lsst/summit/utils/guiders/detection.py | 74 +++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/python/lsst/summit/utils/guiders/detection.py b/python/lsst/summit/utils/guiders/detection.py index a2c95bdc..8578c182 100644 --- a/python/lsst/summit/utils/guiders/detection.py +++ b/python/lsst/summit/utils/guiders/detection.py @@ -87,6 +87,11 @@ class GuiderStarTrackerConfig: Aperture size in arcseconds for star detection. gain : `float` Gain factor for the guider data, used in flux calculations. + nFallbackStamps : `int` + Number of individual stamps to sample when coadd detection fails. + peakSnrThreshold : `float` + Minimum peak pixel SNR to consider a stamp non-empty during + single-stamp fallback detection. """ minSnr: float = 10.0 @@ -96,6 +101,8 @@ class GuiderStarTrackerConfig: cutOutSize: int = 50 aperSizeArcsec: float = 5.0 gain: float = 1.0 + nFallbackStamps: int = 10 + peakSnrThreshold: float = 5.0 def trackStarAcrossStamp( @@ -597,6 +604,69 @@ def makeEllipticalGaussianStar( return image +def _detectOnSingleStamps( + guiderData: GuiderData, + guiderName: str, + config: GuiderStarTrackerConfig, + apertureRadius: int, + log: logging.Logger, +) -> pd.DataFrame: + """Try source detection on individual stamps when coadd detection fails. + + Samples up to ``config.nFallbackStamps`` stamps evenly spread across + the sequence. Returns the detection with the highest SNR, or an empty + DataFrame if no star is found on any stamp. + + A stamp is considered truly empty if its peak pixel SNR + (max - median) / std is below ``config.peakSnrThreshold``. + """ + nStamps = len(guiderData) + if nStamps == 0: + return pd.DataFrame(columns=DEFAULT_COLUMNS) + + nSample = min(config.nFallbackStamps, nStamps) + indices = np.linspace(0, nStamps - 1, nSample, dtype=int) + + bestSources = pd.DataFrame(columns=DEFAULT_COLUMNS) + bestSnr = 0.0 + nTrulyEmpty = 0 + + for idx in indices: + stamp = guiderData[guiderName, idx].astype(np.float32) + arr = stamp - np.nanmin(stamp) + + # Quick peak-SNR check before running full detection + med = np.nanmedian(arr) + std = np.nanstd(arr) + if std > 0: + peakSnr = (np.nanmax(arr) - med) / std + else: + peakSnr = 0.0 + + if peakSnr < config.peakSnrThreshold: + nTrulyEmpty += 1 + continue + + sources = runSourceDetection( + arr, + threshold=config.minSnr, + apertureRadius=apertureRadius, + cutOutSize=config.cutOutSize, + gain=config.gain, + ) + if not sources.empty and sources["snr"].max() > bestSnr: + bestSnr = sources["snr"].max() + bestSources = sources + + if not bestSources.empty: + log.info( + f"Single-stamp fallback recovered source on {guiderName} " + f"(SNR={bestSnr:.1f}, {nTrulyEmpty}/{nSample} stamps empty)" + ) + + return bestSources + + def buildReferenceCatalog( guiderData: GuiderData, log: logging.Logger, @@ -638,6 +708,10 @@ def buildReferenceCatalog( cutOutSize=cutOutSize, gain=gain, ) + if sources.empty: + # Fallback: try detection on individual stamps. The coadd can + # wash out stars that drift between frames or have artifacts. + sources = _detectOnSingleStamps(guiderData, guiderName, config, apertureRadius, log) if sources.empty: log.warning(f"No sources detected in `buildReferenceCatalog`for {guiderName} in {expId}. ") continue From 641d4e232e74f61efb518eea1bfc33aa647f669a Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Tue, 3 Mar 2026 05:10:12 -0800 Subject: [PATCH 2/8] Improve source detection robustness in guider pipeline - Fix NaN handling on cutout window in measureStarOnStamp: replace NaN from border-clipped Cutout2D with median instead of rejecting the whole cutout, recovering stars near stamp edges. - Add GalSim failure warning in runSourceDetection when FootprintSet finds sources but GalSim fails on all of them. - Fix fluxErr sqrt warning in runGalSim: clamp negative flux/gain to 0. - Track detection_method column ("coadd" or "single_stamp") in buildReferenceCatalog. - Replace absolute flux threshold in isBlankImage with peak pixel SNR. - Add module-level logger for detection diagnostics. (DM-54305) --- python/lsst/summit/utils/guiders/detection.py | 45 ++++++++++++++----- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/python/lsst/summit/utils/guiders/detection.py b/python/lsst/summit/utils/guiders/detection.py index 8578c182..456f1638 100644 --- a/python/lsst/summit/utils/guiders/detection.py +++ b/python/lsst/summit/utils/guiders/detection.py @@ -43,6 +43,8 @@ from .reading import GuiderData +log = logging.getLogger(__name__) + _DEFAULT_COLUMNS: str = ( "trackid detector expid elapsed_time dalt daz dtheta dx dy " "fwhm xroi yroi xccd yccd xroi_ref yroi_ref xccd_ref yccd_ref " @@ -337,6 +339,7 @@ def runSourceDetection( if not footprints: return pd.DataFrame(columns=DEFAULT_COLUMNS) + nFootprints = len(footprints.getFootprints()) results = [] for fp in footprints.getFootprints(): # Create a cutout of the image around the footprint @@ -345,6 +348,11 @@ def runSourceDetection( if not star.empty: results.append(star) if not results: + if nFootprints > 0: + log.warning( + f"FootprintSet found {nFootprints} sources but GalSim " + f"failed on all of them (cutOutSize={cutOutSize})." + ) return pd.DataFrame(columns=DEFAULT_COLUMNS) df = pd.concat([sf for sf in results], ignore_index=True) return df @@ -380,11 +388,17 @@ def measureStarOnStamp( StarMeasurement object with populated fields (may be empty on failure). """ cutout = getCutouts(stamp, refCenter, cutoutSize=cutOutSize) - data = cutout.data + data = cutout.data.copy() - if np.all(data == 0) | (not np.isfinite(data).all()): + if np.all(data == 0): return StarMeasurement() + # Replace NaN (from border-clipped cutouts) with the median so + # GalSim can still measure stars near the stamp edge. + nan_mask = ~np.isfinite(data) + if nan_mask.any(): + data[nan_mask] = np.nanmedian(data) + # 1) Subtract the background annulus = (apertureRadius * 1.0, apertureRadius * 2) dataBkgSub, bkgStd = annulusBackgroundSubtraction(data, annulus) @@ -446,7 +460,7 @@ def runGalSim( ellipticity = np.sqrt(e1**2 + e2**2) nEff = 2 * np.pi * sigma**2 * np.sqrt(1 - ellipticity**2) shotNoise = np.sqrt(nEff * bkgStd**2) - fluxErr = np.sqrt(flux / gain + shotNoise**2) + fluxErr = np.sqrt(max(0, flux / gain) + shotNoise**2) snr = flux / (shotNoise + 1e-9) if shotNoise > 0 else 0.0 # Calculate second moments @@ -708,14 +722,18 @@ def buildReferenceCatalog( cutOutSize=cutOutSize, gain=gain, ) + detectionMethod = "coadd" if sources.empty: # Fallback: try detection on individual stamps. The coadd can # wash out stars that drift between frames or have artifacts. sources = _detectOnSingleStamps(guiderData, guiderName, config, apertureRadius, log) + detectionMethod = "single_stamp" if sources.empty: log.warning(f"No sources detected in `buildReferenceCatalog`for {guiderName} in {expId}. ") continue + sources["detection_method"] = detectionMethod + sources.sort_values(by=["snr"], ascending=False, inplace=True) sources.reset_index(drop=True, inplace=True) @@ -755,23 +773,28 @@ def getCutouts(imageArray: np.ndarray, refCenter: tuple[float, float], cutoutSiz return Cutout2D(imageArray, (refX, refY), size=cutoutSize, mode="partial", fill_value=np.nan) -def isBlankImage(image: np.ndarray, fluxMin: int = 300) -> bool: +def isBlankImage(image: np.ndarray, peakSnrMin: float = 5.0) -> bool: """ Returns True if the image has no significant source (e.g., no star). + Uses peak pixel SNR: (max - median) / std. This is more robust than + an absolute flux threshold because it adapts to the noise level. + Parameters ---------- - image : 2D array - Image data (float or int). - fluxMin : float - Minimum deviation from the background median to be considered a source. + image : `np.ndarray` + 2D image data. + peakSnrMin : `float` + Minimum peak pixel SNR to consider the image non-blank. Returns ------- bool True if the image is blank, False otherwise. - (no pixel deviates more than flux_min) """ med = np.nanmedian(image) - diff = np.abs(image - med) - return not np.any(diff > fluxMin) + std = np.nanstd(image) + if std <= 0: + return True + peakSnr = (np.nanmax(image) - med) / std + return peakSnr < peakSnrMin From b5372f93739164ecff672e3f1bb1a0a171526f88 Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Tue, 3 Mar 2026 06:37:35 -0800 Subject: [PATCH 3/8] Use MAD-based robust std in isBlankImage Replace np.nanstd with MAD (median absolute deviation) scaled by 1.4826 for a robust noise estimate. The previous np.nanstd was inflated by streaks and negative artifacts, causing isBlankImage to misclassify guiders with real stars as blank. (DM-54305) --- python/lsst/summit/utils/guiders/detection.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/lsst/summit/utils/guiders/detection.py b/python/lsst/summit/utils/guiders/detection.py index 456f1638..9c7ccfeb 100644 --- a/python/lsst/summit/utils/guiders/detection.py +++ b/python/lsst/summit/utils/guiders/detection.py @@ -793,7 +793,8 @@ def isBlankImage(image: np.ndarray, peakSnrMin: float = 5.0) -> bool: True if the image is blank, False otherwise. """ med = np.nanmedian(image) - std = np.nanstd(image) + mad = np.nanmedian(np.abs(image - med)) + std = 1.4826 * mad # MAD-based robust std estimate if std <= 0: return True peakSnr = (np.nanmax(image) - med) / std From 9f3f0a2f31d67ee25d8aad55fbdcde30f7563aff Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Tue, 3 Mar 2026 06:51:10 -0800 Subject: [PATCH 4/8] Add flux threshold to isBlankImage alongside MAD-based SNR check Restores the original fluxMin=300 ADU absolute threshold as the primary check so bright stars are never misclassified as blank, and keeps the MAD-based peakSNR as a secondary check for fainter stars in low-noise images. --- python/lsst/summit/utils/guiders/detection.py | 22 ++++++++++++++----- 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/python/lsst/summit/utils/guiders/detection.py b/python/lsst/summit/utils/guiders/detection.py index 9c7ccfeb..42b67ac8 100644 --- a/python/lsst/summit/utils/guiders/detection.py +++ b/python/lsst/summit/utils/guiders/detection.py @@ -773,19 +773,22 @@ def getCutouts(imageArray: np.ndarray, refCenter: tuple[float, float], cutoutSiz return Cutout2D(imageArray, (refX, refY), size=cutoutSize, mode="partial", fill_value=np.nan) -def isBlankImage(image: np.ndarray, peakSnrMin: float = 5.0) -> bool: +def isBlankImage(image: np.ndarray, fluxMin: float = 300, peakSnrMin: float = 5.0) -> bool: """ Returns True if the image has no significant source (e.g., no star). - Uses peak pixel SNR: (max - median) / std. This is more robust than - an absolute flux threshold because it adapts to the noise level. + An image is considered non-blank if the peak flux above the median + exceeds ``fluxMin`` OR the peak pixel SNR (using MAD-based robust + std) exceeds ``peakSnrMin``. Parameters ---------- image : `np.ndarray` 2D image data. + fluxMin : `float` + Minimum peak flux above median (ADU) to consider non-blank. peakSnrMin : `float` - Minimum peak pixel SNR to consider the image non-blank. + Minimum peak pixel SNR to consider non-blank. Returns ------- @@ -793,9 +796,16 @@ def isBlankImage(image: np.ndarray, peakSnrMin: float = 5.0) -> bool: True if the image is blank, False otherwise. """ med = np.nanmedian(image) + peakFlux = np.nanmax(image) - med + + # Absolute flux check + if peakFlux > fluxMin: + return False + + # SNR check with MAD-based robust std mad = np.nanmedian(np.abs(image - med)) - std = 1.4826 * mad # MAD-based robust std estimate + std = 1.4826 * mad if std <= 0: return True - peakSnr = (np.nanmax(image) - med) / std + peakSnr = peakFlux / std return peakSnr < peakSnrMin From 931465170e4ec04d14ea05fe22c7d2323023a4da Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Tue, 3 Mar 2026 09:46:08 -0800 Subject: [PATCH 5/8] Fall back to flux-weighted centroid when GalSim fails When FindAdaptiveMom cannot converge, use the 0th moment (flux-weighted center of mass) and aperture photometry instead of discarding the star. Adds galsim_failed flag to StarMeasurement and DEFAULT_COLUMNS for downstream tracking. --- python/lsst/summit/utils/guiders/detection.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/python/lsst/summit/utils/guiders/detection.py b/python/lsst/summit/utils/guiders/detection.py index 42b67ac8..06e7346f 100644 --- a/python/lsst/summit/utils/guiders/detection.py +++ b/python/lsst/summit/utils/guiders/detection.py @@ -51,7 +51,7 @@ "dxfp dyfp xfp yfp alt az xfp_ref yfp_ref alt_ref az_ref " "xerr yerr theta theta_err theta_ref flux flux_err magoffset snr " "ixx iyy ixy e1 e2 e1_altaz e2_altaz " - "ampname timestamp stamp detid filter exptime " + "ampname timestamp stamp detid filter exptime galsim_failed " ) DEFAULT_COLUMNS: tuple[str, ...] = tuple(_DEFAULT_COLUMNS.split()) @@ -220,6 +220,7 @@ class StarMeasurement: flux: float = field(default=np.nan) flux_err: float = field(default=0.0) snr: float = field(default=0.0) + galsim_failed: bool = field(default=False) def toDataFrame(self) -> pd.DataFrame: """ @@ -406,6 +407,17 @@ def measureStarOnStamp( # 2) Track the star across all stamps for this guider star = runGalSim(dataBkgSub, gain=gain, bkgStd=bkgStd) + # If GalSim fails, fall back to flux-weighted centroid + aperture photometry. + # This gives valid x, y, flux for tracking even without shape info. + if not np.isfinite(star.xroi): + weights = np.clip(dataBkgSub, 0, None) + total = np.nansum(weights) + if total > 0: + yy, xx = np.indices(dataBkgSub.shape, dtype=float) + star.xroi = float(np.nansum(xx * weights) / total) + star.yroi = float(np.nansum(yy * weights) / total) + star.galsim_failed = True + # 3) Make aperture photometry measurements # Galsim flux is the normalization of the Gaussian, not w/ fixed aper. star.runAperturePhotometry(dataBkgSub, apertureRadius, gain=gain, bkgStd=bkgStd) From f25f7faa416c3d2a04b828e42386298a2533611e Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Wed, 4 Mar 2026 04:55:11 -0800 Subject: [PATCH 6/8] Revert flux-weighted centroid fallback The fallback in measureStarOnStamp corrupted buildReferenceCatalog by returning non-empty results for GalSim failures, which prevented _detectOnSingleStamps from running. Net result was a regression (208 missed vs 191 without fallback). Reverted entirely. --- python/lsst/summit/utils/guiders/detection.py | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/python/lsst/summit/utils/guiders/detection.py b/python/lsst/summit/utils/guiders/detection.py index 06e7346f..42b67ac8 100644 --- a/python/lsst/summit/utils/guiders/detection.py +++ b/python/lsst/summit/utils/guiders/detection.py @@ -51,7 +51,7 @@ "dxfp dyfp xfp yfp alt az xfp_ref yfp_ref alt_ref az_ref " "xerr yerr theta theta_err theta_ref flux flux_err magoffset snr " "ixx iyy ixy e1 e2 e1_altaz e2_altaz " - "ampname timestamp stamp detid filter exptime galsim_failed " + "ampname timestamp stamp detid filter exptime " ) DEFAULT_COLUMNS: tuple[str, ...] = tuple(_DEFAULT_COLUMNS.split()) @@ -220,7 +220,6 @@ class StarMeasurement: flux: float = field(default=np.nan) flux_err: float = field(default=0.0) snr: float = field(default=0.0) - galsim_failed: bool = field(default=False) def toDataFrame(self) -> pd.DataFrame: """ @@ -407,17 +406,6 @@ def measureStarOnStamp( # 2) Track the star across all stamps for this guider star = runGalSim(dataBkgSub, gain=gain, bkgStd=bkgStd) - # If GalSim fails, fall back to flux-weighted centroid + aperture photometry. - # This gives valid x, y, flux for tracking even without shape info. - if not np.isfinite(star.xroi): - weights = np.clip(dataBkgSub, 0, None) - total = np.nansum(weights) - if total > 0: - yy, xx = np.indices(dataBkgSub.shape, dtype=float) - star.xroi = float(np.nansum(xx * weights) / total) - star.yroi = float(np.nansum(yy * weights) / total) - star.galsim_failed = True - # 3) Make aperture photometry measurements # Galsim flux is the normalization of the Gaussian, not w/ fixed aper. star.runAperturePhotometry(dataBkgSub, apertureRadius, gain=gain, bkgStd=bkgStd) From dcea84d69871012934dda26faf722b1c261a3758 Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Wed, 4 Mar 2026 05:03:35 -0800 Subject: [PATCH 7/8] Add per-candidate rejection diagnostics when quality cuts reject all stars When all reference catalog candidates for a guider are rejected by quality cuts, log the best candidate (highest SNR) with flux, snr, ellipticity, fwhm, and position, plus which cut rejected it. Distinguishes between "no detection" and "detected but rejected". Uses the same mask logic as applyQualityCuts for accuracy. --- python/lsst/summit/utils/guiders/tracking.py | 73 +++++++++++++++++++- 1 file changed, 72 insertions(+), 1 deletion(-) diff --git a/python/lsst/summit/utils/guiders/tracking.py b/python/lsst/summit/utils/guiders/tracking.py index 9b64ed18..738d111d 100644 --- a/python/lsst/summit/utils/guiders/tracking.py +++ b/python/lsst/summit/utils/guiders/tracking.py @@ -95,6 +95,7 @@ def trackGuiderStars(self, refCatalog: None | pd.DataFrame = None) -> pd.DataFra DataFrame with tracked stars and their properties, including positions, fluxes, and residual offsets. """ + _refCatalog = None if refCatalog is None: self.log.info("Using self-generated refcat") _refCatalog = buildReferenceCatalog(self.guiderData, self.log, self.config) @@ -108,7 +109,23 @@ def trackGuiderStars(self, refCatalog: None | pd.DataFrame = None) -> pd.DataFra for guiderName in self.guiderData.guiderNames: refDet = refCatalog.query(f"detector == '{guiderName}'") if refDet.empty: - self.log.warning(f"No stars found in reference catalog for {guiderName} in {self.expid}.") + # Diagnose: check what was in the unfiltered refcat + preCut = _refCatalog.query(f"detector == '{guiderName}'") if _refCatalog is not None else None + if preCut is not None and not preCut.empty: + best = preCut.sort_values("snr", ascending=False).iloc[0] + reason = _diagnoseQualityCutRejections( + preCut.sort_values("snr", ascending=False).iloc[:1], self.shape, self.config + )[0] + self.log.warning( + f"No stars in refcat for {guiderName} in {self.expid}. " + f"Best candidate: flux={best['flux']:.0f} snr={best['snr']:.1f} " + f"e={np.hypot(best['e1'], best['e2']):.2f} fwhm={best['fwhm']:.1f} " + f"pos=({best['xroi']:.0f},{best['yroi']:.0f}). Rejected: {reason}" + ) + else: + self.log.warning( + f"No stars detected for {guiderName} in {self.expid} " f"(empty before quality cuts)." + ) continue table = self._trackStarForOneGuider(refDet, guiderName) @@ -336,6 +353,60 @@ def applyQualityCuts( return stars[mask].copy() +def _diagnoseQualityCutRejections( + stars: pd.DataFrame, shape: tuple[float, float], config: GuiderStarTrackerConfig +) -> list[str]: + """Return per-row rejection reasons using the same masks as + applyQualityCuts. + + Parameters + ---------- + stars : `pd.DataFrame` + Pre-cut DataFrame (rows that were all rejected). + shape : `tuple[float, float]` + ROI shape (rows, cols). + config : `GuiderStarTrackerConfig` + Configuration with quality cut thresholds. + + Returns + ------- + reasons : `list[str]` + One string per row describing which cuts failed. + """ + minSnr = config.minSnr + maxEllipticity = config.maxEllipticity + edgeMargin = config.edgeMargin + roiCols, roiRows = shape + + # Same masks as applyQualityCuts + snrOk = (stars["snr"] >= minSnr) & (stars["flux"] > 0) & (stars["flux_err"] > 0) + + eabs = np.hypot(stars["e1"], stars["e2"]) + ellipOk = (stars["e1"].abs() <= maxEllipticity) & (stars["e1"].abs() <= maxEllipticity) + ellipOk &= eabs <= maxEllipticity + + edgeOk = ( + (stars["xroi"] >= edgeMargin) + & (stars["xroi"] <= roiRows - edgeMargin) + & (stars["yroi"] >= edgeMargin) + & (stars["yroi"] <= roiCols - edgeMargin) + ) + + reasons = [] + for i, (_, row) in enumerate(stars.iterrows()): + r = [] + if not snrOk.iloc[i]: + snr = row["snr"] + r.append(f"snr={snr:.1f}" if np.isfinite(snr) else "snr=NaN") + if not ellipOk.iloc[i]: + e = np.hypot(row["e1"], row["e2"]) + r.append(f"e={e:.2f}" if np.isfinite(e) else "e=NaN") + if not edgeOk.iloc[i]: + r.append(f"edge(x={row['xroi']:.0f},y={row['yroi']:.0f})") + reasons.append("|".join(r) if r else "?") + return reasons + + def setUniqueId(guiderData: GuiderData, stars: pd.DataFrame) -> pd.DataFrame: """ Assign unique IDs to tracked stars. From 3bed6003b3104665c51a9a4a756d0f4576404b0e Mon Sep 17 00:00:00 2001 From: Johnny Esteves Date: Wed, 4 Mar 2026 09:05:36 -0800 Subject: [PATCH 8/8] Use percentile instead of median for column bias subtraction As suggested by Robert, use a low percentile (default 10th) instead of the median for column bias estimation. Stars and trails are very unlikely to live in the 10th percentile, which avoids the dark dip artifacts created by the median being pulled by bright sources. The percentile is slightly noisier but eliminates the dark bands. Added biasPercentile config option to GuiderReader.get() and processStamps() to allow tuning the percentile value. --- python/lsst/summit/utils/guiders/reading.py | 33 ++++++++++++++------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/python/lsst/summit/utils/guiders/reading.py b/python/lsst/summit/utils/guiders/reading.py index 2ca65206..c1799f9c 100644 --- a/python/lsst/summit/utils/guiders/reading.py +++ b/python/lsst/summit/utils/guiders/reading.py @@ -701,6 +701,7 @@ def get( seqNum: int, doSubtractMedian: bool = True, columnMaskK: float = 50.0, + biasPercentile: float = 10.0, scienceDetNum: int = 94, ) -> GuiderData: """ @@ -713,9 +714,12 @@ def get( seqNum : `int` Sequence number. doSubtractMedian : `bool`, optional - If True, subtract median row bias from each stamp. + If True, subtract column bias from each stamp. columnMaskK : `float`, optional Threshold factor for column mask detection. + biasPercentile : `float`, optional + Percentile (0-100) for column bias estimation. + Lower values avoid star/trail contamination. scienceDetNum : `int`, optional Science detector number for WCS reference. @@ -748,6 +752,7 @@ def get( rawStampsDict, doSubtractMedian, columnMaskK, + biasPercentile, ) guiderData = GuiderData( seqNum=seqNum, @@ -797,9 +802,10 @@ def processStamps( rawStampsDict: dict[str, Stamps], doSubtractMedian: bool, columnMaskK: float, + biasPercentile: float = 10.0, ) -> dict[str, Stamps]: """ - Apply median row bias subtraction and per-stamp column masking. + Apply column bias subtraction and per-stamp column masking. Parameters @@ -807,9 +813,11 @@ def processStamps( rawStampsDict : `dict[str, Stamps]` ROI view of stamps from butler `guider_raw`. doSubtractMedian : `bool`, optional - If True, subtract median row bias. + If True, subtract column bias. columnMaskK : `float`, optional Threshold factor for column mask detection. + biasPercentile : `float`, optional + Percentile for column bias estimation. Returns ------- @@ -829,23 +837,28 @@ def processStamps( # Work on a copy - never modify original data = stamps[i].stamp_im.image.array.copy() - # Compute median row bias - medianRows = np.nanmedian(data, axis=0) + # Compute column bias using low percentile to avoid + # star contamination (median pulls dark dips at star cols). + percRows = np.nanpercentile(data, biasPercentile, axis=0) + rPerc = np.nanpercentile(data, biasPercentile) medianValue = np.nanmedian(data) - # Replace NaN medians (fully masked columns) with global median - medianRows = np.where(np.isnan(medianRows), medianValue, medianRows) + # Replace NaN values (fully masked columns) with global median + percRows = np.where(np.isnan(percRows), medianValue, percRows) # Compute column mask on bias-subtracted data - colMask = getColumnMask(data - medianRows[np.newaxis, :], k=columnMaskK) + colMask = getColumnMask( + data - percRows[np.newaxis, :] - (rPerc - medianValue), + k=columnMaskK, + ) # Create mask array for MaskedImageF (BAD=1 for masked cols) nRows, nCols = data.shape maskArray = np.zeros((nRows, nCols), dtype=np.uint32) maskArray[colMask] = 1 # Set BAD bit for masked columns - # Apply median row bias subtraction + # Apply column bias subtraction if doSubtractMedian: - data = data - medianRows[np.newaxis, :] + data = data - percRows[np.newaxis, :] - (rPerc - medianValue) # Fill masked columns with global median if colMask.any():