diff --git a/python/lsst/ip/isr/ampOffset.py b/python/lsst/ip/isr/ampOffset.py index 89fe1b9af..ca248a9d4 100644 --- a/python/lsst/ip/isr/ampOffset.py +++ b/python/lsst/ip/isr/ampOffset.py @@ -28,7 +28,7 @@ from lsst.afw.table import SourceTable from lsst.meas.algorithms import SourceDetectionTask, SubtractBackgroundTask from lsst.pex.config import Config, ConfigurableField, Field -from lsst.pipe.base import Task +from lsst.pipe.base import Task, Struct class AmpOffsetConfig(Config): @@ -73,12 +73,13 @@ def setDefaults(self): dtype=float, default=5.0, ) - ampEdgeWindow = Field( - doc="Pixel size of the sliding window used to generate rolling average amp offset values. It should " - "should be reconfigured for every instrument (HSC, LSSTCam, etc.) and should not exceed the shorter " - "dimension of amps. If not provided, it defaults to the shorter dimension of amps.", - dtype=int, - default=0, + ampEdgeWindowFrac = Field( + doc="Fraction of the amp edge lengths utilized as the sliding window for generating rolling average " + "amp offset values. It should be reconfigured for every instrument (HSC, LSSTCam, etc.) and should " + "not exceed 1. If not provided, it defaults to the fraction that recovers the pixel size of the " + "sliding window used in obs_subaru for compatibility with existing HSC data.", + dtype=float, + default=512 / 4176, ) doBackground = Field( doc="Estimate and subtract background prior to amp offset estimation?", @@ -118,6 +119,8 @@ def __init__(self, *args, **kwargs): self.makeSubtask("background") if self.config.doDetection: self.makeSubtask("detection") + # Initialize all of the instance variables here. + self.shortAmpSide = 0 def run(self, exposure): """Calculate amp offset values, determine corrective pedestals for each @@ -133,9 +136,14 @@ def run(self, exposure): exp = exposure.clone() bitMask = exp.mask.getPlaneBitMask(self.background.config.ignoredPixelMask) amps = exp.getDetector().getAmplifiers() + + # Check that all amps have the same gemotry. + ampDims = [amp.getBBox().getDimensions() for amp in amps] + if not all(dim == ampDims[0] for dim in ampDims): + raise RuntimeError("All amps should have the same geometry.") + # Assuming all the amps have the same geometry. - self.shortAmpSide = np.min(amps[0].getBBox().getDimensions()) - self.longAmpSide = np.max(amps[0].getBBox().getDimensions()) + self.shortAmpSide = np.min(ampDims[0]) # Fit and subtract background. if self.config.doBackground: @@ -165,18 +173,18 @@ def run(self, exposure): # Safety check: do any pixels remain for amp offset estimation? if (exp.mask.array & bitMask).all(): self.log.warning("All pixels masked: cannot calculate any amp offset corrections.") + pedestals = np.zeros(len(amps)) else: # Set up amp offset inputs. im = exp.image im.array[(exp.mask.array & bitMask) > 0] = np.nan - if self.config.ampEdgeWindow > self.shortAmpSide: + if self.config.ampEdgeWindowFrac > 1: raise RuntimeError( - f"The sliding window ({self.config.ampEdgeWindow} pixels) is larger than the smaller " - f"dimension in amplifiers ({self.shortAmpSide} pixels). This can lead to downstream " - "issues after convolution in the `getSideAmpOffset()` method. Please adjust " - f"`ampEdgeWindow` in the config to a value equal to or less than {self.shortAmpSide} " - "pixels and rerun." + f"The specified fraction (`ampEdgeWindowFrac`={self.config.ampEdgeWindowFrac}) of the " + "edge length exceeds 1. This leads to complications downstream, after convolution in " + "the `getSideAmpOffset()` method. Please modify the `ampEdgeWindowFrac` value in the " + "config to be 1 or less and rerun." ) # Determine amplifier geometry. @@ -193,17 +201,20 @@ def run(self, exposure): # If least-squares minimization fails, convert NaNs to zeroes, # ensuring that no values are erroneously added/subtracted. pedestals = np.nan_to_num(np.linalg.lstsq(A, B, rcond=None)[0]) - metadata = exposure.getMetadata() - for amp, pedestal in zip(amps, pedestals): - ampIm = exposure.image[amp.getBBox()].array - ampIm -= pedestal - ampName = amp.getName() - metadata.set( - f"LSST ISR AMPOFFSET PEDESTAL {ampName}", - float(pedestal), - f"Pedestal level subtracted from amp {ampName}" - ) - self.log.info(f"amp pedestal values: {', '.join([f'{x:.4f}' for x in pedestals])}") + + metadata = exposure.getMetadata() + for amp, pedestal in zip(amps, pedestals): + ampIm = exposure.image[amp.getBBox()].array + ampIm -= pedestal + ampName = amp.getName() + metadata.set( + f"LSST ISR AMPOFFSET PEDESTAL {ampName}", + float(pedestal), + f"Pedestal level subtracted from amp {ampName}", + ) + self.log.info(f"amp pedestal values: {', '.join([f'{x:.4f}' for x in pedestals])}") + + return Struct(pedestals=pedestals) def getAmpAssociations(self, amps): """Determine amp geometry and amp associations from a list of @@ -303,66 +314,6 @@ def getNeighbors(self, ampIds, ampId): sides.append(side) return neighbors, sides - def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB): - """Calculate the amp offset for a given interface between two - amplifiers. - - Parameters - ---------- - ampIdA : int - ID of the first amplifier. - ampIdB : int - ID of the second amplifier. - edgeA : numpy.ndarray - Amp edge for the first amplifier. - edgeB : numpy.ndarray - Amp edge for the second amplifier. - - Returns - ------- - interfaceOffset : float - The calculated amp offset value for the given interface between - amps A and B. - """ - sctrl = StatisticsControl() - # NOTE: Taking the difference with the order below fixes the sign flip - # in the B matrix. - edgeDiff = edgeA - edgeB - if self.config.ampEdgeWindow == 0: - if len(edgeDiff) > self.shortAmpSide: - window = self.shortAmpSide - else: - window = int(self.shortAmpSide**2 / self.longAmpSide) - self.log.info( - f"ampEdgeWindow automatically set to {window} pixels, for the " - f"amp {ampIdA} interfacing amp {ampIdB}, as value of 0 or no " - "input is provided for ampEdgeWindow." - ) - else: - window = self.config.ampEdgeWindow - # Compute rolling averages. - edgeDiffSum = np.convolve(np.nan_to_num(edgeDiff), np.ones(window), "same") - edgeDiffNum = np.convolve(~np.isnan(edgeDiff), np.ones(window), "same") - edgeDiffAvg = edgeDiffSum / np.clip(edgeDiffNum, 1, None) - edgeDiffAvg[np.isnan(edgeDiff)] = np.nan - # Take clipped mean of rolling average data as amp offset value. - interfaceOffset = makeStatistics(edgeDiffAvg, MEANCLIP, sctrl).getValue() - # Perform a couple of do-no-harm safety checks: - # a) The fraction of unmasked pixel rows is > ampEdgeMinFrac, - # b) The absolute offset ADU value is < ampEdgeMaxOffset. - ampEdgeGoodFrac = 1 - (np.sum(np.isnan(edgeDiffAvg)) / len(edgeDiffAvg)) - minFracFail = ampEdgeGoodFrac < self.config.ampEdgeMinFrac - maxOffsetFail = np.abs(interfaceOffset) > self.config.ampEdgeMaxOffset - if minFracFail or maxOffsetFail: - interfaceOffset = 0 - interfaceId = f"{ampIdA}{ampIdB}" - self.log.debug( - f"amp interface {interfaceId} : " - f"viable edge difference frac = {ampEdgeGoodFrac}, " - f"interface offset = {interfaceOffset:.3f}" - ) - return interfaceOffset - def getAmpOffsets(self, im, amps, associations, sides): """Calculate the amp offsets for all amplifiers. @@ -448,3 +399,57 @@ def getAmpEdges(self, im, amps, ampSides): warnings.filterwarnings("ignore", r"All-NaN (slice|axis) encountered") ampEdges[ampId][ampSide] = np.nanmedian(strip, axis=ampSide % 2) # 1D medianified strip return ampEdges + + def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB): + """Calculate the amp offset for a given interface between two + amplifiers. + + Parameters + ---------- + ampIdA : int + ID of the first amplifier. + ampIdB : int + ID of the second amplifier. + edgeA : numpy.ndarray + Amp edge for the first amplifier. + edgeB : numpy.ndarray + Amp edge for the second amplifier. + + Returns + ------- + interfaceOffset : float + The calculated amp offset value for the given interface between + amps A and B. + """ + interfaceId = f"{ampIdA}{ampIdB}" + sctrl = StatisticsControl() + # NOTE: Taking the difference with the order below fixes the sign flip + # in the B matrix. + edgeDiff = edgeA - edgeB + window = int(self.config.ampEdgeWindowFrac * len(edgeDiff)) + # Compute rolling averages. + edgeDiffSum = np.convolve(np.nan_to_num(edgeDiff), np.ones(window), "same") + edgeDiffNum = np.convolve(~np.isnan(edgeDiff), np.ones(window), "same") + edgeDiffAvg = edgeDiffSum / np.clip(edgeDiffNum, 1, None) + edgeDiffAvg[np.isnan(edgeDiff)] = np.nan + # Take clipped mean of rolling average data as amp offset value. + interfaceOffset = makeStatistics(edgeDiffAvg, MEANCLIP, sctrl).getValue() + # Perform a couple of do-no-harm safety checks: + # a) The fraction of unmasked pixel rows is > ampEdgeMinFrac, + # b) The absolute offset ADU value is < ampEdgeMaxOffset. + ampEdgeGoodFrac = 1 - (np.sum(np.isnan(edgeDiffAvg)) / len(edgeDiffAvg)) + minFracFail = ampEdgeGoodFrac < self.config.ampEdgeMinFrac + maxOffsetFail = np.abs(interfaceOffset) > self.config.ampEdgeMaxOffset + if minFracFail or maxOffsetFail: + interfaceOffset = 0 + self.log.warning( + f"The fraction of unmasked pixels for amp interface {interfaceId} is below the threshold " + f"({self.config.ampEdgeMinFrac}) or the absolute offset value exceeds the limit " + f"({self.config.ampEdgeMaxOffset} ADU). Setting the interface offset to 0." + ) + self.log.debug( + f"amp interface {interfaceId} : " + f"viable edge difference frac = {ampEdgeGoodFrac}, " + f"interface offset = {interfaceOffset:.3f}" + ) + return interfaceOffset diff --git a/tests/test_ampOffset.py b/tests/test_ampOffset.py index 1ecc08108..69d8e4bc5 100644 --- a/tests/test_ampOffset.py +++ b/tests/test_ampOffset.py @@ -22,79 +22,128 @@ import numpy as np import lsst.utils.tests -import lsst.afw.image as afwImage -from lsst.obs.lsst import LsstCam -from lsst.obs.subaru import HyperSuprimeCam from lsst.ip.isr.ampOffset import AmpOffsetConfig, AmpOffsetTask +from lsst.ip.isr.isrMock import IsrMock class AmpOffsetTest(lsst.utils.tests.TestCase): def setUp(self): - self.camera = None # `lsst.afw.cameraGeom.Camera` object - self.backgroundFractionSample = None # not used for now - self.ampEdgeWindow = None # not used for now + # Test with a single detector featuring 4x2 amplifiers to ensure + # functionality in a general 2-dimensional scenario. + config = IsrMock.ConfigClass() + config.isLsstLike = True + config.doAddBias = False + config.doAddDark = False + config.doAddFlat = False + config.doAddFringe = False + config.doGenerateImage = True + config.doGenerateData = True + config.doGenerateAmpDict = True + self.mock = IsrMock(config=config) + self.measuredValuesBackground = [ + 0.13336708, + 0.06299729, + 0.34589145, + 0.27551733, + -0.27551731, + -0.34589142, + -0.06299731, + -0.1333671, + ] + self.measuredValuesRampBackground = [ + 0.19094621, + 0.21447723, + 0.19441773, + 0.21793666, + -0.21794104, + -0.1944167, + -0.21447818, + -0.19094191, + ] + self.measuredSigmaBackground = 1.7861559064582404 + self.measuredSigmaRampBackground = 1.8148713749702268 def tearDown(self): - del self.camera + del self.mock def buildExposure(self, addBackground=False): - exp = afwImage.ExposureF(self.camera[0].getBBox()) - exp.setDetector(self.camera[0]) - amps = exp.getDetector().getAmplifiers() + exp = self.mock.getExposure() + detector = exp.getDetector() + amps = detector.getAmplifiers() self.values = np.linspace(-2.5, 2.5, len(amps)) - for i, amp in enumerate(amps): - exp.image[amp.getBBox()] = self.values[i] + for amp, value in zip(amps, self.values): + exp.image[amp.getBBox()] = value if addBackground: exp.image.array += 100 return exp def testAmpOffset(self): - if self.camera is None: - raise unittest.SkipTest("Skipping due to no camera set.") exp = self.buildExposure(addBackground=False) - amps = exp.getDetector().getAmplifiers() config = AmpOffsetConfig() config.doBackground = False config.doDetection = False task = AmpOffsetTask(config=config) - task.run(exp) - meta = exp.getMetadata().toDict() + pedestals = task.run(exp).pedestals self.assertEqual(np.sum(exp.image.array), 0) - for amp, value in zip(amps, self.values): - ampName = amp.getName() - self.assertAlmostEqual(meta[f"LSST ISR AMPOFFSET PEDESTAL {ampName}"], value, 6) + for pedestal, value in zip(pedestals, self.values): + self.assertAlmostEqual(pedestal, value, 6) def testAmpOffsetWithBackground(self): - if self.camera is None: - raise unittest.SkipTest("Skipping due to no camera set.") exp = self.buildExposure(addBackground=True) amps = exp.getDetector().getAmplifiers() config = AmpOffsetConfig() config.doBackground = True config.doDetection = True task = AmpOffsetTask(config=config) - task.run(exp) - meta = exp.getMetadata().toDict() + pedestals = task.run(exp).pedestals nAmps = len(amps) - ampNames = [] - for amp in amps: - ampNames.append(amp.getName()) for i in range(nAmps // 2): - self.assertAlmostEqual( - meta[f"LSST ISR AMPOFFSET PEDESTAL {ampNames[i]}"], - -meta[f"LSST ISR AMPOFFSET PEDESTAL {ampNames[nAmps-i-1]}"], - 7, - ) - - -class AmpOffsetTestSubaru(AmpOffsetTest): - def setUp(self): - self.camera = HyperSuprimeCam().getCamera() - - -class AmpOffsetTestRubin(AmpOffsetTest): - def setUp(self): - self.camera = LsstCam().getCamera() + self.assertAlmostEqual(pedestals[i], -pedestals[nAmps - i - 1], 5) + for pedestal, value in zip(pedestals, self.measuredValuesBackground): + self.assertAlmostEqual(pedestal, value, 5) + # If we are getting it wrong, let's not get it wrong by more than some + # DN. + self.assertLessEqual(np.std(pedestals - self.values), self.measuredSigmaBackground, 5) + + def testAmpOffsetWithRampBackground(self): + exp = self.buildExposure(addBackground=True) + amps = exp.getDetector().getAmplifiers() + yscale = 100.0 + xscale = 0.0 + # Add a gradient. + self.amplifierAddYGradient(exp.image, 0.0, yscale) + # Add another gradient to the other direction. + self.amplifierAddXGradient(exp.image, 0.0, xscale) + config = AmpOffsetConfig() + config.doBackground = True + config.doDetection = True + task = AmpOffsetTask(config=config) + pedestals = task.run(exp).pedestals + nAmps = len(amps) + for i in range(nAmps // 2): + self.assertAlmostEqual(pedestals[i], -pedestals[nAmps - i - 1], 5) + for pedestal, value in zip(pedestals, self.measuredValuesRampBackground): + self.assertAlmostEqual(pedestal, value, 5) + self.assertLessEqual(np.std(pedestals - self.values), self.measuredSigmaRampBackground, 5) + + # The two static methods below are taken from ip_isr/isrMock. + @staticmethod + def amplifierAddYGradient(ampData, start, end): + nPixY = ampData.getDimensions().getY() + ampArr = ampData.array + ampArr[:] = ampArr[:] + ( + np.interp(range(nPixY), (0, nPixY - 1), (start, end)).reshape(nPixY, 1) + + np.zeros(ampData.getDimensions()).transpose() + ) + + @staticmethod + def amplifierAddXGradient(ampData, start, end): + nPixX = ampData.getDimensions().getX() + ampArr = ampData.array + ampArr[:] = ampArr[:] + ( + np.interp(range(nPixX), (0, nPixX - 1), (start, end)).reshape(1, nPixX) + + np.zeros(ampData.getDimensions()).transpose() + ) class TestMemory(lsst.utils.tests.MemoryTestCase):