Fix FootprintSet STDEV detection failure in runSourceDetection

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,

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:
         """

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)