DM-49364: Switch LATISS/bestEffortIsr to IsrTaskLSST

pipelines/quickLook.yaml

@@ -4,3 +4,5 @@ tasks:
     class: lsst.summit.utils.quickLook.QuickLookIsrTask
     config:
       doRepairCosmics: true
+      doCorrectGains: false
+      doDeferredCharge: false

python/lsst/summit/utils/bestEffort.py

@@ -26,7 +26,7 @@
 import lsst.afw.image as afwImage
 import lsst.daf.butler as dafButler
 from lsst.daf.butler.registry import ConflictingDefinitionError
-from lsst.ip.isr import IsrTask
+from lsst.ip.isr import IsrTaskLSST
 from lsst.pex.config import Config
 from lsst.summit.utils.butlerUtils import getLatissDefaultCollections
 from lsst.summit.utils.quickLook import QuickLookIsrTask
@@ -58,7 +58,7 @@ class BestEffortIsr:
         Extra collections to add to the butler init. Collections are prepended.
     defaultExtraIsrOptions : `dict`, optional
         A dict of extra isr config options to apply. Each key should be an
-    attribute of an isrTaskConfigClass.
+    attribute of an isrTaskLSSTConfigClass.
     doRepairCosmics : `bool`, optional
         Repair cosmic ray hits?
     doWrite : `bool`, optional
