Coverage for python/lsst/ip/isr/isrFunctions.py: 87%
548 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 22:35 +0000
1#
2# LSST Data Management System
3# Copyright 2008, 2009, 2010 LSST Corporation.
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 <http://www.lsstcorp.org/LegalNotices/>.
21#
23__all__ = [
24 "applyGains",
25 "attachTransmissionCurve",
26 "biasCorrection",
27 "checkFilter",
28 "compareCameraKeywords",
29 "countMaskedPixels",
30 "createPsf",
31 "darkCorrection",
32 "flatCorrection",
33 "gainContext",
34 "getPhysicalFilter",
35 "growMasks",
36 "maskDECamEdgeBleed",
37 "maskE2VEdgeBleed",
38 "maskITLEdgeBleed",
39 "maskITLSatSag",
40 "maskITLDip",
41 "illuminationCorrection",
42 "interpolateDefectList",
43 "interpolateFromMask",
44 "makeThresholdMask",
45 "saturationCorrection",
46 "setBadRegions",
47 "transposeMaskedImage",
48 "trimToMatchCalibBBox",
49 "updateVariance",
50 "widenSaturationTrails",
51 "getExposureGains",
52 "getExposureReadNoises",
53]
55import logging
56import math
57import numpy
59import lsst.geom
60import lsst.afw.image as afwImage
61import lsst.afw.detection as afwDetection
62import lsst.afw.math as afwMath
63import lsst.meas.algorithms as measAlg
64import lsst.afw.cameraGeom as camGeom
66from lsst.afw.geom import SpanSet, Stencil
67from lsst.meas.algorithms.detection import SourceDetectionTask
69from contextlib import contextmanager
71from .defects import Defects
74def createPsf(fwhm):
75 """Make a double Gaussian PSF.
77 Parameters
78 ----------
79 fwhm : scalar
80 FWHM of double Gaussian smoothing kernel.
82 Returns
83 -------
84 psf : `lsst.meas.algorithms.DoubleGaussianPsf`
85 The created smoothing kernel.
86 """
87 ksize = 4*int(fwhm) + 1
88 return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
91def transposeMaskedImage(maskedImage):
92 """Make a transposed copy of a masked image.
94 Parameters
95 ----------
96 maskedImage : `lsst.afw.image.MaskedImage`
97 Image to process.
99 Returns
100 -------
101 transposed : `lsst.afw.image.MaskedImage`
102 The transposed copy of the input image.
103 """
104 transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
105 transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
106 transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
107 transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
108 return transposed
111def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None,
112 maskNameList=None, useLegacyInterp=True):
113 """Interpolate over defects specified in a defect list.
115 Parameters
116 ----------
117 maskedImage : `lsst.afw.image.MaskedImage`
118 Image to process.
119 defectList : `lsst.meas.algorithms.Defects`
120 List of defects to interpolate over.
121 fwhm : `float`
122 FWHM of double Gaussian smoothing kernel.
123 fallbackValue : scalar, optional
124 Fallback value if an interpolated value cannot be determined.
125 If None, then the clipped mean of the image is used.
126 maskNameList : `list [string]`
127 List of the defects to interpolate over (used for GP interpolator).
128 useLegacyInterp : `bool`
129 Use the legacy interpolation (polynomial interpolation) if True. Use
130 Gaussian Process interpolation if False.
132 Notes
133 -----
134 The ``fwhm`` parameter is used to create a PSF, but the underlying
135 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
136 not currently make use of this information in legacy Interpolation, but use
137 if for the Gaussian Process as an estimation of the correlation lenght.
138 """
139 psf = createPsf(fwhm)
140 if fallbackValue is None: 140 ↛ 142line 140 didn't jump to line 142 because the condition on line 140 was always true
141 fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
142 if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict(): 142 ↛ 143line 142 didn't jump to line 143 because the condition on line 142 was never true
143 maskedImage.getMask().addMaskPlane('INTRP')
145 # Hardcoded fwhm value. PSF estimated latter in step1,
146 # not in ISR.
147 if useLegacyInterp:
148 kwargs = {}
149 fwhm = fwhm
150 else:
151 # tested on a dozens of images and looks a good set of
152 # hyperparameters, but cannot guarrenty this is optimal,
153 # need further testing.
154 kwargs = {"bin_spacing": 20,
155 "threshold_dynamic_binning": 2000,
156 "threshold_subdivide": 20000}
157 fwhm = 15
159 measAlg.interpolateOverDefects(maskedImage, psf, defectList,
160 fallbackValue=fallbackValue,
161 useFallbackValueAtEdge=True,
162 fwhm=fwhm,
163 useLegacyInterp=useLegacyInterp,
164 maskNameList=maskNameList, **kwargs)
165 return maskedImage
168def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
169 """Mask pixels based on threshold detection.
171 Parameters
172 ----------
173 maskedImage : `lsst.afw.image.MaskedImage`
174 Image to process. Only the mask plane is updated.
175 threshold : scalar
176 Detection threshold.
177 growFootprints : scalar, optional
178 Number of pixels to grow footprints of detected regions.
179 maskName : str, optional
180 Mask plane name, or list of names to convert
182 Returns
183 -------
184 defectList : `lsst.meas.algorithms.Defects`
185 Defect list constructed from pixels above the threshold.
186 """
187 # find saturated regions
188 thresh = afwDetection.Threshold(threshold)
189 fs = afwDetection.FootprintSet(maskedImage, thresh)
191 if growFootprints > 0:
192 fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
193 fpList = fs.getFootprints()
195 # set mask
196 mask = maskedImage.getMask()
197 bitmask = mask.getPlaneBitMask(maskName)
198 afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
200 return Defects.fromFootprintList(fpList)
203def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
204 """Grow a mask by an amount and add to the requested plane.
206 Parameters
207 ----------
208 mask : `lsst.afw.image.Mask`
209 Mask image to process.
210 radius : scalar
211 Amount to grow the mask.
212 maskNameList : `str` or `list` [`str`]
213 Mask names that should be grown.
214 maskValue : `str`
215 Mask plane to assign the newly masked pixels to.
216 """
217 if radius > 0: 217 ↛ exitline 217 didn't return from function 'growMasks' because the condition on line 217 was always true
218 spans = SpanSet.fromMask(mask, mask.getPlaneBitMask(maskNameList))
219 # Use MANHATTAN for equivalence with 'isotropic=False` footprint grows,
220 # but CIRCLE is probably better and might be just as fast.
221 spans = spans.dilated(radius, Stencil.MANHATTAN)
222 spans = spans.clippedTo(mask.getBBox())
223 spans.setMask(mask, mask.getPlaneBitMask(maskValue))
226def maskE2VEdgeBleed(exposure, e2vEdgeBleedSatMinArea=10000,
227 e2vEdgeBleedSatMaxArea=100000,
228 e2vEdgeBleedYMax=350,
229 saturatedMaskName="SAT", log=None):
230 """Mask edge bleeds in E2V detectors.
232 Parameters
233 ----------
234 exposure : `lsst.afw.image.Exposure`
235 Exposure to apply masking to.
236 e2vEdgeBleedSatMinArea : `int`, optional
237 Minimum limit of saturated cores footprint area.
238 e2vEdgeBleedSatMaxArea : `int`, optional
239 Maximum limit of saturated cores footprint area.
240 e2vEdgeBleedYMax: `float`, optional
241 Height of edge bleed masking.
242 saturatedMaskName : `str`, optional
243 Mask name for saturation.
244 log : `logging.Logger`, optional
245 Logger to handle messages.
246 """
248 log = log if log else logging.getLogger(__name__)
250 maskedImage = exposure.maskedImage
251 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
253 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
255 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
257 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
258 largeAreas, = numpy.where((satAreas >= e2vEdgeBleedSatMinArea)
259 & (satAreas < e2vEdgeBleedSatMaxArea))
260 for largeAreasIndex in largeAreas:
261 fpCore = fpList[largeAreasIndex]
262 xCore, yCore = fpCore.getCentroid()
263 xCore = int(xCore)
264 yCore = int(yCore)
266 for amp in exposure.getDetector():
267 if amp.getBBox().contains(xCore, yCore):
268 ampName = amp.getName()
269 if ampName[:2] == 'C0': 269 ↛ 266line 269 didn't jump to line 266 because the condition on line 269 was always true
270 # Check that the footprint reaches the bottom of the
271 # amplifier.
272 if fpCore.getBBox().getMinY() == 0: 272 ↛ 266line 272 didn't jump to line 266 because the condition on line 272 was always true
273 # This is a large saturation footprint that hits the
274 # edge, and is thus classified as an edge bleed.
276 # TODO DM-50587: Optimize number of rows to mask by
277 # looking at the median signal level as a function of
278 # row number on the right side of the saturation trail.
280 log.info("Found E2V edge bleed in amp %s, column %d.", ampName, xCore)
281 maskedImage.mask[amp.getBBox()].array[:e2vEdgeBleedYMax, :] |= saturatedBit
284def maskDECamEdgeBleed(exposure, satMinArea=10000, satMaxArea=100000, approachRows=20,
285 nSigma=5.0, nRowsCheck=20, minLowPixelsPerRow=30, minLowPixelsExtent=10,
286 marginFraction=0.125, saturatedMaskName="SAT", log=None):
287 """Mask DECam-style edge bleeds with low rows next to the read register.
289 For each amplifier that contains a large saturated footprint reaching
290 within ``approachRows`` of its read edge, this function confirms the dip
291 from the per-row count of pixels more than ``nSigma`` below the clipped
292 sky, measures how far the dip extends inward, and marks those rows
293 (plus a margin) with the saturation mask plane.
295 Parameters
296 ----------
297 exposure : `lsst.afw.image.Exposure`
298 Assembled exposure to mask.
299 satMinArea : `int`, optional
300 Minimum area (pixels) of a saturated footprint to be considered.
301 satMaxArea : `int`, optional
302 Maximum area (pixels) of a saturated footprint to be considered.
303 approachRows : `int`, optional
304 A footprint must come within this many rows of the read edge.
305 nSigma : `float`, optional
306 A pixel is "low" if it is more than this many sigma below sky.
307 nRowsCheck : `int`, optional
308 Number of rows from the read edge used to confirm the bleed.
309 minLowPixelsPerRow : `int`, optional
310 Mean number of low pixels per row over the check rows required
311 to confirm the bleed.
312 minLowPixelsExtent : `int`, optional
313 Number of low pixels a row must exceed to count toward the bleed
314 height; the scan stops after five consecutive rows that do not.
315 Should be smaller than ``minLowPixelsPerRow``.
316 marginFraction : `float`, optional
317 Extra rows masked beyond the measured height, as a fraction of
318 that height (plus one row).
319 saturatedMaskName : `str`, optional
320 Mask plane name for saturation.
321 log : `logging.Logger`, optional
322 Logger to handle messages.
324 Raises
325 ------
326 RuntimeError
327 Raised if the exposure has no detector with amplifiers.
328 """
329 log = log if log else logging.getLogger(__name__)
331 detector = exposure.getDetector()
332 if detector is None or len(detector) == 0:
333 raise RuntimeError("Cannot mask DECam edge bleeds: exposure has no detector with amplifiers.")
335 maskedImage = exposure.maskedImage
336 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
337 unusableBits = maskedImage.mask.getPlaneBitMask([saturatedMaskName, "BAD", "NO_DATA"])
338 statsBits = unusableBits | maskedImage.mask.getPlaneBitMask("SUSPECT")
340 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
341 fpList = afwDetection.FootprintSet(maskedImage.mask, thresh).getFootprints()
342 candidates = []
343 for fp in fpList:
344 if satMinArea <= fp.getArea() < satMaxArea:
345 xCore, yCore = fp.getCentroid()
346 candidates.append((int(xCore), int(yCore), fp.getBBox()))
347 if not candidates:
348 log.debug("No saturated footprints in the DECam edge bleed area range.")
349 return
351 statsCtrl = afwMath.StatisticsControl()
352 statsCtrl.setAndMask(statsBits)
353 # Number of consecutive rows without a dip that ends the height scan.
354 nRowsStop = 5
356 for amp in detector:
357 ampBBox = amp.getBBox()
358 if not exposure.getBBox().contains(ampBBox):
359 continue
360 readTop = amp.getReadoutCorner() in (camGeom.ReadoutCorner.UL, camGeom.ReadoutCorner.UR)
362 if readTop:
363 reachesReadEdge = any(ampBBox.contains(x, y)
364 and fpBBox.getMaxY() >= ampBBox.getMaxY() - approachRows
365 for x, y, fpBBox in candidates)
366 else:
367 reachesReadEdge = any(ampBBox.contains(x, y)
368 and fpBBox.getMinY() <= ampBBox.getMinY() + approachRows
369 for x, y, fpBBox in candidates)
370 if not reachesReadEdge:
371 continue
373 ampImage = maskedImage[ampBBox]
374 stats = afwMath.makeStatistics(ampImage, afwMath.MEANCLIP | afwMath.STDEVCLIP | afwMath.NPOINT,
375 statsCtrl)
376 if stats.getValue(afwMath.NPOINT) < 1000: 376 ↛ 377line 376 didn't jump to line 377 because the condition on line 376 was never true
377 log.debug("Skipping DECam edge bleed check in amp %s: too few unmasked pixels.",
378 amp.getName())
379 continue
380 sky = stats.getValue(afwMath.MEANCLIP)
381 sigma = stats.getValue(afwMath.STDEVCLIP)
383 usable = (ampImage.mask.array & unusableBits) == 0
384 low = (ampImage.image.array < sky - nSigma*sigma) & usable
385 # Row profiles with index 0 at the read edge, increasing inward.
386 profile = low.sum(axis=1)
387 nUsable = usable.sum(axis=1)
388 if readTop:
389 profile = profile[::-1]
390 nUsable = nUsable[::-1]
392 # Skip edge rows that cannot hold enough low pixels to confirm a dip.
393 enoughUsable = numpy.nonzero(nUsable >= minLowPixelsPerRow)[0]
394 if len(enoughUsable) == 0: 394 ↛ 395line 394 didn't jump to line 395 because the condition on line 394 was never true
395 log.debug("Skipping DECam edge bleed check in amp %s: no rows with usable pixels.",
396 amp.getName())
397 continue
398 start = enoughUsable[0]
400 if profile[start:start + nRowsCheck].mean() <= minLowPixelsPerRow:
401 log.debug("Saturated footprint reaches read edge of amp %s but no dip found.",
402 amp.getName())
403 continue
405 height = 0
406 nBelow = 0
407 for i in range(start, len(profile)): 407 ↛ 415line 407 didn't jump to line 415 because the loop on line 407 didn't complete
408 if profile[i] > minLowPixelsExtent:
409 height = i + 1
410 nBelow = 0
411 else:
412 nBelow += 1
413 if nBelow >= nRowsStop:
414 break
415 height += int(height*marginFraction) + 1
416 height = min(height, ampBBox.getHeight())
418 if readTop:
419 maskedImage.mask[ampBBox].array[-height:, :] |= saturatedBit
420 else:
421 maskedImage.mask[ampBBox].array[:height, :] |= saturatedBit
422 log.info("Found DECam edge bleed in amp %s at the %s read edge; masked %d rows.",
423 amp.getName(), "top" if readTop else "bottom", height)
426def maskITLEdgeBleed(ccdExposure, badAmpDict,
427 fpCore, itlEdgeBleedSatMinArea=10000,
428 itlEdgeBleedSatMaxArea=100000,
429 itlEdgeBleedThreshold=5000.,
430 itlEdgeBleedModelConstant=0.02,
431 saturatedMaskName="SAT", log=None):
432 """Mask edge bleeds in ITL detectors.
434 Parameters
435 ----------
436 ccdExposure : `lsst.afw.image.Exposure`
437 Exposure to apply masking to.
438 badAmpDict : `dict` [`str`, `bool`]
439 Dictionary of amplifiers, keyed by name, value is True if
440 amplifier is fully masked.
441 fpCore : `lsst.afw.detection._detection.Footprint`
442 Footprint of saturated core.
443 itlEdgeBleedThreshold : `float`, optional
444 Threshold above median sky background for edge bleed detection
445 (electron units).
446 itlEdgeBleedModelConstant : `float`, optional
447 Constant in the decaying exponential in the edge bleed masking.
448 saturatedMaskName : `str`, optional
449 Mask name for saturation.
450 log : `logging.Logger`, optional
451 Logger to handle messages.
452 """
454 log = log if log else logging.getLogger(__name__)
456 # Get median of amplifier saturation level
457 satLevel = numpy.nanmedian([ccdExposure.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
458 for amp in ccdExposure.getDetector() if not badAmpDict[amp.getName()]])
460 # 1. we check if there are several cores in the footprint:
461 # Get centroid of saturated core
462 xCore, yCore = fpCore.getCentroid()
463 # Turn the Y detector coordinate into Y footprint coordinate
464 yCoreFP = int(yCore) - fpCore.getBBox().getMinY()
465 # Now test if there is one or more cores by checking if the slice at the
466 # center is full of saturated pixels or has several segments of saturated
467 # columns (i.e. several cores with trails)
468 checkCoreNbRow = fpCore.getSpans().asArray()[yCoreFP, :]
469 nbCore = 0
470 indexSwitchTrue = []
471 indexSwitchFalse = []
472 if checkCoreNbRow[0]:
473 # If the slice starts with saturated pixels
474 inSatSegment = True
475 nbCore = 1
476 indexSwitchTrue.append(0)
477 else:
478 # If the slice starts with non saturated pixels
479 inSatSegment = False
481 for i, value in enumerate(checkCoreNbRow):
482 if value:
483 if not inSatSegment:
484 indexSwitchTrue.append(i)
485 # nbCore is the number of detected cores.
486 nbCore += 1
487 inSatSegment = True
488 elif inSatSegment:
489 indexSwitchFalse.append(i)
490 inSatSegment = False
492 # 1. we look for edge bleed in saturated cores in the footprint
493 if nbCore == 2:
494 # we now estimate the x coordinates of the edges of the subfootprint
495 # for each core
496 xEdgesCores = [0]
497 xEdgesCores.append(int((indexSwitchTrue[1] + indexSwitchFalse[0])/2))
498 xEdgesCores.append(fpCore.getSpans().asArray().shape[1])
499 # Get the X and Y footprint coordinates of the cores
500 for i in range(nbCore):
501 subfp = fpCore.getSpans().asArray()[:, xEdgesCores[i]:xEdgesCores[i+1]]
502 xCoreFP = int(xEdgesCores[i] + numpy.argmax(numpy.sum(subfp, axis=0)))
503 # turn into X coordinate in detector space
504 xCore = xCoreFP + fpCore.getBBox().getMinX()
505 # get Y footprint coordinate of the core
506 # by trimming the edges where edge bleeds are potentially dominant
507 if subfp.shape[0] <= 200: 507 ↛ 510line 507 didn't jump to line 510 because the condition on line 507 was always true
508 yCoreFP = int(numpy.argmax(numpy.sum(subfp, axis=1)))
509 else:
510 yCoreFP = int(numpy.argmax(numpy.sum(subfp[100:-100, :],
511 axis=1)))
512 yCoreFP = 100+yCoreFP
514 # Estimate the width of the saturated core
515 widthSat = numpy.sum(subfp[int(yCoreFP), :])
517 subfpArea = numpy.sum(subfp)
518 if subfpArea > itlEdgeBleedSatMinArea and subfpArea < itlEdgeBleedSatMaxArea: 518 ↛ 500line 518 didn't jump to line 500 because the condition on line 518 was always true
519 _applyMaskITLEdgeBleed(ccdExposure, xCore,
520 satLevel, widthSat,
521 itlEdgeBleedThreshold,
522 itlEdgeBleedModelConstant,
523 saturatedMaskName, log)
524 elif nbCore > 2: 524 ↛ 526line 524 didn't jump to line 526 because the condition on line 524 was never true
525 # TODO DM-49736: support N cores in saturated footprint
526 log.warning(
527 "Too many (%d) cores in saturated footprint to mask edge bleeds.",
528 nbCore,
529 )
530 else:
531 # Get centroid of saturated core
532 xCore, yCore = fpCore.getCentroid()
533 # Turn the Y detector coordinate into Y footprint coordinate
534 yCoreFP = yCore - fpCore.getBBox().getMinY()
535 # Get the number of saturated columns around the centroid
536 widthSat = numpy.sum(fpCore.getSpans().asArray()[int(yCoreFP), :])
537 _applyMaskITLEdgeBleed(ccdExposure, xCore,
538 satLevel, widthSat, itlEdgeBleedThreshold,
539 itlEdgeBleedModelConstant, saturatedMaskName, log)
542def _applyMaskITLEdgeBleed(ccdExposure, xCore,
543 satLevel, widthSat,
544 itlEdgeBleedThreshold=5000.,
545 itlEdgeBleedModelConstant=0.03,
546 saturatedMaskName="SAT", log=None):
547 """Apply ITL edge bleed masking model.
549 Parameters
550 ----------
551 ccdExposure : `lsst.afw.image.Exposure`
552 Exposure to apply masking to.
553 xCore: `int`
554 X coordinate of the saturated core.
555 satLevel: `float`
556 Minimum saturation level of the detector.
557 widthSat: `float`
558 Width of the saturated core.
559 itlEdgeBleedThreshold : `float`, optional
560 Threshold above median sky background for edge bleed detection
561 (electron units).
562 itlEdgeBleedModelConstant : `float`, optional
563 Constant in the decaying exponential in the edge bleed masking.
564 saturatedMaskName : `str`, optional
565 Mask name for saturation.
566 log : `logging.Logger`, optional
567 Logger to handle messages.
568 """
569 log = log if log else logging.getLogger(__name__)
571 maskedImage = ccdExposure.maskedImage
572 xmax = maskedImage.image.array.shape[1]
573 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
575 for amp in ccdExposure.getDetector():
576 # Select the 2 top and bottom amplifiers around the saturated
577 # core with a potential edge bleed by selecting the amplifiers
578 # that have the same X coordinate as the saturated core.
579 # As we don't care about the Y coordinate, we set it to the
580 # center of the BBox.
581 yBox = amp.getBBox().getCenter()[1]
582 if amp.getBBox().contains(xCore, yBox):
584 # Get the amp name
585 ampName = amp.getName()
587 # Because in ITLs the edge bleed happens on both edges
588 # of the detector, we make a cutout around
589 # both the top and bottom
590 # edge bleed candidates around the saturated core.
591 # We flip the cutout of the top amplifier
592 # to then work with the same coordinates for both.
593 # The way of selecting top vs bottom amp
594 # is very specific to ITL.
595 if ampName[:2] == 'C1':
596 sliceImage = maskedImage.image.array[:200, :]
597 sliceMask = maskedImage.mask.array[:200, :]
598 elif ampName[:2] == 'C0': 598 ↛ 612line 598 didn't jump to line 612 because the condition on line 598 was always true
599 sliceImage = numpy.flipud(maskedImage.image.array[-200:, :])
600 sliceMask = numpy.flipud(maskedImage.mask.array[-200:, :])
602 # The middle columns of edge bleeds often have
603 # high counts, so we check there is an edge bleed
604 # by looking at a small image up to 50 pixels from the edge
605 # and around the saturated columns
606 # of the saturated core, and checking its median is
607 # above the sky background by itlEdgeBleedThreshold
609 # If the centroid is too close to the edge of the detector
610 # (within 5 pixels), we set the limit to the mean check
611 # to the edge of the detector
612 lowerRangeSmall = int(xCore)-5
613 upperRangeSmall = int(xCore)+5
614 if lowerRangeSmall < 0:
615 lowerRangeSmall = 0
616 if upperRangeSmall > xmax: 616 ↛ 617line 616 didn't jump to line 617 because the condition on line 616 was never true
617 upperRangeSmall = xmax
618 ampImageBG = numpy.median(maskedImage[amp.getBBox()].image.array)
619 edgeMedian = numpy.median(sliceImage[:50, lowerRangeSmall:upperRangeSmall])
620 if edgeMedian > (ampImageBG + itlEdgeBleedThreshold):
622 log.info("Found ITL edge bleed in amp %s, column %d.", ampName, xCore)
624 # We need an estimate of the maximum width
625 # of the edge bleed for our masking model
626 # so we now estimate it by measuring the width of
627 # areas above 60 percent of the saturation level
628 # close to the edge,
629 # in a cutout up to 100 pixels from the edge,
630 # with a width of around the width of an amplifier.
631 subImageXMin = int(xCore)-250
632 subImageXMax = int(xCore)+250
633 if subImageXMin < 0:
634 subImageXMin = 0
635 elif subImageXMax > xmax:
636 subImageXMax = xmax
638 subImage = sliceImage[:100, subImageXMin:subImageXMax]
639 maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
640 axis=1))
642 # Mask edge bleed with a decaying exponential model
643 for y in range(200):
644 edgeBleedHalfWidth = \
645 int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
646 + widthSat)/2.)
647 lowerRange = int(xCore)-edgeBleedHalfWidth
648 upperRange = int(xCore)+edgeBleedHalfWidth
649 # If the edge bleed model goes outside the detector
650 # we set the limit for the masking
651 # to the edge of the detector
652 if lowerRange < 0:
653 lowerRange = 0
654 if upperRange > xmax:
655 upperRange = xmax
656 sliceMask[y, lowerRange:upperRange] |= saturatedBit
659def maskITLSatSag(ccdExposure, fpCore, saturatedMaskName="SAT"):
660 """Mask columns presenting saturation sag in saturated footprints in
661 ITL detectors.
663 Parameters
664 ----------
665 ccdExposure : `lsst.afw.image.Exposure`
666 Exposure to apply masking to.
667 fpCore : `lsst.afw.detection._detection.Footprint`
668 Footprint of saturated core.
669 saturatedMaskName : `str`, optional
670 Mask name for saturation.
671 """
673 # TODO DM-49736: add a flux level check to apply masking
675 maskedImage = ccdExposure.maskedImage
676 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
678 cc = numpy.sum(fpCore.getSpans().asArray(), axis=0)
679 # Mask full columns that have 20 percent of the height of the footprint
680 # saturated
681 columnsToMaskFP = numpy.where(cc > fpCore.getSpans().asArray().shape[0]/5.)
683 columnsToMask = [x + int(fpCore.getBBox().getMinX()) for x in columnsToMaskFP]
684 maskedImage.mask.array[:, columnsToMask] |= saturatedBit
687def maskITLDip(exposure, detectorConfig, maskPlaneNames=["SUSPECT", "ITL_DIP"], log=None):
688 """Add mask bits according to the ITL dip model.
690 Parameters
691 ----------
692 exposure : `lsst.afw.image.Exposure`
693 Exposure to do ITL dip masking.
694 detectorConfig : `lsst.ip.isr.overscanAmpConfig.OverscanDetectorConfig`
695 Configuration for this detector.
696 maskPlaneNames : `list [`str`], optional
697 Name of the ITL Dip mask planes.
698 log : `logging.Logger`, optional
699 If not set, a default logger will be used.
700 """
701 if detectorConfig.itlDipBackgroundFraction == 0.0:
702 # Nothing to do.
703 return
705 if log is None: 705 ↛ 708line 705 didn't jump to line 708 because the condition on line 705 was always true
706 log = logging.getLogger(__name__)
708 thresh = afwDetection.Threshold(
709 exposure.mask.getPlaneBitMask("SAT"),
710 afwDetection.Threshold.BITMASK,
711 )
712 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
714 heights = numpy.asarray([fp.getBBox().getHeight() for fp in fpList])
716 largeHeights, = numpy.where(heights >= detectorConfig.itlDipMinHeight)
718 if len(largeHeights) == 0: 718 ↛ 719line 718 didn't jump to line 719 because the condition on line 718 was never true
719 return
721 # Get the approximate image background.
722 approxBackground = numpy.median(exposure.image.array)
723 maskValue = exposure.mask.getPlaneBitMask(maskPlaneNames)
725 maskBak = exposure.mask.array.copy()
726 nMaskedCols = 0
728 for index in largeHeights:
729 fp = fpList[index]
730 center = fp.getCentroid()
732 nSat = numpy.sum(fp.getSpans().asArray(), axis=0)
733 width = numpy.sum(nSat > detectorConfig.itlDipMinHeight)
735 if width < detectorConfig.itlDipMinWidth:
736 continue
738 width = numpy.clip(width, None, detectorConfig.itlDipMaxWidth)
740 dipMax = detectorConfig.itlDipBackgroundFraction * approxBackground * width
742 # Assume sky-noise dominated; we could add in read noise here.
743 if dipMax < detectorConfig.itlDipMinBackgroundNoiseFraction * numpy.sqrt(approxBackground):
744 continue
746 minCol = int(center.getX() - (detectorConfig.itlDipWidthScale * width) / 2.)
747 maxCol = int(center.getX() + (detectorConfig.itlDipWidthScale * width) / 2.)
748 minCol = numpy.clip(minCol, 0, None)
749 maxCol = numpy.clip(maxCol, None, exposure.mask.array.shape[1] - 1)
751 log.info(
752 "Found ITL dip (width %d; bkg %.2f); masking column %d to %d.",
753 width,
754 approxBackground,
755 minCol,
756 maxCol,
757 )
759 exposure.mask.array[:, minCol: maxCol + 1] |= maskValue
761 nMaskedCols += (maxCol - minCol + 1)
763 if nMaskedCols > detectorConfig.itlDipMaxColsPerImage:
764 log.warning(
765 "Too many (%d) columns would be masked on this image from dip masking; restoring original mask.",
766 nMaskedCols,
767 )
768 exposure.mask.array[:, :] = maskBak
771def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
772 maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
773 """Interpolate over defects identified by a particular set of mask planes.
775 Parameters
776 ----------
777 maskedImage : `lsst.afw.image.MaskedImage`
778 Image to process.
779 fwhm : `float`
780 FWHM of double Gaussian smoothing kernel.
781 growSaturatedFootprints : scalar, optional
782 Number of pixels to grow footprints for saturated pixels.
783 maskNameList : `List` of `str`, optional
784 Mask plane name.
785 fallbackValue : scalar, optional
786 Value of last resort for interpolation.
788 Notes
789 -----
790 The ``fwhm`` parameter is used to create a PSF, but the underlying
791 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
792 not currently make use of this information.
793 """
794 mask = maskedImage.getMask()
796 if growSaturatedFootprints > 0 and "SAT" in maskNameList:
797 # If we are interpolating over an area larger than the original masked
798 # region, we need to expand the original mask bit to the full area to
799 # explain why we interpolated there.
800 growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
802 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
803 fpSet = afwDetection.FootprintSet(mask, thresh)
804 defectList = Defects.fromFootprintList(fpSet.getFootprints())
806 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
807 maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
809 return maskedImage
812def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
813 fallbackValue=None, useLegacyInterp=True):
814 """Mark saturated pixels and optionally interpolate over them
816 Parameters
817 ----------
818 maskedImage : `lsst.afw.image.MaskedImage`
819 Image to process.
820 saturation : scalar
821 Saturation level used as the detection threshold.
822 fwhm : `float`
823 FWHM of double Gaussian smoothing kernel.
824 growFootprints : scalar, optional
825 Number of pixels to grow footprints of detected regions.
826 interpolate : Bool, optional
827 If True, saturated pixels are interpolated over.
828 maskName : str, optional
829 Mask plane name.
830 fallbackValue : scalar, optional
831 Value of last resort for interpolation.
833 Notes
834 -----
835 The ``fwhm`` parameter is used to create a PSF, but the underlying
836 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
837 not currently make use of this information.
838 """
839 defectList = makeThresholdMask(
840 maskedImage=maskedImage,
841 threshold=saturation,
842 growFootprints=growFootprints,
843 maskName=maskName,
844 )
845 if interpolate:
846 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
847 maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
849 return maskedImage
852def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
853 """Compute number of edge trim pixels to match the calibration data.
855 Use the dimension difference between the raw exposure and the
856 calibration exposure to compute the edge trim pixels. This trim
857 is applied symmetrically, with the same number of pixels masked on
858 each side.
860 Parameters
861 ----------
862 rawMaskedImage : `lsst.afw.image.MaskedImage`
863 Image to trim.
864 calibMaskedImage : `lsst.afw.image.MaskedImage`
865 Calibration image to draw new bounding box from.
867 Returns
868 -------
869 replacementMaskedImage : `lsst.afw.image.MaskedImage`
870 ``rawMaskedImage`` trimmed to the appropriate size.
872 Raises
873 ------
874 RuntimeError
875 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
876 match ``calibMaskedImage``.
877 """
878 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
879 if nx != ny: 879 ↛ 880line 879 didn't jump to line 880 because the condition on line 879 was never true
880 raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
881 if nx % 2 != 0: 881 ↛ 882line 881 didn't jump to line 882 because the condition on line 881 was never true
882 raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
883 if nx < 0: 883 ↛ 884line 883 didn't jump to line 884 because the condition on line 883 was never true
884 raise RuntimeError("Calibration maskedImage is larger than raw data.")
886 nEdge = nx//2
887 if nEdge > 0:
888 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
889 SourceDetectionTask.setEdgeBits(
890 rawMaskedImage,
891 replacementMaskedImage.getBBox(),
892 rawMaskedImage.getMask().getPlaneBitMask("EDGE")
893 )
894 else:
895 replacementMaskedImage = rawMaskedImage
897 return replacementMaskedImage
900def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
901 """Apply bias correction in place.
903 Parameters
904 ----------
905 maskedImage : `lsst.afw.image.MaskedImage`
906 Image to process. The image is modified by this method.
907 biasMaskedImage : `lsst.afw.image.MaskedImage`
908 Bias image of the same size as ``maskedImage``
909 trimToFit : `Bool`, optional
910 If True, raw data is symmetrically trimmed to match
911 calibration size.
913 Raises
914 ------
915 RuntimeError
916 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
917 the same size.
919 """
920 if trimToFit:
921 maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
923 if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
924 raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
925 (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
926 maskedImage -= biasMaskedImage
929def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
930 """Apply dark correction in place.
932 Parameters
933 ----------
934 maskedImage : `lsst.afw.image.MaskedImage`
935 Image to process. The image is modified by this method.
936 darkMaskedImage : `lsst.afw.image.MaskedImage`
937 Dark image of the same size as ``maskedImage``.
938 expScale : scalar
939 Dark exposure time for ``maskedImage``.
940 darkScale : scalar
941 Dark exposure time for ``darkMaskedImage``.
942 invert : `Bool`, optional
943 If True, re-add the dark to an already corrected image.
944 trimToFit : `Bool`, optional
945 If True, raw data is symmetrically trimmed to match
946 calibration size.
948 Raises
949 ------
950 RuntimeError
951 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
952 the same size.
954 Notes
955 -----
956 The dark correction is applied by calculating:
957 maskedImage -= dark * expScaling / darkScaling
958 """
959 if trimToFit:
960 maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
962 if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
963 raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
964 (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
966 scale = expScale / darkScale
967 if not invert:
968 maskedImage.scaledMinus(scale, darkMaskedImage)
969 else:
970 maskedImage.scaledPlus(scale, darkMaskedImage)
973def updateVariance(maskedImage, gain, readNoise, replace=True):
974 """Set the variance plane based on the image plane.
976 The maskedImage must have units of `adu` (if gain != 1.0) or
977 electron (if gain == 1.0). This routine will always produce a
978 variance plane in the same units as the image.
980 Parameters
981 ----------
982 maskedImage : `lsst.afw.image.MaskedImage`
983 Image to process. The variance plane is modified.
984 gain : scalar
985 The amplifier gain in electron/adu.
986 readNoise : scalar
987 The amplifier read noise in electron/pixel.
988 replace : `bool`, optional
989 Replace the current variance? If False, the image
990 variance will be added to the current variance plane.
991 """
992 var = maskedImage.variance
993 if replace:
994 var[:, :] = maskedImage.image
995 else:
996 var[:, :] += maskedImage.image
997 var /= gain
998 var += (readNoise/gain)**2
1001def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
1002 """Apply flat correction in place.
1004 Parameters
1005 ----------
1006 maskedImage : `lsst.afw.image.MaskedImage`
1007 Image to process. The image is modified.
1008 flatMaskedImage : `lsst.afw.image.MaskedImage`
1009 Flat image of the same size as ``maskedImage``
1010 scalingType : str
1011 Flat scale computation method. Allowed values are 'MEAN',
1012 'MEDIAN', or 'USER'.
1013 userScale : scalar, optional
1014 Scale to use if ``scalingType='USER'``.
1015 invert : `Bool`, optional
1016 If True, unflatten an already flattened image.
1017 trimToFit : `Bool`, optional
1018 If True, raw data is symmetrically trimmed to match
1019 calibration size.
1021 Raises
1022 ------
1023 RuntimeError
1024 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
1025 the same size or if ``scalingType`` is not an allowed value.
1026 """
1027 if trimToFit:
1028 maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
1030 if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
1031 raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
1032 (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
1034 # Figure out scale from the data
1035 # Ideally the flats are normalized by the calibration product pipeline,
1036 # but this allows some flexibility in the case that the flat is created by
1037 # some other mechanism.
1038 if scalingType in ('MEAN', 'MEDIAN'):
1039 scalingType = afwMath.stringToStatisticsProperty(scalingType)
1040 flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
1041 elif scalingType == 'USER':
1042 flatScale = userScale
1043 else:
1044 raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
1046 if not invert:
1047 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
1048 else:
1049 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
1052def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
1053 """Apply illumination correction in place.
1055 Parameters
1056 ----------
1057 maskedImage : `lsst.afw.image.MaskedImage`
1058 Image to process. The image is modified.
1059 illumMaskedImage : `lsst.afw.image.MaskedImage`
1060 Illumination correction image of the same size as ``maskedImage``.
1061 illumScale : scalar
1062 Scale factor for the illumination correction.
1063 trimToFit : `Bool`, optional
1064 If True, raw data is symmetrically trimmed to match
1065 calibration size.
1067 Raises
1068 ------
1069 RuntimeError
1070 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
1071 the same size.
1072 """
1073 if trimToFit:
1074 maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
1076 if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
1077 raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
1078 (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
1080 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
1083@contextmanager
1084def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
1085 """Context manager that applies and removes gain.
1087 Parameters
1088 ----------
1089 exp : `lsst.afw.image.Exposure`
1090 Exposure to apply/remove gain.
1091 image : `lsst.afw.image.Image`
1092 Image to apply/remove gain.
1093 apply : `bool`
1094 If True, apply and remove the amplifier gain.
1095 gains : `dict` [`str`, `float`], optional
1096 A dictionary, keyed by amplifier name, of the gains to use.
1097 If gains is None, the nominal gains in the amplifier object are used.
1098 invert : `bool`, optional
1099 Invert the gains (e.g. convert electrons to adu temporarily)?
1100 isTrimmed : `bool`, optional
1101 Is this a trimmed exposure?
1103 Yields
1104 ------
1105 exp : `lsst.afw.image.Exposure`
1106 Exposure with the gain applied.
1107 """
1108 # check we have all of them if provided because mixing and matching would
1109 # be a real mess
1110 if gains and apply is True:
1111 ampNames = [amp.getName() for amp in exp.getDetector()]
1112 for ampName in ampNames:
1113 if ampName not in gains.keys(): 1113 ↛ 1114line 1113 didn't jump to line 1114 because the condition on line 1113 was never true
1114 raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
1116 if apply:
1117 ccd = exp.getDetector()
1118 for amp in ccd:
1119 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1120 if gains:
1121 gain = gains[amp.getName()]
1122 else:
1123 gain = amp.getGain()
1124 if invert: 1124 ↛ 1125line 1124 didn't jump to line 1125 because the condition on line 1124 was never true
1125 sim /= gain
1126 else:
1127 sim *= gain
1129 try:
1130 yield exp
1131 finally:
1132 if apply:
1133 ccd = exp.getDetector()
1134 for amp in ccd:
1135 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1136 if gains:
1137 gain = gains[amp.getName()]
1138 else:
1139 gain = amp.getGain()
1140 if invert: 1140 ↛ 1141line 1140 didn't jump to line 1141 because the condition on line 1140 was never true
1141 sim *= gain
1142 else:
1143 sim /= gain
1146def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
1147 sensorTransmission=None, atmosphereTransmission=None):
1148 """Attach a TransmissionCurve to an Exposure, given separate curves for
1149 different components.
1151 Parameters
1152 ----------
1153 exposure : `lsst.afw.image.Exposure`
1154 Exposure object to modify by attaching the product of all given
1155 ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
1156 Must have a valid ``Detector`` attached that matches the detector
1157 associated with sensorTransmission.
1158 opticsTransmission : `lsst.afw.image.TransmissionCurve`
1159 A ``TransmissionCurve`` that represents the throughput of the optics,
1160 to be evaluated in focal-plane coordinates.
1161 filterTransmission : `lsst.afw.image.TransmissionCurve`
1162 A ``TransmissionCurve`` that represents the throughput of the filter
1163 itself, to be evaluated in focal-plane coordinates.
1164 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1165 A ``TransmissionCurve`` that represents the throughput of the sensor
1166 itself, to be evaluated in post-assembly trimmed detector coordinates.
1167 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1168 A ``TransmissionCurve`` that represents the throughput of the
1169 atmosphere, assumed to be spatially constant.
1171 Returns
1172 -------
1173 combined : `lsst.afw.image.TransmissionCurve`
1174 The TransmissionCurve attached to the exposure.
1176 Notes
1177 -----
1178 All ``TransmissionCurve`` arguments are optional; if none are provided, the
1179 attached ``TransmissionCurve`` will have unit transmission everywhere.
1180 """
1181 combined = afwImage.TransmissionCurve.makeIdentity()
1182 if atmosphereTransmission is not None:
1183 combined *= atmosphereTransmission
1184 if opticsTransmission is not None:
1185 combined *= opticsTransmission
1186 if filterTransmission is not None:
1187 combined *= filterTransmission
1188 detector = exposure.getDetector()
1189 fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
1190 toSys=camGeom.PIXELS)
1191 combined = combined.transformedBy(fpToPix)
1192 if sensorTransmission is not None:
1193 combined *= sensorTransmission
1194 exposure.getInfo().setTransmissionCurve(combined)
1195 return combined
1198def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
1199 """Scale an exposure by the amplifier gains.
1201 Parameters
1202 ----------
1203 exposure : `lsst.afw.image.Exposure`
1204 Exposure to process. The image is modified.
1205 normalizeGains : `Bool`, optional
1206 If True, then amplifiers are scaled to force the median of
1207 each amplifier to equal the median of those medians.
1208 ptcGains : `dict`[`str`], optional
1209 Dictionary keyed by amp name containing the PTC gains.
1210 isTrimmed : `bool`, optional
1211 Is the input image trimmed?
1212 """
1213 ccd = exposure.getDetector()
1214 ccdImage = exposure.getMaskedImage()
1216 medians = []
1217 for amp in ccd:
1218 if isTrimmed:
1219 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1220 else:
1221 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1222 if ptcGains: 1222 ↛ 1225line 1222 didn't jump to line 1225 because the condition on line 1222 was always true
1223 sim *= ptcGains[amp.getName()]
1224 else:
1225 sim *= amp.getGain()
1227 if normalizeGains: 1227 ↛ 1228line 1227 didn't jump to line 1228 because the condition on line 1227 was never true
1228 medians.append(numpy.median(sim.getImage().getArray()))
1230 if normalizeGains: 1230 ↛ 1231line 1230 didn't jump to line 1231 because the condition on line 1230 was never true
1231 median = numpy.median(numpy.array(medians))
1232 for index, amp in enumerate(ccd):
1233 if isTrimmed:
1234 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1235 else:
1236 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1237 if medians[index] != 0.0:
1238 sim *= median/medians[index]
1241def widenSaturationTrails(mask):
1242 """Grow the saturation trails by an amount dependent on the width of the
1243 trail.
1245 Parameters
1246 ----------
1247 mask : `lsst.afw.image.Mask`
1248 Mask which will have the saturated areas grown.
1249 """
1251 extraGrowDict = {}
1252 for i in range(1, 6):
1253 extraGrowDict[i] = 0
1254 for i in range(6, 8):
1255 extraGrowDict[i] = 1
1256 for i in range(8, 10):
1257 extraGrowDict[i] = 3
1258 extraGrowMax = 4
1260 if extraGrowMax <= 0: 1260 ↛ 1261line 1260 didn't jump to line 1261 because the condition on line 1260 was never true
1261 return
1263 saturatedBit = mask.getPlaneBitMask("SAT")
1265 xmin, ymin = mask.getBBox().getMin()
1266 width = mask.getWidth()
1268 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1269 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1271 for fp in fpList: 1271 ↛ 1272line 1271 didn't jump to line 1272 because the loop on line 1271 never started
1272 for s in fp.getSpans():
1273 x0, x1 = s.getX0(), s.getX1()
1275 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1276 if extraGrow > 0:
1277 y = s.getY() - ymin
1278 x0 -= xmin + extraGrow
1279 x1 -= xmin - extraGrow
1281 if x0 < 0:
1282 x0 = 0
1283 if x1 >= width - 1:
1284 x1 = width - 1
1286 mask.array[y, x0:x1+1] |= saturatedBit
1289def setBadRegions(exposure, badStatistic="MEDIAN"):
1290 """Set all BAD areas of the chip to the average of the rest of the exposure
1292 Parameters
1293 ----------
1294 exposure : `lsst.afw.image.Exposure`
1295 Exposure to mask. The exposure mask is modified.
1296 badStatistic : `str`, optional
1297 Statistic to use to generate the replacement value from the
1298 image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
1300 Returns
1301 -------
1302 badPixelCount : scalar
1303 Number of bad pixels masked.
1304 badPixelValue : scalar
1305 Value substituted for bad pixels.
1307 Raises
1308 ------
1309 RuntimeError
1310 Raised if `badStatistic` is not an allowed value.
1311 """
1312 if badStatistic == "MEDIAN":
1313 statistic = afwMath.MEDIAN
1314 elif badStatistic == "MEANCLIP":
1315 statistic = afwMath.MEANCLIP
1316 else:
1317 raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
1319 mi = exposure.getMaskedImage()
1320 mask = mi.getMask()
1321 BAD = mask.getPlaneBitMask("BAD")
1322 INTRP = mask.getPlaneBitMask("INTRP")
1324 sctrl = afwMath.StatisticsControl()
1325 sctrl.setAndMask(BAD)
1326 value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
1328 maskArray = mask.getArray()
1329 imageArray = mi.getImage().getArray()
1330 badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
1331 imageArray[:] = numpy.where(badPixels, value, imageArray)
1333 return badPixels.sum(), value
1336def checkFilter(exposure, filterList, log):
1337 """Check to see if an exposure is in a filter specified by a list.
1339 The goal of this is to provide a unified filter checking interface
1340 for all filter dependent stages.
1342 Parameters
1343 ----------
1344 exposure : `lsst.afw.image.Exposure`
1345 Exposure to examine.
1346 filterList : `list` [`str`]
1347 List of physical_filter names to check.
1348 log : `logging.Logger`
1349 Logger to handle messages.
1351 Returns
1352 -------
1353 result : `bool`
1354 True if the exposure's filter is contained in the list.
1355 """
1356 if len(filterList) == 0:
1357 return False
1358 thisFilter = exposure.getFilter()
1359 if thisFilter is None: 1359 ↛ 1360line 1359 didn't jump to line 1360 because the condition on line 1359 was never true
1360 log.warning("No FilterLabel attached to this exposure!")
1361 return False
1363 thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
1364 if thisPhysicalFilter in filterList: 1364 ↛ 1366line 1364 didn't jump to line 1366 because the condition on line 1364 was always true
1365 return True
1366 elif thisFilter.bandLabel in filterList:
1367 if log:
1368 log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
1369 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1370 return True
1371 else:
1372 return False
1375def getPhysicalFilter(filterLabel, log):
1376 """Get the physical filter label associated with the given filterLabel.
1378 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1379 associated with the given ``filterLabel``, the returned label will be
1380 "Unknown".
1382 Parameters
1383 ----------
1384 filterLabel : `lsst.afw.image.FilterLabel`
1385 The `lsst.afw.image.FilterLabel` object from which to derive the
1386 physical filter label.
1387 log : `logging.Logger`
1388 Logger to handle messages.
1390 Returns
1391 -------
1392 physicalFilter : `str`
1393 The value returned by the physicalLabel attribute of ``filterLabel`` if
1394 it exists, otherwise set to \"Unknown\".
1395 """
1396 if filterLabel is None:
1397 physicalFilter = "Unknown"
1398 log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
1399 else:
1400 try:
1401 physicalFilter = filterLabel.physicalLabel
1402 except RuntimeError:
1403 log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
1404 physicalFilter = "Unknown"
1405 return physicalFilter
1408def countMaskedPixels(maskedIm, maskPlane):
1409 """Count the number of pixels in a given mask plane.
1411 Parameters
1412 ----------
1413 maskedIm : `~lsst.afw.image.MaskedImage`
1414 Masked image to examine.
1415 maskPlane : `str`
1416 Name of the mask plane to examine.
1418 Returns
1419 -------
1420 nPix : `int`
1421 Number of pixels in the requested mask plane.
1422 """
1423 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1424 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1425 return nPix
1428def getExposureGains(exposure):
1429 """Get the per-amplifier gains used for this exposure.
1431 Parameters
1432 ----------
1433 exposure : `lsst.afw.image.Exposure`
1434 The exposure to find gains for.
1436 Returns
1437 -------
1438 gains : `dict` [`str` `float`]
1439 Dictionary of gain values, keyed by amplifier name.
1440 Returns empty dict when detector is None.
1441 """
1442 det = exposure.getDetector()
1443 if det is None: 1443 ↛ 1444line 1443 didn't jump to line 1444 because the condition on line 1443 was never true
1444 return dict()
1446 metadata = exposure.getMetadata()
1447 gains = {}
1448 for amp in det:
1449 ampName = amp.getName()
1450 # The key may use the new LSST ISR or the old LSST prefix
1451 if (key1 := f"LSST ISR GAIN {ampName}") in metadata:
1452 gains[ampName] = metadata[key1]
1453 elif (key2 := f"LSST GAIN {ampName}") in metadata:
1454 gains[ampName] = metadata[key2]
1455 else:
1456 gains[ampName] = amp.getGain()
1457 return gains
1460def getExposureReadNoises(exposure):
1461 """Get the per-amplifier read noise used for this exposure.
1463 Parameters
1464 ----------
1465 exposure : `lsst.afw.image.Exposure`
1466 The exposure to find read noise for.
1468 Returns
1469 -------
1470 readnoises : `dict` [`str` `float`]
1471 Dictionary of read noise values, keyed by amplifier name.
1472 Returns empty dict when detector is None.
1473 """
1474 det = exposure.getDetector()
1475 if det is None: 1475 ↛ 1476line 1475 didn't jump to line 1476 because the condition on line 1475 was never true
1476 return dict()
1478 metadata = exposure.getMetadata()
1479 readnoises = {}
1480 for amp in det:
1481 ampName = amp.getName()
1482 # The key may use the new LSST ISR or the old LSST prefix
1483 if (key1 := f"LSST ISR READNOISE {ampName}") in metadata:
1484 readnoises[ampName] = metadata[key1]
1485 elif (key2 := f"LSST READNOISE {ampName}") in metadata:
1486 readnoises[ampName] = metadata[key2]
1487 else:
1488 readnoises[ampName] = amp.getReadNoise()
1489 return readnoises
1492def isTrimmedExposure(exposure):
1493 """Check if the unused pixels (pre-/over-scan pixels) have
1494 been trimmed from an exposure.
1496 Parameters
1497 ----------
1498 exposure : `lsst.afw.image.Exposure`
1499 The exposure to check.
1501 Returns
1502 -------
1503 result : `bool`
1504 True if the image is trimmed, else False.
1505 """
1506 return exposure.getDetector().getBBox() == exposure.getBBox()
1509def isTrimmedImage(image, detector):
1510 """Check if the unused pixels (pre-/over-scan pixels) have
1511 been trimmed from an image
1513 Parameters
1514 ----------
1515 image : `lsst.afw.image.Image`
1516 The image to check.
1517 detector : `lsst.afw.cameraGeom.Detector`
1518 The detector associated with the image.
1520 Returns
1521 -------
1522 result : `bool`
1523 True if the image is trimmed, else False.
1524 """
1525 return detector.getBBox() == image.getBBox()
1528def compareCameraKeywords(
1529 doRaiseOnCalibMismatch,
1530 cameraKeywordsToCompare,
1531 exposureMetadata,
1532 calib,
1533 calibName,
1534 log=None,
1535):
1536 """Compare header keywords to confirm camera states match.
1538 Parameters
1539 ----------
1540 doRaiseOnCalibMismatch : `bool`
1541 Raise on calibration mismatch? Otherwise, log a warning.
1542 cameraKeywordsToCompare : `list` [`str`]
1543 List of camera keywords to compare.
1544 exposureMetadata : `lsst.daf.base.PropertyList`
1545 Header for the exposure being processed.
1546 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1547 Calibration to be applied.
1548 calibName : `str`
1549 Calib type for log message.
1550 log : `logging.Logger`, optional
1551 Logger to handle messages.
1552 """
1553 try:
1554 calibMetadata = calib.metadata
1555 except AttributeError:
1556 return
1558 log = log if log else logging.getLogger(__name__)
1560 missingKeywords = []
1561 for keyword in cameraKeywordsToCompare:
1562 exposureValue = exposureMetadata.get(keyword, None)
1563 if exposureValue is None:
1564 log.debug("Sequencer keyword %s not found in exposure metadata.", keyword)
1565 continue
1567 calibValue = calibMetadata.get(keyword, None)
1569 # We don't log here if there is a missing keyword.
1570 if calibValue is None:
1571 missingKeywords.append(keyword)
1572 continue
1574 if exposureValue != calibValue:
1575 if doRaiseOnCalibMismatch:
1576 raise RuntimeError(
1577 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1578 calibName,
1579 keyword,
1580 exposureValue,
1581 calibValue,
1582 )
1583 else:
1584 log.warning(
1585 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1586 calibName,
1587 keyword,
1588 exposureValue,
1589 calibValue,
1590 )
1591 exposureMetadata[f"ISR {calibName.upper()} SEQUENCER MISMATCH"] = True
1593 if missingKeywords:
1594 log.info(
1595 "Calibration %s missing keywords %s, which were not checked.",
1596 calibName,
1597 ",".join(missingKeywords),
1598 )
1601def symmetrize(inputArray):
1602 """ Copy array over 4 quadrants prior to convolution.
1604 Parameters
1605 ----------
1606 inputarray : `numpy.array`
1607 Input array to symmetrize.
1609 Returns
1610 -------
1611 aSym : `numpy.array`
1612 Symmetrized array.
1613 """
1614 targetShape = list(inputArray.shape)
1615 r1, r2 = inputArray.shape[-1], inputArray.shape[-2]
1616 targetShape[-1] = 2*r1-1
1617 targetShape[-2] = 2*r2-1
1618 aSym = numpy.ndarray(tuple(targetShape))
1619 aSym[..., r2-1:, r1-1:] = inputArray
1620 aSym[..., r2-1:, r1-1::-1] = inputArray
1621 aSym[..., r2-1::-1, r1-1::-1] = inputArray
1622 aSym[..., r2-1::-1, r1-1:] = inputArray
1624 return aSym