Coverage for tests / test_isrFunctions.py: 12%
537 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:10 +0000
1#
2# LSST Data Management System
3# Copyright 2008-2017 AURA/LSST.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
21#
23import unittest
24import logging
25import numpy as np
27from lsst.daf.base import PropertyList
28import lsst.geom as geom
29import lsst.afw.image as afwImage
30import lsst.afw.detection as afwDetection
31import lsst.utils.tests
32import lsst.ip.isr as ipIsr
33import lsst.ip.isr.isrFunctions as ipIsrFunctions
34import lsst.ip.isr.isrMock as isrMock
37def computeImageMedianAndStd(image):
38 """Function to calculate median and std of image data.
40 Parameters
41 ----------
42 image : `lsst.afw.image.Image`
43 Image to measure statistics on.
45 Returns
46 -------
47 median : `float`
48 Image median.
49 std : `float`
50 Image stddev.
51 """
52 median = np.nanmedian(image.getArray())
53 std = np.nanstd(image.getArray())
55 return (median, std)
58def _setExposureSatCore(exposure, x, y, halfWidthX, halfWidthY,
59 satVal, satMaskBit):
60 """Helper function to make a rectangular saturated core.
62 Parameters
63 ----------
64 exposure : `lsst.afw.image.Exposure`
65 Input exposure.
66 x : `int`
67 X coordinate for the center of the rectangular saturated core.
68 y : `int`
69 Y coordinate for the center of the rectangular saturated core.
70 halfWidthX : `int`
71 Half the width of the rectangular saturated core along the x axis.
72 halfWidthY : `int`
73 Half the width of the rectangular saturated core along the y axis.
74 satVal : `float`
75 Value the saturated core is set to.
76 satMaskBit : `int`
77 Saturated mask bit the saturated core mask is set to.
78 """
79 xmax = exposure.image.array.shape[1]
80 lowerRange = x-halfWidthX
81 upperRange = x+halfWidthX
82 if lowerRange < 0:
83 lowerRange = 0
84 if upperRange > xmax:
85 upperRange = xmax
86 exposure.image.array[y-halfWidthY:y+halfWidthY,
87 lowerRange:upperRange] = satVal
88 exposure.mask.array[y-halfWidthY:y+halfWidthY,
89 lowerRange:upperRange] = satMaskBit
92def _setExposureSatColumns(exposure, x, y, halfWidthX, limY,
93 satVal, satMaskBit, isTop=False):
94 """Helper function to add saturated columns.
95 We need to add these around a saturated core to test a typical edge bleed.
97 Parameters
98 ----------
99 exposure : `lsst.afw.image.Exposure`
100 Input exposure.
101 x : `int`
102 X coordinate for the central saturated column.
103 y : `int`
104 Y coordinate for the central saturated column.
105 halfWidthX : `int`
106 Half the number of saturated columns.
107 limY : `int`
108 Y coordinate up to which the saturated columns go to.
109 satVal : `float`
110 Value the saturated columns are set to.
111 satMaskBit : `int`
112 Saturated mask bit the saturated columns mask is set to.
113 isTop : `bool`
114 True if the saturated columns go toward the top of the detector.
115 """
116 if isTop:
117 exposure.image.array[y:limY,
118 x-halfWidthX:x+halfWidthX] = satVal
119 exposure.mask.array[y:limY,
120 x-halfWidthX:x+halfWidthX] = satMaskBit
121 else:
122 exposure.image.array[limY:y,
123 x-halfWidthX:x+halfWidthX] = satVal
124 exposure.mask.array[limY:y,
125 x-halfWidthX:x+halfWidthX] = satMaskBit
128def _makeEdgeBleed(exposure, x, extentY, edgeBleedWidth,
129 edgeBleedConstant, edgeBleedWidthLimit,
130 saturationLevel, saturationFrac,
131 satVal, satMaskBit, isTop=False):
132 """Helper function to add an edge bleed with a decaying exponential
133 model.
135 Parameters
136 ----------
137 exposure : `lsst.afw.image.Exposure`
138 Input exposure.
139 x : `int`
140 X coordinate for the center of the edge bleed (generally the same
141 as the center of the saturated core and columns).
142 extentY : `int`
143 Y coordinate up to which the edge bleed goes to (generally the same
144 as limY of the saturated columns).
145 edgeBleedWidth : `int`
146 Width of the edge bleed at the edge of the detector.
147 edgeBleedConstant : `float`
148 Constant for the decaying exponential model.
149 edgeBleedWidthLimit : `int`
150 The width the edge bleed goes to away from the edge (generally the
151 width of the saturated columns).
152 saturationLevel : `float`
153 Saturation level.
154 saturationFrac : `float`
155 The inside of the edge bleed is set to a value equal to this fraction
156 of the saturation level.
157 satVal : `float`
158 The value the contour of the edge bleed is set to, greater or equal to
159 the saturation level.
160 satMaskBit : `int`
161 Saturated mask bit the contour of the edge bleed mask is set to.
162 isTop : `bool`
163 True if the edge bleed is a the top edge of the detector.
164 """
165 xmax = exposure.image.array.shape[1]
166 for y in range(extentY):
167 edgeBleedWidthY = edgeBleedWidth*np.exp(-edgeBleedConstant*y) \
168 + edgeBleedWidthLimit
169 if isTop:
170 # For edge bleed in top amplifier
171 y = exposure.image.array.shape[0]-1-y
172 lowerRange = x-int(edgeBleedWidthY/2.)
173 upperRange = x+int(edgeBleedWidthY/2.)
174 if lowerRange < 0:
175 lowerRange = 0
176 if upperRange >= xmax:
177 upperRange = xmax-1
178 exposure.image.array[y, lowerRange] = satVal
179 exposure.mask.array[y, lowerRange] = satMaskBit
180 exposure.image.array[y, upperRange] = satVal
181 exposure.mask.array[y, upperRange] = satMaskBit
182 exposure.image.array[y,
183 lowerRange+1:upperRange
184 ] = saturationFrac*saturationLevel
187class MockE2VAmp:
188 def __init__(self, name, bbox):
189 self._name = name
190 self._bbox = bbox
192 def getName(self):
193 return self._name
195 def getBBox(self):
196 return self._bbox
198 def __repr__(self):
199 return f"MockE2VAmp({self._name})"
202class MockE2VDetector(list):
203 def __init__(self):
204 amps = []
205 for i in range(8):
206 name = f"C1{i}"
207 bbox = geom.Box2I(corner=geom.Point2I(i*512, 2002), dimensions=geom.Extent2I(512, 2002))
208 amps.append(MockITLAmp(name, bbox))
209 for i in reversed(range(8)):
210 name = f"C0{i}"
211 bbox = geom.Box2I(corner=geom.Point2I(i*512, 0), dimensions=geom.Extent2I(512, 2002))
212 amps.append(MockITLAmp(name, bbox))
214 super().__init__(amps)
216 def getBBox(self):
217 return geom.Box2I(corner=geom.Point2I(0, 0), dimensions=geom.Extent2I(4096, 4004))
220class MockE2VExposure(afwImage.ExposureF):
221 def __init__(self, *args, **kwargs):
222 super().__init__(*args, **kwargs)
224 def setDetector(self, detector):
225 self._detector = detector
227 def fillSaturationMetadata(self, saturationLevel):
228 if not self._detector:
229 raise RuntimeError("No detector set.")
230 for amp in self.getDetector():
231 self.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"] = saturationLevel
233 def getDetector(self):
234 return self._detector
237class MockITLAmp:
238 def __init__(self, name, bbox):
239 self._name = name
240 self._bbox = bbox
242 def getName(self):
243 return self._name
245 def getBBox(self):
246 return self._bbox
248 def __repr__(self):
249 return f"MockITLAmp({self._name})"
252class MockITLDetector(list):
253 def __init__(self):
254 amps = []
255 for i in range(8):
256 name = f"C1{i}"
257 bbox = geom.Box2I(corner=geom.Point2I(i*509, 2000), dimensions=geom.Extent2I(509, 2000))
258 amps.append(MockITLAmp(name, bbox))
259 for i in reversed(range(8)):
260 name = f"C0{i}"
261 bbox = geom.Box2I(corner=geom.Point2I(i*509, 0), dimensions=geom.Extent2I(509, 2000))
262 amps.append(MockITLAmp(name, bbox))
264 super().__init__(amps)
266 def getBBox(self):
267 return geom.Box2I(corner=geom.Point2I(0, 0), dimensions=geom.Extent2I(4072, 4000))
270class MockITLExposure(afwImage.ExposureF):
271 def __init__(self, *args, **kwargs):
272 super().__init__(*args, **kwargs)
274 # Monkey-patch our fake detector.
275 def setDetector(self, detector):
276 self._detector = detector
278 def fillSaturationMetadata(self, saturationLevel):
279 if not self._detector:
280 raise RuntimeError("No detector set.")
281 for amp in self.getDetector():
282 self.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"] = saturationLevel
284 def getDetector(self):
285 return self._detector
288class IsrFunctionsCases(lsst.utils.tests.TestCase):
289 """Test that functions for ISR produce expected outputs.
290 """
291 def setUp(self):
292 self.inputExp = isrMock.TrimmedRawMock().run()
293 self.mi = self.inputExp.getMaskedImage()
295 def test_transposeMaskedImage(self):
296 """Expect height and width to be exchanged.
297 """
298 transposed = ipIsr.transposeMaskedImage(self.mi)
299 self.assertEqual(transposed.getImage().getBBox().getHeight(),
300 self.mi.getImage().getBBox().getWidth())
301 self.assertEqual(transposed.getImage().getBBox().getWidth(),
302 self.mi.getImage().getBBox().getHeight())
304 def test_interpolateDefectList(self):
305 """Expect number of interpolated pixels to be non-zero.
306 """
307 defectList = isrMock.DefectMock().run()
308 self.assertEqual(len(defectList), 1)
310 for fallbackValue in (None, -999.0):
311 for haveMask in (True, False):
312 with self.subTest(fallbackValue=fallbackValue, haveMask=haveMask):
313 if haveMask is False:
314 if 'INTRP' in self.mi.getMask().getMaskPlaneDict():
315 self.mi.getMask().removeAndClearMaskPlane('INTRP')
316 else:
317 if 'INTRP' not in self.mi.getMask().getMaskPlaneDict():
318 self.mi.getMask().addMaskPlane('INTRP')
319 numBit = ipIsr.countMaskedPixels(self.mi, "INTRP")
320 self.assertEqual(numBit, 0)
322 def test_transposeDefectList(self):
323 """Expect bbox dimension values to flip.
324 """
325 defectList = isrMock.DefectMock().run()
326 transposed = defectList.transpose()
328 for d, t in zip(defectList, transposed):
329 self.assertEqual(d.getBBox().getDimensions().getX(), t.getBBox().getDimensions().getY())
330 self.assertEqual(d.getBBox().getDimensions().getY(), t.getBBox().getDimensions().getX())
332 def test_makeThresholdMask(self):
333 """Expect list of defects to have elements.
334 """
335 defectList = ipIsr.makeThresholdMask(self.mi, 200,
336 growFootprints=2,
337 maskName='SAT')
339 self.assertEqual(len(defectList), 1)
341 def test_E2VEdgeBleedMask(self):
342 """Test expected number of masked pixels according to e2v edge bleed
343 masking.
344 """
345 detector = MockE2VDetector()
346 exposure = MockE2VExposure(detector.getBBox())
347 exposure.setDetector(detector)
348 exposure.mask.array[:, :] = 0
349 satMaskBit = exposure.mask.getPlaneBitMask('SAT')
351 saturationLevel = 12000.
352 exposure.fillSaturationMetadata(saturationLevel)
354 # add rectangular edge bleed to mock exposure
355 _setExposureSatCore(exposure, x=1600, y=1000, halfWidthX=200, halfWidthY=100, satVal=saturationLevel,
356 satMaskBit=satMaskBit)
357 _setExposureSatColumns(exposure, x=1600, y=1000, halfWidthX=4, limY=0, satVal=saturationLevel,
358 satMaskBit=satMaskBit, isTop=False)
360 exposure.image.array[:100, 1600:1800] = 10000.
362 numPixSatBottomEdgeBefore = len(np.where(exposure.mask.array[0, :] == satMaskBit)[0])
364 ipIsrFunctions.maskE2VEdgeBleed(exposure, e2vEdgeBleedSatMinArea=20000, e2vEdgeBleedSatMaxArea=100000,
365 e2vEdgeBleedYMax=350,
366 saturatedMaskName='SAT')
368 numPixSatBottomEdgeAfter = len(np.where(exposure.mask.array[0, :] == satMaskBit)[0])
370 # Check the number of saturated pixels
371 self.assertEqual(ipIsr.countMaskedPixels(exposure, 'SAT'), 263600)
372 # Check there are more pixels along the bottom edge after masking
373 self.assertGreater(numPixSatBottomEdgeAfter, numPixSatBottomEdgeBefore)
375 def test_ITLEdgeBleedMask(self):
376 """Expect number of masked pixels according to edge bleed masking.
377 """
378 detector = MockITLDetector()
379 exposure = MockITLExposure(detector.getBBox())
380 exposure.setDetector(detector)
381 exposure.mask.array[:, :] = 0
382 satMaskBit = exposure.mask.getPlaneBitMask('SAT')
384 saturationLevel = 12000.
385 exposure.fillSaturationMetadata(saturationLevel)
387 badAmpDict = {}
388 for amp in exposure.getDetector():
389 if amp.getName() == 'C10':
390 badAmpDict[amp.getName()] = True
391 else:
392 badAmpDict[amp.getName()] = False
394 # Case 1:
395 # Create two edge bleeds that fall into the same footprint.
397 # A saturated core in the bottom amplifier
398 # that has a bottom edge bleed
399 satVal = 15000.
400 x = 70
401 y = 120
402 halfWidthX = 40
403 halfWidthY = 70
404 halfWidthColX = 3
405 limY = 100
406 extentY = 100
407 _setExposureSatCore(exposure, x=x, y=y,
408 halfWidthX=halfWidthX, halfWidthY=halfWidthY,
409 satVal=satVal, satMaskBit=satMaskBit)
410 _setExposureSatColumns(exposure, x=x, y=y, halfWidthX=halfWidthColX,
411 limY=limY, satVal=satVal, satMaskBit=satMaskBit,
412 isTop=False)
413 edgeBleedWidth = 60
414 edgeBleedConstant = 0.04 # we make the edge bleed a bit smaller than the masking model
415 edgeBleedWidthLimit = 2*halfWidthColX
416 saturationFrac = 0.85 # higher than default in masking model
417 _makeEdgeBleed(exposure, x=x, extentY=extentY,
418 edgeBleedWidth=edgeBleedWidth,
419 edgeBleedConstant=edgeBleedConstant,
420 edgeBleedWidthLimit=edgeBleedWidthLimit,
421 saturationLevel=saturationLevel,
422 saturationFrac=saturationFrac,
423 satVal=satVal, satMaskBit=satMaskBit, isTop=False)
425 # A saturated core near the left edge of the detector, close to the
426 # previous edge bleed
427 x = 10
428 y = 90
429 halfWidthX = 14
430 halfWidthY = 50
431 halfWidthColX = 3
432 _setExposureSatCore(exposure, x=x, y=y, halfWidthX=halfWidthX,
433 halfWidthY=halfWidthY, satVal=satVal,
434 satMaskBit=satMaskBit)
435 # Bottom saturated columns to connect core and edge bleed
436 limY = 100
437 _setExposureSatColumns(exposure, x=x, y=y, halfWidthX=halfWidthColX,
438 limY=limY, satVal=satVal, satMaskBit=satMaskBit,
439 isTop=False)
440 edgeBleedWidth = 70 # this would go beyond the edge of the detector
441 edgeBleedWidthLimit = 2*halfWidthColX
442 _makeEdgeBleed(exposure, x=x, extentY=extentY,
443 edgeBleedWidth=edgeBleedWidth,
444 edgeBleedConstant=edgeBleedConstant,
445 edgeBleedWidthLimit=edgeBleedWidthLimit,
446 saturationLevel=saturationLevel,
447 saturationFrac=saturationFrac,
448 satVal=satVal, satMaskBit=satMaskBit, isTop=False)
450 thresh = afwDetection.Threshold(exposure.mask.getPlaneBitMask("SAT"),
451 afwDetection.Threshold.BITMASK
452 )
453 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
455 satAreas = np.asarray([fp.getArea() for fp in fpList])
456 largeAreas, = np.where((satAreas >= 1000)
457 & (satAreas < 100000))
459 # We select the only footprint
460 fpCore = fpList[largeAreas[0]]
462 numPixSatBottomEdgeBeforeCase1 = len(np.where(exposure.mask.array[0, :] == satMaskBit)[0])
463 ipIsr.maskITLEdgeBleed(exposure, badAmpDict,
464 fpCore,
465 itlEdgeBleedSatMinArea=2400,
466 itlEdgeBleedSatMaxArea=100000,
467 itlEdgeBleedThreshold=5000.,
468 itlEdgeBleedModelConstant=0.02,
469 saturatedMaskName='SAT')
470 numPixSatBottomEdgeAfterCase1 = len(np.where(exposure.mask.array[0, :] == satMaskBit)[0])
471 # Check the number of saturated pixels
472 self.assertEqual(ipIsr.countMaskedPixels(exposure, 'SAT'), 22967)
473 # Check there are more pixels along the bottom edge after masking
474 self.assertGreater(numPixSatBottomEdgeAfterCase1, numPixSatBottomEdgeBeforeCase1)
476 # Case 2:
477 # A saturated core in the top amplifier
478 # that has a top edge bleed
479 x = 2100
480 y = 3100
481 halfWidthX = 50
482 halfWidthY = 80
483 halfWidthColX = 4
484 limY = exposure.image.array.shape[0]-100
485 _setExposureSatCore(exposure, x=x, y=y, halfWidthX=halfWidthX,
486 halfWidthY=halfWidthY, satVal=satVal,
487 satMaskBit=satMaskBit)
488 _setExposureSatColumns(exposure, x=x, y=y, halfWidthX=halfWidthColX,
489 limY=limY, satVal=satVal, satMaskBit=satMaskBit,
490 isTop=True)
491 edgeBleedWidth = 100
492 edgeBleedWidthLimit = 2*halfWidthColX
493 _makeEdgeBleed(exposure, x=x, extentY=extentY,
494 edgeBleedWidth=edgeBleedWidth,
495 edgeBleedConstant=edgeBleedConstant,
496 edgeBleedWidthLimit=edgeBleedWidthLimit,
497 saturationLevel=saturationLevel,
498 saturationFrac=saturationFrac,
499 satVal=satVal, satMaskBit=satMaskBit, isTop=True)
501 # Add a saturation sag column around the saturated core
502 xSatSag = x
503 yminSatSag = y - halfWidthY - 100
504 ymaxSatSag = y - halfWidthY
505 exposure.image.array[yminSatSag:ymaxSatSag, xSatSag] = satVal*0.7
507 # Re-measure the saturated footprints
508 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
510 satAreas = np.asarray([fp.getArea() for fp in fpList])
511 largeAreas, = np.where((satAreas >= 10000)
512 & (satAreas < 100000))
514 for largeAreasIndex in largeAreas:
515 # We select the new edge bleed footprint
516 if fpList[largeAreasIndex].getBBox().contains(x, y):
517 fpCore = fpList[largeAreasIndex]
519 numPixSatTopEdgeBeforeCase2 = len(np.where(exposure.mask.array[-1, :] == satMaskBit)[0])
520 ipIsr.maskITLEdgeBleed(exposure, badAmpDict,
521 fpCore, itlEdgeBleedSatMinArea=10000,
522 itlEdgeBleedSatMaxArea=100000,
523 itlEdgeBleedThreshold=5000.,
524 itlEdgeBleedModelConstant=0.02,
525 saturatedMaskName='SAT')
526 numPixSatTopEdgeAfterCase2 = len(np.where(exposure.mask.array[-1, :] == satMaskBit)[0])
527 # Check the number of saturated pixels
528 self.assertEqual(ipIsr.countMaskedPixels(exposure, 'SAT'), 50745)
529 # Check there are more pixels along the bottom edge
530 # after masking
531 self.assertGreater(numPixSatTopEdgeAfterCase2, numPixSatTopEdgeBeforeCase2)
533 # This will mask the whole column containing the saturation sag
534 ipIsrFunctions.maskITLSatSag(exposure, fpCore, saturatedMaskName='SAT')
535 numPixColumnMask = len(np.where(exposure.mask.array[:, xSatSag] == satMaskBit)[0])
536 # Test that the saturation sag column is masked
537 self.assertEqual(numPixColumnMask, exposure.image.array.shape[0])
539 # Case 3:
540 # A saturated core with an edge bleed on both edges
541 x = 3020
542 y = 1000
543 halfWidthX = 70
544 halfWidthY = 150
545 halfWidthColX = 7
546 _setExposureSatCore(exposure, x=x, y=y, halfWidthX=halfWidthX,
547 halfWidthY=halfWidthY, satVal=satVal,
548 satMaskBit=satMaskBit)
549 # Top saturated columns to connect core and edge bleed
550 limY = exposure.image.array.shape[0]-100
551 _setExposureSatColumns(exposure, x=x, y=y, halfWidthX=halfWidthColX,
552 limY=limY, satVal=satVal, satMaskBit=satMaskBit,
553 isTop=True)
554 # Bottom saturated columns to connect core and edge bleed
555 limY = 100
556 _setExposureSatColumns(exposure, x=x, y=y, halfWidthX=halfWidthColX,
557 limY=limY, satVal=satVal, satMaskBit=satMaskBit,
558 isTop=False)
559 edgeBleedWidth = 150
560 edgeBleedWidthLimit = 2*halfWidthColX
561 _makeEdgeBleed(exposure, x=x, extentY=extentY,
562 edgeBleedWidth=edgeBleedWidth,
563 edgeBleedConstant=edgeBleedConstant,
564 edgeBleedWidthLimit=edgeBleedWidthLimit,
565 saturationLevel=saturationLevel,
566 saturationFrac=saturationFrac,
567 satVal=satVal, satMaskBit=satMaskBit, isTop=True)
568 _makeEdgeBleed(exposure, x=x, extentY=extentY,
569 edgeBleedWidth=edgeBleedWidth,
570 edgeBleedConstant=edgeBleedConstant,
571 edgeBleedWidthLimit=edgeBleedWidthLimit,
572 saturationLevel=saturationLevel,
573 saturationFrac=saturationFrac,
574 satVal=satVal, satMaskBit=satMaskBit, isTop=False)
576 # Re-measure the saturated footprints
577 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
579 satAreas = np.asarray([fp.getArea() for fp in fpList])
580 largeAreas, = np.where((satAreas >= 10000)
581 & (satAreas < 100000))
583 for largeAreasIndex in largeAreas:
584 # We select the new edge bleed footprint
585 if fpList[largeAreasIndex].getBBox().contains(x, y):
586 fpCore = fpList[largeAreasIndex]
588 # Number of saturated pixels at the bottom edge
589 # before applying masking
590 numPixSatBottomEdgeBefore = len(np.where(exposure.mask.array[0, :] == satMaskBit)[0])
591 # Number of saturated pixels at the top edge
592 # before applying masking
593 numPixSatTopEdgeBefore = len(np.where(exposure.mask.array[-1, :] == satMaskBit)[0])
594 # Apply edge bleed masking
595 ipIsr.maskITLEdgeBleed(exposure, badAmpDict,
596 fpCore,
597 itlEdgeBleedSatMinArea=10000,
598 itlEdgeBleedSatMaxArea=100000,
599 itlEdgeBleedThreshold=5000.,
600 itlEdgeBleedModelConstant=0.02,
601 saturatedMaskName='SAT')
602 # Number of saturated pixels at the bottom edge
603 # after applying edge bleed masking
604 numPixSatBottomEdgeAfter = len(np.where(exposure.mask.array[0, :] == satMaskBit)[0])
605 numPixSatTopEdgeAfter = len(np.where(exposure.mask.array[-1, :] == satMaskBit)[0])
607 # Check the number of saturated pixels
608 self.assertEqual(ipIsr.countMaskedPixels(exposure, 'SAT'), 184669)
609 # Check there are more pixels along the bottom edge
610 # after masking
611 self.assertGreater(numPixSatBottomEdgeAfter, numPixSatBottomEdgeBefore)
612 # Same with top edge
613 self.assertGreater(numPixSatTopEdgeAfter, numPixSatTopEdgeBefore)
615 # Case 4:
616 # A saturated core near the right edge of the detector
617 x = 4060
618 y = 3000
619 halfWidthX = 30
620 halfWidthY = 500
621 halfWidthColX = 3
622 _setExposureSatCore(exposure, x=x, y=y, halfWidthX=halfWidthX,
623 halfWidthY=halfWidthY, satVal=satVal,
624 satMaskBit=satMaskBit)
625 # Top saturated columns to connect core and edge bleed
626 limY = exposure.image.array.shape[0]-100
627 _setExposureSatColumns(exposure, x=x, y=y, halfWidthX=halfWidthColX,
628 limY=limY, satVal=satVal, satMaskBit=satMaskBit,
629 isTop=True)
630 edgeBleedWidth = 100 # this would go beyond the edge of the detector
631 edgeBleedWidthLimit = 2*halfWidthColX
632 _makeEdgeBleed(exposure, x=x, extentY=extentY,
633 edgeBleedWidth=edgeBleedWidth,
634 edgeBleedConstant=edgeBleedConstant,
635 edgeBleedWidthLimit=edgeBleedWidthLimit,
636 saturationLevel=saturationLevel,
637 saturationFrac=saturationFrac,
638 satVal=satVal, satMaskBit=satMaskBit, isTop=True)
640 # Re-measure the saturated footprints
641 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
643 satAreas = np.asarray([fp.getArea() for fp in fpList])
644 largeAreas, = np.where((satAreas >= 10000)
645 & (satAreas < 100000))
647 for largeAreasIndex in largeAreas:
648 # We select the new edge bleed footprint
649 if fpList[largeAreasIndex].getBBox().contains(x, y):
650 fpCore = fpList[largeAreasIndex]
652 # Number of saturated pixels at the bottom edge
653 # before applying masking
654 numPixSatTopEdgeBefore = len(np.where(exposure.mask.array[-1, :] == satMaskBit)[0])
655 # Apply edge bleed masking
656 ipIsr.maskITLEdgeBleed(exposure, badAmpDict,
657 fpCore,
658 itlEdgeBleedSatMinArea=10000,
659 itlEdgeBleedSatMaxArea=100000,
660 itlEdgeBleedThreshold=5000.,
661 itlEdgeBleedModelConstant=0.02,
662 saturatedMaskName='SAT')
663 # Number of saturated pixels at the bottom edge
664 # after applying edge bleed masking
665 numPixSatTopEdgeAfter = len(np.where(exposure.mask.array[-1, :] == satMaskBit)[0])
667 # Check the number of saturated pixels
668 self.assertEqual(ipIsr.countMaskedPixels(exposure, 'SAT'), 238379)
669 # Check there are more pixels along the bottom edge
670 # after masking
671 self.assertGreater(numPixSatTopEdgeAfter, numPixSatTopEdgeBefore)
673 def test_itlDipMasking(self):
674 """Test the ITL dip masking."""
675 detector = MockITLDetector()
676 exposure = MockITLExposure(detector.getBBox())
677 exposure.setDetector(detector)
678 exposure.mask.array[:, :] = 0
679 exposure.mask.addMaskPlane("ITL_DIP")
680 dipMaskValue = exposure.mask.getPlaneBitMask(["SUSPECT", "ITL_DIP"])
682 detectorConfig = ipIsr.overscanAmpConfig.OverscanDetectorConfig()
683 detectorConfig.itlDipBackgroundFraction = 0.0025
685 # Set the background to 1000 electrons.
686 exposure.image.array[:, :] = 1000.0
688 # Add some saturated masks, below and above the trigger.
689 satMaskValue = exposure.mask.getPlaneBitMask("SAT")
691 # Above the threshold in width/height.
692 exposure.mask.array[
693 500: 500 + detectorConfig.itlDipMinHeight * 2,
694 1000: 1000 + detectorConfig.itlDipMinWidth * 2
695 ] |= satMaskValue
697 # At the threshold in width/height.
698 exposure.mask.array[
699 1200: 1200 + detectorConfig.itlDipMinHeight + 1,
700 1500: 1500 + detectorConfig.itlDipMinWidth + 1
701 ] |= satMaskValue
703 # Below the threshold in width.
704 exposure.mask.array[
705 1600: 1600 + detectorConfig.itlDipMinHeight + 1,
706 2500: 2500 + detectorConfig.itlDipMinWidth // 2
707 ] |= satMaskValue
709 # Below the threshold in height.
710 exposure.mask.array[
711 2500: 2500 + detectorConfig.itlDipMinHeight // 2,
712 3000: 3000 + detectorConfig.itlDipMinWidth + 1
713 ] |= satMaskValue
715 with self.assertLogs(level=logging.INFO) as cm:
716 ipIsr.isrFunctions.maskITLDip(exposure, detectorConfig)
717 self.assertEqual(len(cm[1]), 2)
718 self.assertIn("Found ITL dip (width 30; bkg 1000.00); masking column 992 to 1037", cm[1][0])
719 self.assertIn("Found ITL dip (width 16; bkg 1000.00); masking column 1495 to 1519", cm[1][1])
721 # This includes the scaled edges
722 np.testing.assert_array_equal(exposure.mask.array[:, 992: 1038] & dipMaskValue, dipMaskValue)
723 np.testing.assert_array_equal(exposure.mask.array[:, 991] & dipMaskValue, 0)
724 np.testing.assert_array_equal(exposure.mask.array[:, 1038] & dipMaskValue, 0)
726 np.testing.assert_array_equal(exposure.mask.array[:, 1495: 1520] & dipMaskValue, dipMaskValue)
727 np.testing.assert_array_equal(exposure.mask.array[:, 1494] & dipMaskValue, 0)
728 np.testing.assert_array_equal(exposure.mask.array[:, 1520] & dipMaskValue, 0)
730 # The other two should not be masked.
731 np.testing.assert_array_equal(exposure.mask.array[:, 2500] & dipMaskValue, 0)
732 np.testing.assert_array_equal(exposure.mask.array[:, 3000] & dipMaskValue, 0)
734 # Change the background to a much lower value, only 1 should be masked.
735 exposure.mask.array[:, :] &= ~dipMaskValue
736 exposure.image.array[:, :] = 50.0
738 with self.assertLogs(level=logging.INFO) as cm:
739 ipIsr.isrFunctions.maskITLDip(exposure, detectorConfig)
740 self.assertEqual(len(cm[1]), 1)
741 self.assertIn("Found ITL dip (width 30; bkg 50.00); masking column 992 to 1037", cm[1][0])
743 # This one should be masked.
744 np.testing.assert_array_equal(exposure.mask.array[:, 992: 1038] & dipMaskValue, dipMaskValue)
745 np.testing.assert_array_equal(exposure.mask.array[:, 991] & dipMaskValue, 0)
746 np.testing.assert_array_equal(exposure.mask.array[:, 1038] & dipMaskValue, 0)
748 # The other three should not be masked.
749 np.testing.assert_array_equal(exposure.mask.array[:, 1500] & dipMaskValue, 0)
750 np.testing.assert_array_equal(exposure.mask.array[:, 2500] & dipMaskValue, 0)
751 np.testing.assert_array_equal(exposure.mask.array[:, 3000] & dipMaskValue, 0)
753 # And blow things up so that masking is not applied.
754 # This needs to be in several sections because of the max width.
755 # We additionally avoid the edge. Note that there are rounding
756 # offsets that occur due to the scaling such that this setting
757 # will mask 42 columns per block or a total of 504.
758 exposure.mask.array[:, :] = 0
759 detectorConfig.itlDipWidthScale = 1.0
760 for i in range(12):
761 exposure.mask.array[:, (i + 1)*100: (i + 1)*100 + 41] |= satMaskValue
763 maskBak = exposure.mask.array.copy()
764 with self.assertLogs(level=logging.WARNING) as cm:
765 ipIsr.isrFunctions.maskITLDip(exposure, detectorConfig)
766 self.assertEqual(len(cm[1]), 1)
767 self.assertIn("Too many (504) columns would be masked", cm[1][0])
768 np.testing.assert_array_equal(exposure.mask.array, maskBak)
770 def test_interpolateFromMask(self):
771 """Expect number of interpolated pixels to be non-zero.
772 """
773 ipIsr.makeThresholdMask(self.mi, 200, growFootprints=2,
774 maskName='SAT')
775 for growFootprints in range(0, 3):
776 for useLegacyInterp in (False, True):
777 interpMaskedImage = ipIsr.interpolateFromMask(self.mi, 2.0,
778 growSaturatedFootprints=growFootprints,
779 maskNameList=['SAT'],
780 useLegacyInterp=useLegacyInterp)
781 numBit = ipIsr.countMaskedPixels(interpMaskedImage, "INTRP")
782 if growFootprints == 0 and not useLegacyInterp:
783 # All pixel need to be interpolated over. There is
784 # no external information to interpolate over. In
785 # the GP code in this case, it is not doing
786 # anything.
787 self.assertEqual(numBit, 0,
788 msg=f"interpolateFromMask with growFootprints={growFootprints}")
789 else:
790 self.assertEqual(numBit, 40800,
791 msg=f"interpolateFromMask with growFootprints={growFootprints}")
793 def test_saturationCorrectionInterpolate(self):
794 """Expect number of mask pixels with SAT marked to be non-zero.
795 """
796 corrMaskedImage = ipIsr.saturationCorrection(self.mi, 200, 2.0,
797 growFootprints=2, interpolate=True,
798 maskName='SAT')
799 numBit = ipIsr.countMaskedPixels(corrMaskedImage, "SAT")
800 self.assertEqual(numBit, 40800)
802 def test_saturationCorrectionNoInterpolate(self):
803 """Expect number of mask pixels with SAT marked to be non-zero.
804 """
805 corrMaskedImage = ipIsr.saturationCorrection(self.mi, 200, 2.0,
806 growFootprints=2, interpolate=False,
807 maskName='SAT')
808 numBit = ipIsr.countMaskedPixels(corrMaskedImage, "SAT")
809 self.assertEqual(numBit, 40800)
811 def test_trimToMatchCalibBBox(self):
812 """Expect bounding boxes to match.
813 """
814 darkExp = isrMock.DarkMock().run()
815 darkMi = darkExp.getMaskedImage()
817 nEdge = 2
818 darkMi = darkMi[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
819 newInput = ipIsr.trimToMatchCalibBBox(self.mi, darkMi)
821 self.assertEqual(newInput.getImage().getBBox(), darkMi.getImage().getBBox())
823 def test_darkCorrection(self):
824 """Expect round-trip application to be equal.
825 Expect RuntimeError if sizes are different.
826 """
827 darkExp = isrMock.DarkMock().run()
828 darkMi = darkExp.getMaskedImage()
830 mi = self.mi.clone()
832 # The `invert` parameter controls the direction of the
833 # application. This will apply, and un-apply the dark.
834 ipIsr.darkCorrection(self.mi, darkMi, 1.0, 1.0, trimToFit=True)
835 ipIsr.darkCorrection(self.mi, darkMi, 1.0, 1.0, trimToFit=True, invert=True)
837 self.assertMaskedImagesAlmostEqual(self.mi, mi, atol=1e-3)
839 darkMi = darkMi[1:-1, 1:-1, afwImage.LOCAL]
840 with self.assertRaises(RuntimeError):
841 ipIsr.darkCorrection(self.mi, darkMi, 1.0, 1.0, trimToFit=False)
843 def test_biasCorrection(self):
844 """Expect smaller median image value after.
845 Expect RuntimeError if sizes are different.
846 """
847 biasExp = isrMock.BiasMock().run()
848 biasMi = biasExp.getMaskedImage()
850 mi = self.mi.clone()
851 ipIsr.biasCorrection(self.mi, biasMi, trimToFit=True)
852 self.assertLess(computeImageMedianAndStd(self.mi.getImage())[0],
853 computeImageMedianAndStd(mi.getImage())[0])
855 biasMi = biasMi[1:-1, 1:-1, afwImage.LOCAL]
856 with self.assertRaises(RuntimeError):
857 ipIsr.biasCorrection(self.mi, biasMi, trimToFit=False)
859 def test_flatCorrection(self):
860 """Expect round-trip application to be equal.
861 Expect RuntimeError if sizes are different.
862 """
863 flatExp = isrMock.FlatMock().run()
864 flatMi = flatExp.getMaskedImage()
866 mi = self.mi.clone()
867 for scaling in ('USER', 'MEAN', 'MEDIAN'):
868 # The `invert` parameter controls the direction of the
869 # application. This will apply, and un-apply the flat.
870 ipIsr.flatCorrection(self.mi, flatMi, scaling, userScale=1.0, trimToFit=True)
871 ipIsr.flatCorrection(self.mi, flatMi, scaling, userScale=1.0,
872 trimToFit=True, invert=True)
874 self.assertMaskedImagesAlmostEqual(self.mi, mi, atol=1e-3,
875 msg=f"flatCorrection with scaling {scaling}")
877 flatMi = flatMi[1:-1, 1:-1, afwImage.LOCAL]
878 with self.assertRaises(RuntimeError):
879 ipIsr.flatCorrection(self.mi, flatMi, 'USER', userScale=1.0, trimToFit=False)
881 def test_flatCorrectionUnknown(self):
882 """Raise if an unknown scaling is used.
884 The `scaling` parameter must be a known type. If not, the
885 flat correction will raise a RuntimeError.
886 """
887 flatExp = isrMock.FlatMock().run()
888 flatMi = flatExp.getMaskedImage()
890 with self.assertRaises(RuntimeError):
891 ipIsr.flatCorrection(self.mi, flatMi, "UNKNOWN", userScale=1.0, trimToFit=True)
893 def test_illumCorrection(self):
894 """Expect larger median value after.
895 Expect RuntimeError if sizes are different.
896 """
897 flatExp = isrMock.FlatMock().run()
898 flatMi = flatExp.getMaskedImage()
900 mi = self.mi.clone()
901 ipIsr.illuminationCorrection(self.mi, flatMi, 1.0)
902 self.assertGreater(computeImageMedianAndStd(self.mi.getImage())[0],
903 computeImageMedianAndStd(mi.getImage())[0])
905 flatMi = flatMi[1:-1, 1:-1, afwImage.LOCAL]
906 with self.assertRaises(RuntimeError):
907 ipIsr.illuminationCorrection(self.mi, flatMi, 1.0, trimToFit=False)
909 def test_brighterFatterCorrection(self):
910 """Expect smoother image/smaller std before.
911 """
912 bfKern = isrMock.BfKernelMock().run()
914 before = computeImageMedianAndStd(self.inputExp.getImage())
915 ipIsr.brighterFatterCorrection(self.inputExp, bfKern, 10, 1e-2, False)
916 after = computeImageMedianAndStd(self.inputExp.getImage())
918 self.assertLess(before[1], after[1])
920 def test_gainContext(self):
921 """Expect image to be unmodified before and after
922 """
923 mi = self.inputExp.getMaskedImage().clone()
924 with ipIsr.gainContext(self.inputExp, self.inputExp.getImage(), apply=True):
925 pass
927 self.assertIsNotNone(mi)
928 self.assertMaskedImagesAlmostEqual(self.inputExp.getMaskedImage(), mi)
930 def test_widenSaturationTrails(self):
931 """Expect more mask pixels with SAT set after.
932 """
933 numBitBefore = ipIsr.countMaskedPixels(self.mi, "SAT")
935 ipIsr.widenSaturationTrails(self.mi.getMask())
936 numBitAfter = ipIsr.countMaskedPixels(self.mi, "SAT")
938 self.assertGreaterEqual(numBitAfter, numBitBefore)
940 def test_setBadRegions(self):
941 """Expect RuntimeError if improper statistic given.
942 Expect a float value otherwise.
943 """
944 for badStatistic in ('MEDIAN', 'MEANCLIP', 'UNKNOWN'):
945 if badStatistic == 'UNKNOWN':
946 with self.assertRaises(RuntimeError,
947 msg=f"setBadRegions did not fail for stat {badStatistic}"):
948 nBad, value = ipIsr.setBadRegions(self.inputExp, badStatistic=badStatistic)
949 else:
950 nBad, value = ipIsr.setBadRegions(self.inputExp, badStatistic=badStatistic)
951 self.assertGreaterEqual(abs(value), 0.0,
952 msg=f"setBadRegions did not find valid value for stat {badStatistic}")
954 def test_attachTransmissionCurve(self):
955 """Expect no failure and non-None output from attachTransmissionCurve.
956 """
957 curve = isrMock.TransmissionMock().run()
958 combined = ipIsr.attachTransmissionCurve(self.inputExp,
959 opticsTransmission=curve,
960 filterTransmission=curve,
961 sensorTransmission=curve,
962 atmosphereTransmission=curve)
963 # DM-19707: ip_isr functionality not fully tested by unit tests
964 self.assertIsNotNone(combined)
966 def test_attachTransmissionCurve_None(self):
967 """Expect no failure and non-None output from attachTransmissionCurve.
968 """
969 curve = None
970 combined = ipIsr.attachTransmissionCurve(self.inputExp,
971 opticsTransmission=curve,
972 filterTransmission=curve,
973 sensorTransmission=curve,
974 atmosphereTransmission=curve)
975 # DM-19707: ip_isr functionality not fully tested by unit tests
976 self.assertIsNotNone(combined)
978 def test_countMaskedPixels(self):
979 mockImageConfig = isrMock.IsrMock.ConfigClass()
981 # flatDrop is not really relevant as we replace the data
982 # but good to note it in case we change how this image is made
983 mockImageConfig.flatDrop = 0.99999
984 mockImageConfig.isTrimmed = True
986 flatExp = isrMock.FlatMock(config=mockImageConfig).run()
987 (shapeY, shapeX) = flatExp.getDimensions()
989 rng = np.random.RandomState(0)
990 flatMean = 1000
991 flatWidth = np.sqrt(flatMean)
992 flatData = rng.normal(flatMean, flatWidth, (shapeX, shapeY))
993 flatExp.image.array[:] = flatData
995 exp = flatExp.clone()
996 mi = exp.maskedImage
997 self.assertEqual(ipIsr.countMaskedPixels(mi, 'NO_DATA'), 0)
998 self.assertEqual(ipIsr.countMaskedPixels(mi, 'BAD'), 0)
1000 NODATABIT = mi.mask.getPlaneBitMask("NO_DATA")
1001 noDataBox = geom.Box2I(geom.Point2I(31, 49), geom.Extent2I(3, 6))
1002 mi.mask[noDataBox] |= NODATABIT
1004 self.assertEqual(ipIsr.countMaskedPixels(mi, 'NO_DATA'), noDataBox.getArea())
1005 self.assertEqual(ipIsr.countMaskedPixels(mi, 'BAD'), 0)
1007 mi.mask[noDataBox] ^= NODATABIT # XOR to reset what we did
1008 self.assertEqual(ipIsr.countMaskedPixels(mi, 'NO_DATA'), 0)
1010 def test_getExposureGains(self):
1011 exposure = self.inputExp.clone()
1012 metadata = exposure.metadata
1013 amps = exposure.getDetector().getAmplifiers()
1015 # Check the default values
1016 values = ipIsr.getExposureGains(exposure).values()
1017 self.assertEqual(list(values), len(amps)*[1.0])
1019 # Set values using old ISR keys
1020 for amp in amps:
1021 metadata[f"LSST GAIN {amp.getName()}"] = 2.0
1022 values = ipIsr.getExposureGains(exposure).values()
1023 self.assertEqual(list(values), len(amps)*[2.0])
1025 # Set values using new ISR keys
1026 for amp in amps:
1027 metadata[f"LSST ISR GAIN {amp.getName()}"] = 3.0
1028 values = ipIsr.getExposureGains(exposure).values()
1029 self.assertEqual(list(values), len(amps)*[3.0])
1031 def test_getExposureReadNoises(self):
1032 exposure = self.inputExp.clone()
1033 metadata = exposure.metadata
1034 amps = exposure.getDetector().getAmplifiers()
1036 # Check the default values
1037 values = ipIsr.getExposureReadNoises(exposure).values()
1038 self.assertEqual(list(values), len(amps)*[10.0])
1040 # Set values using old ISR keys
1041 for amp in amps:
1042 metadata[f"LSST READNOISE {amp.getName()}"] = 2.0
1043 values = ipIsr.getExposureReadNoises(exposure).values()
1044 self.assertEqual(list(values), len(amps)*[2.0])
1046 # Set values using new ISR keys
1047 for amp in amps:
1048 metadata[f"LSST ISR READNOISE {amp.getName()}"] = 3.0
1049 values = ipIsr.getExposureReadNoises(exposure).values()
1050 self.assertEqual(list(values), len(amps)*[3.0])
1052 def test_compareCameraKeywords(self):
1053 keywords = ["SEQFILE", "SEQNAME", "VOLTAGE"]
1055 # Test success, with with warnings and doRaise.
1056 exposureMetadata = PropertyList()
1057 exposureMetadata["SEQFILE"] = "filename"
1058 exposureMetadata["SEQNAME"] = "v30"
1059 exposureMetadata["VOLTAGE"] = 100.0
1061 calib = ipIsr.IsrCalib()
1062 calib._metadata["SEQFILE"] = "filename"
1063 calib._metadata["SEQNAME"] = "v30"
1064 calib._metadata["VOLTAGE"] = 100.0
1066 with self.assertNoLogs(level=logging.WARNING):
1067 ipIsrFunctions.compareCameraKeywords(False, keywords, exposureMetadata, calib, "test")
1069 ipIsrFunctions.compareCameraKeywords(True, keywords, exposureMetadata, calib, "test")
1071 # Test warning on missing exposure keyword.
1072 exposureMetadata["VOLTAGE"] = None
1073 with self.assertLogs(level=logging.DEBUG) as cm:
1074 ipIsrFunctions.compareCameraKeywords(False, keywords, exposureMetadata, calib, "test")
1075 self.assertEqual(len(cm[1]), 1)
1076 self.assertIn("Sequencer keyword VOLTAGE not found", cm[1][0])
1078 # Test no warning + info on missing calib keywords.
1079 exposureMetadata["VOLTAGE"] = 100.0
1080 calib._metadata["SEQFILE"] = None
1081 calib._metadata["SEQNAME"] = None
1082 with self.assertNoLogs(level=logging.WARNING):
1083 ipIsrFunctions.compareCameraKeywords(False, keywords, exposureMetadata, calib, "test")
1084 with self.assertLogs(level=logging.INFO) as cm:
1085 ipIsrFunctions.compareCameraKeywords(False, keywords, exposureMetadata, calib, "test")
1086 self.assertEqual(len(cm[1]), 1)
1087 self.assertIn("Calibration test missing keywords SEQFILE,SEQNAME", cm[1][0])
1089 # Test mismatch, with warnings and doRaise.
1090 calib._metadata["VOLTAGE"] = 101.0
1091 with self.assertLogs(level=logging.WARNING) as cm:
1092 ipIsrFunctions.compareCameraKeywords(False, keywords, exposureMetadata, calib, "test")
1093 self.assertEqual(len(cm[1]), 1)
1094 self.assertIn("Sequencer mismatch", cm[1][0])
1095 self.assertIn("ISR TEST SEQUENCER MISMATCH", exposureMetadata)
1096 self.assertEqual(exposureMetadata["ISR TEST SEQUENCER MISMATCH"], True)
1098 with self.assertRaises(RuntimeError):
1099 ipIsrFunctions.compareCameraKeywords(True, keywords, exposureMetadata, calib, "test")
1102class MemoryTester(lsst.utils.tests.MemoryTestCase):
1103 pass
1106def setup_module(module):
1107 lsst.utils.tests.init()
1110if __name__ == "__main__": 1110 ↛ 1111line 1110 didn't jump to line 1111 because the condition on line 1110 was never true
1111 lsst.utils.tests.init()
1112 unittest.main()