diff --git a/python/lsst/summit/utils/guiders/detection.py b/python/lsst/summit/utils/guiders/detection.py index 4f17b3c2..a2c95bdc 100644 --- a/python/lsst/summit/utils/guiders/detection.py +++ b/python/lsst/summit/utils/guiders/detection.py @@ -37,8 +37,9 @@ from astropy.nddata import Cutout2D from astropy.stats import sigma_clipped_stats +import lsst.afw.detection as afwDetect from lsst.afw.image import ExposureF, ImageF, MaskedImageF -from lsst.summit.utils.utils import detectObjectsInExp +from lsst.afw.math import STDEVCLIP, makeStatistics from .reading import GuiderData @@ -277,6 +278,7 @@ def runSourceDetection( cutOutSize: int = 25, apertureRadius: int = 5, gain: float = 1.0, + nPixMin: int = 10, ) -> pd.DataFrame: """ Detect sources in an image and measure their properties. @@ -293,6 +295,8 @@ def runSourceDetection( Aperture radius in pixels for photometry. gain : `float` Detector gain (e-/ADU). + nPixMin : `int` + Minimum number of pixels in a footprint for detection. Returns ------- @@ -302,13 +306,26 @@ def runSourceDetection( # Step 1: Convert numpy image to MaskedImage and Exposure exposure = ExposureF(MaskedImageF(ImageF(image))) - # Step 2: Detect sources - # we assume that we have bright stars - # filter out stamps with no stars + # Step 2: Detect sources using STDEVCLIP for the background noise. + # The input coadd images are dithered (see GuiderData.getStampArrayCoadd) + # to prevent integer quantization from collapsing the pixel distribution + # and causing STDEVCLIP to return 0. (See DM-54263.) + footprints = None if not isBlankImage(image): - footprints = detectObjectsInExp(exposure, nSigma=threshold) - else: - footprints = None + median = np.nanmedian(image) + exposure.image -= median + imageStd = float(makeStatistics(exposure.getMaskedImage(), STDEVCLIP).getValue(STDEVCLIP)) + if imageStd <= 0: + # Fallback: sigma68 is robust to quantization and bright stars. + p16, p84 = np.nanpercentile(image, [16, 84]) + imageStd = (p84 - p16) / 2.0 + if imageStd <= 0: + exposure.image += median + return pd.DataFrame(columns=DEFAULT_COLUMNS) + absThreshold = threshold * imageStd + thresh = afwDetect.Threshold(absThreshold, afwDetect.Threshold.VALUE) + footprints = afwDetect.FootprintSet(exposure.getMaskedImage(), thresh, "DETECTED", nPixMin) + exposure.image += median if not footprints: return pd.DataFrame(columns=DEFAULT_COLUMNS) @@ -613,7 +630,6 @@ def buildReferenceCatalog( apertureRadius = int(config.aperSizeArcsec / pixelScale) array = guiderData.getStampArrayCoadd(guiderName) - # array = np.where(array < 0, 0, array) # Ensure no negative values array = array - np.nanmin(array) # Ensure no negative values sources = runSourceDetection( array, diff --git a/python/lsst/summit/utils/guiders/reading.py b/python/lsst/summit/utils/guiders/reading.py index d671fb5f..2ca65206 100644 --- a/python/lsst/summit/utils/guiders/reading.py +++ b/python/lsst/summit/utils/guiders/reading.py @@ -406,10 +406,15 @@ def getStampArrayCoadd(self, detName: str) -> np.ndarray: if len(stamps) == 0: raise ValueError(f"No stamps found for detector {detName!r}") - # Collect arrays, with optional bias subtraction + # Collect arrays, with optional bias subtraction. arrList = [self[detName, idx] for idx in range(len(stamps))] - stack = np.nanmedian(arrList, axis=0) - return stack + # Add a uniform [0, 1) dither per stamp before the median to break + # integer quantization. Without this, the coadd pixel distribution + # can be so narrow that STDEVCLIP returns 0. (See DM-54263.) + stack = np.array(arrList, dtype=np.float32) + rng = np.random.default_rng(seed=0) + stack += rng.uniform(0, 1, size=stack.shape).astype(np.float32) + return np.nanmedian(stack, axis=0) def getGuiderAmpName(self, detName: str) -> str: """ diff --git a/python/lsst/summit/utils/utils.py b/python/lsst/summit/utils/utils.py index 471ef6ea..ec8750e1 100644 --- a/python/lsst/summit/utils/utils.py +++ b/python/lsst/summit/utils/utils.py @@ -488,6 +488,7 @@ def _getAltAzZenithsFromSeqNum( for seqNum in seqNumList: md = butler.get("raw.metadata", day_obs=dayObs, seq_num=seqNum, detector=0) obsInfo = ObservationInfo(md) + assert obsInfo.altaz_begin is not None, f"altaz_begin is None for seqNum={seqNum}" alt = obsInfo.altaz_begin.alt.value az = obsInfo.altaz_begin.az.value elevations.append(alt)