@@ -238,7 +238,7 @@ def getExposure(
             raise RuntimeError(f"Failed to retrieve raw for exp {dataId}") from None
 
         # default options that are probably good for most engineering time
-        isrConfig = IsrTask.ConfigClass()
+        isrConfig = IsrTaskLSST.ConfigClass()
         with importlib.resources.path("lsst.summit.utils", "resources/config/quickLookIsr.py") as cfgPath:
             isrConfig.load(cfgPath)
 

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/quickLook.py

@@ -23,23 +23,21 @@
 import importlib.resources
 from typing import Any
 
-import numpy as np
-
 import lsst.afw.cameraGeom as camGeom
 import lsst.afw.image as afwImage
 import lsst.ip.isr as ipIsr
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
 import lsst.pipe.base.connectionTypes as cT
-from lsst.ip.isr import IsrTask
-from lsst.ip.isr.isrTask import IsrTaskConnections
+from lsst.ip.isr import IsrTaskLSST
+from lsst.ip.isr.isrTaskLSST import IsrTaskLSSTConnections
 from lsst.meas.algorithms.installGaussianPsf import InstallGaussianPsfTask
 from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask
 
 __all__ = ["QuickLookIsrTask", "QuickLookIsrTaskConfig"]
 
 
-class QuickLookIsrTaskConnections(IsrTaskConnections):
+class QuickLookIsrTaskConnections(IsrTaskLSSTConnections):
     """Copy isrTask's connections, changing prereq min values to zero.
 
     Copy all the connections directly for IsrTask, keeping ccdExposure as
@@ -50,7 +48,9 @@ class QuickLookIsrTaskConnections(IsrTaskConnections):
     def __init__(self, *, config: Any = None):
         # programatically clone all of the connections from isrTask
         # setting minimum values to zero for everything except the ccdExposure
-        super().__init__(config=IsrTask.ConfigClass())  # need a dummy config, isn't used other than for ctor
+        super().__init__(
+            config=IsrTaskLSST.ConfigClass()
+        )  # need a dummy config, isn't used other than for ctor
         for name, connection in self.allConnections.items():
             if hasattr(connection, "minimum"):
                 setattr(
@@ -95,13 +95,16 @@ class QuickLookIsrTask(pipeBase.PipelineTask):
     config: QuickLookIsrTaskConfig
     _DefaultName = "quickLook"
 
-    def __init__(self, isrTask: IsrTask = IsrTask, **kwargs: Any):
+    def __init__(self, isrTask: IsrTaskLSST = IsrTaskLSST, **kwargs: Any):
         super().__init__(**kwargs)
         # Pass in IsrTask so that we can modify it slightly for unit tests.
         # Note that this is not an instance of the IsrTask class, but the class
         # itself, which is then instantiated later on, in the run() method,
         # with the dynamically generated config.
-        self.isrTask = IsrTask
+        # import pdb; pdb.set_trace()
+        if IsrTaskLSST._DefaultName != "isrLSST":
+            raise RuntimeError("QuickLookIsrTask should now always use IsrTaskLSST for processing.")
+        self.isrTask = IsrTaskLSST
 
     def run(
         self,
@@ -111,22 +114,14 @@ def run(
         bias: afwImage.Exposure | None = None,
         dark: afwImage.Exposure | None = None,
         flat: afwImage.Exposure | None = None,
-        fringes: afwImage.Exposure | None = None,
         defects: ipIsr.Defects | None = None,
         linearizer: ipIsr.linearize.LinearizeBase | None = None,
         crosstalk: ipIsr.crosstalk.CrosstalkCalib | None = None,
-        bfKernel: np.ndarray | None = None,
-        newBFKernel: ipIsr.BrighterFatterKernel | None = None,
+        bfKernel: ipIsr.BrighterFatterKernel | None = None,
         ptc: ipIsr.PhotonTransferCurveDataset | None = None,
-        crosstalkSources: list | None = None,
-        isrBaseConfig: ipIsr.IsrTaskConfig | None = None,
-        filterTransmission: afwImage.TransmissionCurve | None = None,
-        opticsTransmission: afwImage.TransmissionCurve | None = None,
-        strayLightData: Any | None = None,
-        sensorTransmission: afwImage.TransmissionCurve | None = None,
-        atmosphereTransmission: afwImage.TransmissionCurve | None = None,
+        isrBaseConfig: ipIsr.IsrTaskLSSTConfig | None = None,
         deferredChargeCalib: Any | None = None,
-        illumMaskedImage: afwImage.MaskedImage | None = None,
+        gainCorrection: ipIsr.IsrCalib | None = None,
     ) -> pipeBase.Struct:
         """Run isr and cosmic ray repair using, doing as much isr as possible.
 
@@ -160,16 +155,14 @@ def run(
             Functor for linearization.
         crosstalk : `lsst.ip.isr.crosstalk.CrosstalkCalib`, optional
             Calibration for crosstalk.
-        bfKernel : `numpy.ndarray`, optional
-            Brighter-fatter kernel.
-        newBFKernel : `ipIsr.BrighterFatterKernel`, optional
+        bfKernel : `ipIsr.BrighterFatterKernel`, optional
             New Brighter-fatter kernel.
         ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
             Photon transfer curve dataset, with, e.g., gains
             and read noise.
-        crosstalkSources : `list`, optional
-            List of possible crosstalk sources.
-        isrBaseConfig : `lsst.ip.isr.IsrTaskConfig`, optional
+        cti : `lsst.ip.isr.DeferredChargeCalib`, optional
+            Charge transfer inefficiency correction calibration.
+        isrBaseConfig : `lsst.ip.isr.IsrTaskLSSTConfig`, optional
             An isrTask config to act as the base configuration. Options which
             involve applying a calibration product are ignored, but this allows
             for the configuration of e.g. the number of overscan columns.
@@ -195,6 +188,9 @@ def run(
             Gains used to override the detector's nominal gains for the
             brighter-fatter correction. A dict keyed by amplifier name for
             the detector in question.
+        gainCorrection : `lsst.ip.isr.IsrCalib`, optional
+            Gain correction calibration, to account for small
+            time-dependent shifts.
 
         Returns
         -------
@@ -204,7 +200,7 @@ def run(
                 The ISRed and cosmic-ray-repaired exposure.
         """
         if not isrBaseConfig:
-            isrConfig = IsrTask.ConfigClass()
+            isrConfig = IsrTaskLSST.ConfigClass()
             with importlib.resources.path("lsst.summit.utils", "resources/config/quickLookIsr.py") as cfgPath:
                 isrConfig.load(cfgPath)
         else:
@@ -213,12 +209,12 @@ def run(
         isrConfig.doBias = False
         isrConfig.doDark = False
         isrConfig.doFlat = False
-        isrConfig.doFringe = False
         isrConfig.doDefect = False
         isrConfig.doLinearize = False
         isrConfig.doCrosstalk = False
         isrConfig.doBrighterFatter = False
-        isrConfig.usePtcGains = False
+        isrConfig.doDeferredCharge = False
+        isrConfig.doCorrectGains = False
 
         if bias:
             isrConfig.doBias = True
@@ -232,10 +228,6 @@ def run(
             isrConfig.doFlat = True
             self.log.info("Running with flat correction")
 
-        if fringes:
-            isrConfig.doFringe = True
-            self.log.info("Running with fringe correction")
-
         if defects:
             isrConfig.doDefect = True
             self.log.info("Running with defect correction")
@@ -248,53 +240,46 @@ def run(
             isrConfig.doCrosstalk = True
             self.log.info("Running with crosstalk correction")
 
-        if newBFKernel is not None:
-            bfGains = newBFKernel.gain
-            isrConfig.doBrighterFatter = True
-            self.log.info("Running with new brighter-fatter correction")
-        else:
-            bfGains = None
-
-        if bfKernel is not None and bfGains is None:
+        if bfKernel is not None:
             isrConfig.doBrighterFatter = True
             self.log.info("Running with brighter-fatter correction")
 
-        if ptc:
-            isrConfig.usePtcGains = True
-            self.log.info("Running with ptc correction")
+        if deferredChargeCalib is not None:
+            isrConfig.doDeferredCharge = True
+            self.log.info("Running with CTI correction")
+
+        if gainCorrection is not None:
+            isrConfig.doCorrectGains = True
+            self.log.info("Running with Gain corrections")
+
+        if ptc is None:
+            raise RuntimeError("IsrTaskLSST requires a PTC.")
 
-        isrConfig.doWrite = False
         isrTask = self.isrTask(config=isrConfig)
 
-        if fringes:
-            # Must be run after isrTask is instantiated.
-            isrTask.fringe.loadFringes(
-                fringes,
-                expId=ccdExposure.info.id,
-                assembler=isrTask.assembleCcd if isrConfig.doAssembleIsrExposures else None,
-            )
+        # DM-47959: TODO Add fringe correction to IsrTaskLSST.
+        # if fringes:
+        #     # Must be run after isrTask is instantiated.
+        #     isrTask.fringe.loadFringes(
+        #         fringes,
+        #         expId=ccdExposure.info.id,
+        #         assembler=isrTask.assembleCcd if
+        # isrConfig.doAssembleIsrExposures else None,
+        #     )
 
         result = isrTask.run(
             ccdExposure,
             camera=camera,
             bias=bias,
             dark=dark,
             flat=flat,
-            fringes=fringes,
             defects=defects,
             linearizer=linearizer,
             crosstalk=crosstalk,
             bfKernel=bfKernel,
-            bfGains=bfGains,
             ptc=ptc,
-            crosstalkSources=crosstalkSources,
-            filterTransmission=filterTransmission,
-            opticsTransmission=opticsTransmission,
-            sensorTransmission=sensorTransmission,
-            atmosphereTransmission=atmosphereTransmission,
-            strayLightData=strayLightData,
             deferredChargeCalib=deferredChargeCalib,
-            illumMaskedImage=illumMaskedImage,
+            gainCorrection=gainCorrection,
         )
 
         postIsr = result.exposure

python/lsst/summit/utils/resources/config/quickLookIsr.py

@@ -1,8 +1,9 @@
 # mypy: disable-error-code="name-defined"
 
-config.doWrite = False  # this task writes separately, no need for this
 config.doSaturation = True  # saturation very important for roundness measurement in qfm
-config.doSaturationInterpolation = True
-config.overscan.fitType = "MEDIAN_PER_ROW"
-config.overscan.doParallelOverscan = True
 config.brighterFatterMaxIter = 2  # Uncomment this to remove test warning
+config.doDeferredCharge = False  # no calib for this yet
+config.doBootstrap = False
+config.doApplyGains = True
+config.doSuspect = False
+config.defaultSaturationSource = "CAMERAMODEL"

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)

tests/test_bestEffortIsr.py

@@ -40,7 +40,7 @@ def setUpClass(cls):
         # NCSA - LATISS/raw/all
         # TTS - LATISS-test-data-tts
         # summit - LATISS_test_data
-        cls.dataId = {"day_obs": 20210121, "seq_num": 743, "detector": 0}
+        cls.dataId = {"day_obs": 20250902, "seq_num": 10, "detector": 0}
 
     def test_getExposure(self):
         # in most locations this will load a pre-made image