From 011fb4bbd1e99cc4c700b1470e377603618bf2d7 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Wed, 9 Oct 2024 08:30:03 -0700 Subject: [PATCH 1/3] Add ability to undo bgModel1 in SkyCorrectionTask --- python/lsst/pipe/tasks/skyCorrection.py | 94 ++++++++++++++++++------- 1 file changed, 70 insertions(+), 24 deletions(-) diff --git a/python/lsst/pipe/tasks/skyCorrection.py b/python/lsst/pipe/tasks/skyCorrection.py index df6f00f1c..6292dd23a 100644 --- a/python/lsst/pipe/tasks/skyCorrection.py +++ b/python/lsst/pipe/tasks/skyCorrection.py @@ -178,6 +178,11 @@ class SkyCorrectionConfig(PipelineTaskConfig, pipelineConnections=SkyCorrectionC dtype=FocalPlaneBackgroundConfig, doc="Initial background model, prior to sky frame subtraction", ) + undoBgModel1 = Field( + dtype=bool, + default=False, + doc="If True, adds back initial background model after sky and removes bgModel1 from the list", + ) sky = ConfigurableField( target=SkyMeasurementTask, doc="Sky measurement", @@ -210,6 +215,11 @@ def setDefaults(self): self.bgModel2.ySize = 256 self.bgModel2.smoothScale = 1.0 + def validate(self): + super().validate() + if self.undoBgModel1 and not self.doSky and not self.doBgModel2: + raise ValueError("If undoBgModel1 is True, task requires at least one of doSky or doBgModel2.") + class SkyCorrectionTask(PipelineTask): """Perform a full focal plane sky correction.""" @@ -278,11 +288,11 @@ def run(self, calExps, calBkgs, skyFrames, camera): Parameters ---------- - calExps : `list` [`lsst.afw.image.exposure.ExposureF`] + calExps : `list` [`lsst.afw.image.ExposureF`] Detector calibrated exposure images for the visit. calBkgs : `list` [`lsst.afw.math.BackgroundList`] Detector background lists matching the calibrated exposures. - skyFrames : `list` [`lsst.afw.image.exposure.ExposureF`] + skyFrames : `list` [`lsst.afw.image.ExposureF`] Sky frame calibration data for the input detectors. camera : `lsst.afw.cameraGeom.Camera` Camera matching the input data to process. @@ -292,24 +302,29 @@ def run(self, calExps, calBkgs, skyFrames, camera): results : `Struct` containing: skyCorr : `list` [`lsst.afw.math.BackgroundList`] Detector-level sky correction background lists. - calExpMosaic : `lsst.afw.image.exposure.ExposureF` + calExpMosaic : `lsst.afw.image.ExposureF` Visit-level mosaic of the sky corrected data, binned. Analogous to `calexp - skyCorr`. - calBkgMosaic : `lsst.afw.image.exposure.ExposureF` + calBkgMosaic : `lsst.afw.image.ExposureF` Visit-level mosaic of the sky correction background, binned. Analogous to `calexpBackground + skyCorr`. """ # Restore original backgrounds in-place; optionally refine mask maps numOrigBkgElements = [len(calBkg) for calBkg in calBkgs] - _ = self._restoreBackgroundRefineMask(calExps, calBkgs) + _ = self._restoreOriginalBackgroundRefineMask(calExps, calBkgs) # Bin exposures, generate full-fp bg, map to CCDs and subtract in-place _ = self._subtractVisitBackground(calExps, calBkgs, camera, self.config.bgModel1) + initialBackgroundIndex = len(calBkgs[0]._backgrounds) - 1 # Subtract a scaled sky frame from all input exposures if self.config.doSky: self._subtractSkyFrame(calExps, skyFrames, calBkgs) + # Adds full-fp bg back onto exposures, removes it from list + if self.config.undoBgModel1: + _ = self._undoInitialBackground(calExps, calBkgs, initialBackgroundIndex) + # Bin exposures, generate full-fp bg, map to CCDs and subtract in-place if self.config.doBgModel2: _ = self._subtractVisitBackground(calExps, calBkgs, camera, self.config.bgModel2) @@ -328,7 +343,7 @@ def run(self, calExps, calBkgs, skyFrames, camera): return Struct(skyCorr=calBkgs, calExpMosaic=calExpMosaic, calBkgMosaic=calBkgMosaic) - def _restoreBackgroundRefineMask(self, calExps, calBkgs): + def _restoreOriginalBackgroundRefineMask(self, calExps, calBkgs): """Restore original background to each calexp and invert the related background model; optionally refine the mask plane. @@ -350,17 +365,17 @@ def _restoreBackgroundRefineMask(self, calExps, calBkgs): Parameters ---------- - calExps : `lsst.afw.image.exposure.ExposureF` + calExps : `lsst.afw.image.ExposureF` Detector level calexp images to process. - calBkgs : `lsst.afw.math._backgroundList.BackgroundList` + calBkgs : `lsst.afw.math.BackgroundList` Detector level background lists associated with the calexps. Returns ------- - calExps : `lsst.afw.image.exposure.ExposureF` - The calexps with the initially subtracted background restored. - skyCorrBases : `lsst.afw.math._backgroundList.BackgroundList` - The inverted initial background models; the genesis for skyCorrs. + calExps : `lsst.afw.image.ExposureF` + The calexps with the originally subtracted background restored. + skyCorrBases : `lsst.afw.math.BackgroundList` + The inverted original background models; the genesis for skyCorrs. """ skyCorrBases = [] for calExp, calBkg in zip(calExps, calBkgs): @@ -379,7 +394,7 @@ def _restoreBackgroundRefineMask(self, calExps, calBkgs): stats = np.nanpercentile(skyCorrBase.array, [50, 75, 25]) self.log.info( - "Detector %d: Initial background restored; BG median = %.1f counts, BG IQR = %.1f counts", + "Detector %d: Original background restored; BG median = %.1f counts, BG IQR = %.1f counts", calExp.getDetector().getId(), -stats[0], np.subtract(*stats[1:]), @@ -387,6 +402,37 @@ def _restoreBackgroundRefineMask(self, calExps, calBkgs): skyCorrBases.append(skyCorrBase) return calExps, skyCorrBases + def _undoInitialBackground(self, calExps, calBkgs, initialBackgroundIndex): + """Undo the initial background subtraction (bgModel1) after sky frame + subtraction. + + Parameters + ---------- + calExps : `list` [`lsst.afw.image.ExposureF`] + Calibrated exposures to be background subtracted. + calBkgs : `list` [`lsst.afw.math.BackgroundList`] + Background lists associated with the input calibrated exposures. + initialBackgroundIndex : `int` + Index of the initial background (bgModel1) in the background list. + + Notes + ----- + Inputs are modified in-place. + """ + for calExp, calBkg in zip(calExps, calBkgs): + image = calExp.getMaskedImage() + + # Remove bgModel1 from the background list; restore in the image + initialBackground = calBkg[initialBackgroundIndex][0].getImageF() + image += initialBackground + calBkg._backgrounds.pop(initialBackgroundIndex) + + self.log.info( + "Detector %d: The initial background model prior to sky frame subtraction (bgModel1) has " + "been removed from the background list", + calExp.getDetector().getId(), + ) + def _subtractVisitBackground(self, calExps, calBkgs, camera, config): """Perform a full focal-plane background subtraction for a visit. @@ -400,9 +446,9 @@ def _subtractVisitBackground(self, calExps, calBkgs, camera, config): Parameters ---------- - calExps : `list` [`lsst.afw.image.exposure.ExposureF`] + calExps : `list` [`lsst.afw.image.ExposureF`] Calibrated exposures to be background subtracted. - calBkgs : `list` [`lsst.afw.math._backgroundList.BackgroundList`] + calBkgs : `list` [`lsst.afw.math.BackgroundList`] Background lists associated with the input calibrated exposures. camera : `lsst.afw.cameraGeom.Camera` Camera description. @@ -413,7 +459,7 @@ def _subtractVisitBackground(self, calExps, calBkgs, camera, config): ------- calExps : `list` [`lsst.afw.image.maskedImage.MaskedImageF`] Background subtracted exposures for creating a focal plane image. - calBkgs : `list` [`lsst.afw.math._backgroundList.BackgroundList`] + calBkgs : `list` [`lsst.afw.math.BackgroundList`] Updated background lists with a visit-level model appended. """ # Set up empty full focal plane background model object @@ -478,16 +524,16 @@ def _subtractDetectorBackground(self, calExp, bgModel): Parameters ---------- - calExp : `lsst.afw.image.exposure.ExposureF` + calExp : `lsst.afw.image.ExposureF` Exposure to subtract the background model from. bgModel : `lsst.pipe.tasks.background.FocalPlaneBackground` Full focal plane camera-level background model. Returns ------- - calExp : `lsst.afw.image.exposure.ExposureF` + calExp : `lsst.afw.image.ExposureF` Background subtracted input exposure. - calBkgElement : `lsst.afw.math._backgroundList.BackgroundList` + calBkgElement : `lsst.afw.math.BackgroundList` Detector level realization of the full focal plane bg model. """ image = calExp.getMaskedImage() @@ -510,11 +556,11 @@ def _subtractSkyFrame(self, calExps, skyFrames, calBkgs): Parameters ---------- - calExps : `list` [`lsst.afw.image.exposure.ExposureF`] + calExps : `list` [`lsst.afw.image.ExposureF`] Calibrated exposures to be background subtracted. - skyFrames : `list` [`lsst.afw.image.exposure.ExposureF`] + skyFrames : `list` [`lsst.afw.image.ExposureF`] Sky frame calibration data for the input detectors. - calBkgs : `list` [`lsst.afw.math._backgroundList.BackgroundList`] + calBkgs : `list` [`lsst.afw.math.BackgroundList`] Background lists associated with the input calibrated exposures. """ skyFrameBgModels = [] @@ -549,11 +595,11 @@ def _binAndMosaic(self, exposures, camera, binning, ids=None, refExps=None): Binning size to be applied to input images. ids : `list` [`int`], optional List of detector ids to iterate over. - refExps : `list` [`lsst.afw.image.exposure.ExposureF`], optional + refExps : `list` [`lsst.afw.image.ExposureF`], optional If supplied, mask planes from these reference images will be used. Returns ------- - mosaicImage : `lsst.afw.image.exposure.ExposureF` + mosaicImage : `lsst.afw.image.ExposureF` Mosaicked full focal plane image. """ refExps = np.resize(refExps, len(exposures)) # type: ignore From 6f0c8e1474530778575c65ebe3a98295a1c7a318 Mon Sep 17 00:00:00 2001 From: Lee Kelvin Date: Wed, 18 Dec 2024 05:19:45 -0800 Subject: [PATCH 2/3] Add skyFrameScale as a returned output --- python/lsst/pipe/tasks/skyCorrection.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/python/lsst/pipe/tasks/skyCorrection.py b/python/lsst/pipe/tasks/skyCorrection.py index 6292dd23a..340e39cb7 100644 --- a/python/lsst/pipe/tasks/skyCorrection.py +++ b/python/lsst/pipe/tasks/skyCorrection.py @@ -300,6 +300,8 @@ def run(self, calExps, calBkgs, skyFrames, camera): Returns ------- results : `Struct` containing: + skyFrameScale : `float` + Scale factor applied to the sky frame. skyCorr : `list` [`lsst.afw.math.BackgroundList`] Detector-level sky correction background lists. calExpMosaic : `lsst.afw.image.ExposureF` @@ -318,8 +320,9 @@ def run(self, calExps, calBkgs, skyFrames, camera): initialBackgroundIndex = len(calBkgs[0]._backgrounds) - 1 # Subtract a scaled sky frame from all input exposures + skyFrameScale = None if self.config.doSky: - self._subtractSkyFrame(calExps, skyFrames, calBkgs) + skyFrameScale = self._subtractSkyFrame(calExps, skyFrames, calBkgs) # Adds full-fp bg back onto exposures, removes it from list if self.config.undoBgModel1: @@ -341,7 +344,9 @@ def run(self, calExps, calBkgs, skyFrames, camera): skyCorrExtras, camera, self.config.binning, ids=calExpIds, refExps=calExps ) - return Struct(skyCorr=calBkgs, calExpMosaic=calExpMosaic, calBkgMosaic=calBkgMosaic) + return Struct( + skyFrameScale=skyFrameScale, skyCorr=calBkgs, calExpMosaic=calExpMosaic, calBkgMosaic=calBkgMosaic + ) def _restoreOriginalBackgroundRefineMask(self, calExps, calBkgs): """Restore original background to each calexp and invert the related @@ -562,6 +567,11 @@ def _subtractSkyFrame(self, calExps, skyFrames, calBkgs): Sky frame calibration data for the input detectors. calBkgs : `list` [`lsst.afw.math.BackgroundList`] Background lists associated with the input calibrated exposures. + + Returns + ------- + scale : `float` + Scale factor applied to the sky frame. """ skyFrameBgModels = [] scales = [] @@ -577,6 +587,7 @@ def _subtractSkyFrame(self, calExps, skyFrames, calBkgs): # also updating the calBkg list in-place self.sky.subtractSkyFrame(calExp.getMaskedImage(), skyFrameBgModel, scale, calBkg) self.log.info("Sky frame subtracted with a scale factor of %.5f", scale) + return scale def _binAndMosaic(self, exposures, camera, binning, ids=None, refExps=None): """Bin input exposures and mosaic across the entire focal plane. From 56ab5dde596e9f3adbd970865d54fb4ba753e9c4 Mon Sep 17 00:00:00 2001 From: Lee Kelvin Date: Wed, 18 Dec 2024 05:20:51 -0800 Subject: [PATCH 3/3] Add unit test for SkyCorrectionTask --- tests/test_skyCorrection.py | 160 ++++++++++++++++++++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 tests/test_skyCorrection.py diff --git a/tests/test_skyCorrection.py b/tests/test_skyCorrection.py new file mode 100644 index 000000000..aa7e06541 --- /dev/null +++ b/tests/test_skyCorrection.py @@ -0,0 +1,160 @@ +# This file is part of pipe_tasks. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import unittest +from copy import deepcopy + +import lsst.utils.tests +import numpy as np +from lsst.afw.image import ExposureF +from lsst.afw.math import BackgroundMI +from lsst.obs.base.instrument_tests import DummyCam +from lsst.pipe.tasks.skyCorrection import SkyCorrectionConfig, SkyCorrectionTask + + +class SkyCorrectionTestCase(lsst.utils.tests.TestCase): + + def setUp(self): + dummyCam = DummyCam() + self.camera = dummyCam.getCamera() + bbox = self.camera[0].getBBox() + + # Configs below set to approximate HSC defaults + self.skyCorrectionConfig = SkyCorrectionConfig() + self.skyCorrectionConfig.doMaskObjects = False + # Set bgModel1 size to a single bin for the whole plane (aka constant) + self.skyCorrectionConfig.bgModel1.xSize = 8192 * 0.015 + self.skyCorrectionConfig.bgModel1.ySize = 8192 * 0.015 + self.skyCorrectionConfig.bgModel1.pixelSize = 0.015 + self.skyCorrectionConfig.bgModel2.xSize = 256 * 0.015 + self.skyCorrectionConfig.bgModel2.ySize = 256 * 0.015 + self.skyCorrectionConfig.bgModel2.pixelSize = 0.015 + + # Generate calexp/calexpBackground/sky for all detectors + self.calExps = [] + self.calBkgs = [] + self.skyFrames = [] + self.background_level = 3000 + self.sky_level = 5 + for detector in [0, 1]: + rng = np.random.default_rng(detector) + + # Science image + calexp = ExposureF(bbox) + calexp.maskedImage.set(0.0, 0x0, 650.0) + calexp.setDetector(self.camera[detector]) + # Add a sky frame signature to a subregion of the image + sky_frame_bin_size = 32 + x_start = 32 * sky_frame_bin_size + x_stop = 64 * sky_frame_bin_size + y_start = 31 * sky_frame_bin_size + y_stop = 63 * sky_frame_bin_size + calexp.image.array[:, x_start:x_stop] += self.sky_level + calexp.image.array[y_start:y_stop, :] += self.sky_level + # Add random noise + calexp.image.array += rng.normal(0.0, 25.0, (bbox.getDimensions().y, bbox.getDimensions().x)) + self.calExps.append(calexp) + + # Background image + backgroundImage = ExposureF(bbox) + backgroundImage.maskedImage.set(0.0, 0x0, 1.0) + backgroundImage.setDetector(self.camera[detector]) + backgroundImage.image.array += self.background_level + background = BackgroundMI(bbox, backgroundImage.getMaskedImage()) + calexpBackground = lsst.afw.math.BackgroundList( + ( + background, + lsst.afw.math.Interpolate.CONSTANT, + lsst.afw.math.UndersampleStyle.REDUCE_INTERP_ORDER, + lsst.afw.math.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) + self.calBkgs.append(calexpBackground) + + # Sky frame + sky = ExposureF(128, 125) + sky.maskedImage.set(0.0, 0x0, 1.0) + sky.setDetector(self.camera[detector]) + header = sky.getMetadata() + header.set("BOX.MINX", bbox.getMinX()) + header.set("BOX.MINY", bbox.getMinY()) + header.set("BOX.MAXX", bbox.getMaxX()) + header.set("BOX.MAXY", bbox.getMaxY()) + header.set("ALGORITHM", "NATURAL_SPLINE") + sky.image.array[:, 32:64] += 1 # x + sky.image.array[31:63, :] += 1 # y + # Add random noise + sky.image.array += rng.normal(0.0, 0.1, (125, 128)) + sky.image.array -= np.sum(sky.image.array) / (125 * 128) + self.skyFrames.append(sky) + + def tearDown(self): + del self.camera + del self.calExps + del self.calBkgs + del self.skyFrames + + def testSkyCorrectionDefault(self): + """Test SkyCorrectionTask with mostly default configuration values.""" + + skyCorrectionTask = SkyCorrectionTask(config=self.skyCorrectionConfig) + # Pass in deep copies, as the task modifies the input data + results = skyCorrectionTask.run( + deepcopy(self.calExps), deepcopy(self.calBkgs), self.skyFrames, self.camera + ) + skyFrameScale = results.skyFrameScale + skyCorr = results.skyCorr + self.assertEqual(len(skyCorr), len(self.calExps)) + self.assertAlmostEqual(skyFrameScale, self.sky_level, delta=1e-1) + self.assertAlmostEqual(np.nanmean(results.calExpMosaic.array), 0, delta=1e-2) + + def testSkyCorrectionSkyFrameOnly(self): + """Test SkyCorrectionTask with the config undoBgModel1 set to True.""" + + skyCorrectionConfig = deepcopy(self.skyCorrectionConfig) + skyCorrectionConfig.undoBgModel1 = True + skyCorrectionConfig.doBgModel2 = False + skyCorrectionTask = SkyCorrectionTask(config=skyCorrectionConfig) + # Pass in deep copies, as the task modifies the input data + results = skyCorrectionTask.run( + deepcopy(self.calExps), deepcopy(self.calBkgs), self.skyFrames, self.camera + ) + self.assertAlmostEqual( + np.nanmean(results.calExpMosaic.array), + np.nanmean(self.calExps[0].image.array) + self.background_level, + delta=1e-2, + ) + + +class MemoryTester(lsst.utils.tests.MemoryTestCase): + pass + + +def setup_module(module): + lsst.utils.tests.init() + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main()