Coverage for python/lsst/ip/isr/isrFunctions.py: 87%
560 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-10 09:19 +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, approachRows=20, nSigma=5.0,
285 nRowsCheck=20, nRowsCheckShort=5, minUsablePixelsPerRow=30,
286 minLowFraction=0.03, minLowFractionExtent=0.01, minBlockedRows=5,
287 marginFraction=0.125, saturatedMaskName="SAT", log=None):
288 """Mask DECam-style edge bleeds with low rows next to the read register.
290 For each amplifier that contains a large saturated footprint reaching
291 within ``approachRows`` of its read edge, this function identifies a dip
292 using the per-row fraction of usable pixels more than ``nSigma`` below
293 the clipped sky, measures how far the dip extends inward, and marks
294 those rows (plus a margin) with the saturation mask plane.
296 Rows next to the read edge that are unusable are the signature of a
297 drained serial register, which the overscan correction flags as bad;
298 ``minBlockedRows`` rows are also flagged, despite not having a
299 dip measurement.
301 Parameters
302 ----------
303 exposure : `lsst.afw.image.Exposure`
304 Assembled exposure to mask.
305 satMinArea : `int`, optional
306 Minimum area (pixels) of a saturated footprint to be considered.
307 approachRows : `int`, optional
308 A saturated footprint must come within this many rows of the read edge
309 to trigger the edge bleed check.
310 nSigma : `float`, optional
311 A pixel is "low" if it is more than this many sigma below sky.
312 nRowsCheck : `int`, optional
313 Number of rows from the read edge to check for an edge bleed dip.
314 nRowsCheckShort : `int`, optional
315 A second, shorter window from the read edge to check for a short
316 edge bleed dip.
317 minUsablePixelsPerRow : `int`, optional
318 Minimum number of pixels not saturated, bad or missing to search for
319 dips. Rows with fewer pixels may be "blocked" due to a drained
320 register.
321 minLowFraction : `float`, optional
322 Mean fraction of usable pixels that are low, over either ``nRowsCheck``
323 or ``nRowsCheckShort``, to confirm a dip.
324 minLowFractionExtent : `float`, optional
325 Fraction of usable pixels that are low that each row must exceed
326 to count toward the bleed height; the scan stops after five consecutive
327 rows that do not. Should not exceed ``minLowFraction``.
328 minBlockedRows : `int`, optional
329 Minimum number of blocked rows beyond the normal border at the read
330 edge to be masked if found.
331 marginFraction : `float`, optional
332 Extra rows masked beyond the measured height, as a fraction of
333 that height (plus one row).
334 saturatedMaskName : `str`, optional
335 Mask plane name for saturation.
336 log : `logging.Logger`, optional
337 Logger to handle messages.
339 Raises
340 ------
341 RuntimeError
342 Raised if the exposure has no detector with amplifiers.
343 """
344 log = log if log else logging.getLogger(__name__)
346 detector = exposure.getDetector()
347 if detector is None or len(detector) == 0:
348 raise RuntimeError("Cannot mask DECam edge bleeds: exposure has no detector with amplifiers.")
350 maskedImage = exposure.maskedImage
351 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
352 borderBits = maskedImage.mask.getPlaneBitMask(["BAD", "NO_DATA"])
353 unusableBits = borderBits | saturatedBit
354 excludeStatsBits = unusableBits | maskedImage.mask.getPlaneBitMask("SUSPECT")
356 # find saturated footprints
357 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
358 fpList = afwDetection.FootprintSet(maskedImage.mask, thresh).getFootprints()
359 candidates = []
360 for fp in fpList:
361 if fp.getArea() >= satMinArea:
362 xCore, yCore = fp.getCentroid()
363 candidates.append((int(xCore), int(yCore), fp.getBBox()))
364 if not candidates:
365 log.debug("No saturated footprints large enough for the DECam edge bleed check.")
366 return
368 statsCtrl = afwMath.StatisticsControl()
369 statsCtrl.setAndMask(excludeStatsBits)
371 # End the dip scan after this many consecutive rows without a dip.
372 nRowsStop = 5
374 for amp in detector:
375 ampBBox = amp.getBBox()
376 if not exposure.getBBox().contains(ampBBox):
377 continue
378 readTop = amp.getReadoutCorner() in (camGeom.ReadoutCorner.UL, camGeom.ReadoutCorner.UR)
380 # check to see if saturated footprints reach the readout edge
381 if readTop:
382 reachesReadEdge = any(ampBBox.contains(x, y)
383 and fpBBox.getMaxY() >= ampBBox.getMaxY() - approachRows
384 for x, y, fpBBox in candidates)
385 else:
386 reachesReadEdge = any(ampBBox.contains(x, y)
387 and fpBBox.getMinY() <= ampBBox.getMinY() + approachRows
388 for x, y, fpBBox in candidates)
389 if not reachesReadEdge:
390 continue
392 ampImage = maskedImage[ampBBox]
393 stats = afwMath.makeStatistics(ampImage, afwMath.MEANCLIP | afwMath.STDEVCLIP | afwMath.NPOINT,
394 statsCtrl)
395 if stats.getValue(afwMath.NPOINT) < 1000: 395 ↛ 396line 395 didn't jump to line 396 because the condition on line 395 was never true
396 log.debug("Skipping DECam edge bleed check in amp %s: too few unmasked pixels.",
397 amp.getName())
398 continue
399 sky = stats.getValue(afwMath.MEANCLIP)
400 sigma = stats.getValue(afwMath.STDEVCLIP)
402 usable = (ampImage.mask.array & unusableBits) == 0
403 low = (ampImage.image.array < sky - nSigma*sigma) & usable
404 # Row profiles with index 0 at the read edge, increasing inward.
405 nLow = low.sum(axis=1)
406 nUsable = usable.sum(axis=1)
407 nBorderUsable = ((ampImage.mask.array & borderBits) == 0).sum(axis=1)
408 if readTop:
409 nLow = nLow[::-1]
410 nUsable = nUsable[::-1]
411 else:
412 nBorderUsable = nBorderUsable[::-1]
414 enoughUsable = numpy.nonzero(nUsable >= minUsablePixelsPerRow)[0]
415 if len(enoughUsable) == 0: 415 ↛ 416line 415 didn't jump to line 416 because the condition on line 415 was never true
416 log.debug("Skipping DECam edge bleed check in amp %s: no rows with usable pixels.",
417 amp.getName())
418 continue
419 start = int(enoughUsable[0])
420 borderRows = numpy.nonzero(nBorderUsable >= minUsablePixelsPerRow)[0]
421 border = int(borderRows[0]) if len(borderRows) else 0
422 blocked = start - border
424 fraction = numpy.zeros(len(nLow), dtype=float)
425 good = nUsable >= minUsablePixelsPerRow
426 fraction[good] = nLow[good]/nUsable[good]
428 dipFound = fraction[start:start + nRowsCheck].mean() > minLowFraction
429 if nRowsCheckShort > 0:
430 dipFound |= fraction[start:start + nRowsCheckShort].mean() > minLowFraction
431 if not dipFound and blocked < minBlockedRows:
432 log.debug("Saturated footprint reaches read edge of amp %s but no dip found.",
433 amp.getName())
434 continue
436 height = start if blocked >= minBlockedRows else 0
437 nBelow = 0
438 for i in range(start, len(fraction)): 438 ↛ 446line 438 didn't jump to line 446 because the loop on line 438 didn't complete
439 if fraction[i] > minLowFractionExtent:
440 height = i + 1
441 nBelow = 0
442 else:
443 nBelow += 1
444 if nBelow >= nRowsStop:
445 break
446 height += int(height*marginFraction) + 1
447 height = min(height, ampBBox.getHeight())
449 if readTop:
450 maskedImage.mask[ampBBox].array[-height:, :] |= saturatedBit
451 else:
452 maskedImage.mask[ampBBox].array[:height, :] |= saturatedBit
453 log.info("Found DECam edge bleed in amp %s at the %s read edge (%d blocked rows); masked %d rows.",
454 amp.getName(), "top" if readTop else "bottom", max(blocked, 0), height)
457def maskITLEdgeBleed(ccdExposure, badAmpDict,
458 fpCore, itlEdgeBleedSatMinArea=10000,
459 itlEdgeBleedSatMaxArea=100000,
460 itlEdgeBleedThreshold=5000.,
461 itlEdgeBleedModelConstant=0.02,
462 saturatedMaskName="SAT", log=None):
463 """Mask edge bleeds in ITL detectors.
465 Parameters
466 ----------
467 ccdExposure : `lsst.afw.image.Exposure`
468 Exposure to apply masking to.
469 badAmpDict : `dict` [`str`, `bool`]
470 Dictionary of amplifiers, keyed by name, value is True if
471 amplifier is fully masked.
472 fpCore : `lsst.afw.detection._detection.Footprint`
473 Footprint of saturated core.
474 itlEdgeBleedThreshold : `float`, optional
475 Threshold above median sky background for edge bleed detection
476 (electron units).
477 itlEdgeBleedModelConstant : `float`, optional
478 Constant in the decaying exponential in the edge bleed masking.
479 saturatedMaskName : `str`, optional
480 Mask name for saturation.
481 log : `logging.Logger`, optional
482 Logger to handle messages.
483 """
485 log = log if log else logging.getLogger(__name__)
487 # Get median of amplifier saturation level
488 satLevel = numpy.nanmedian([ccdExposure.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
489 for amp in ccdExposure.getDetector() if not badAmpDict[amp.getName()]])
491 # 1. we check if there are several cores in the footprint:
492 # Get centroid of saturated core
493 xCore, yCore = fpCore.getCentroid()
494 # Turn the Y detector coordinate into Y footprint coordinate
495 yCoreFP = int(yCore) - fpCore.getBBox().getMinY()
496 # Now test if there is one or more cores by checking if the slice at the
497 # center is full of saturated pixels or has several segments of saturated
498 # columns (i.e. several cores with trails)
499 checkCoreNbRow = fpCore.getSpans().asArray()[yCoreFP, :]
500 nbCore = 0
501 indexSwitchTrue = []
502 indexSwitchFalse = []
503 if checkCoreNbRow[0]:
504 # If the slice starts with saturated pixels
505 inSatSegment = True
506 nbCore = 1
507 indexSwitchTrue.append(0)
508 else:
509 # If the slice starts with non saturated pixels
510 inSatSegment = False
512 for i, value in enumerate(checkCoreNbRow):
513 if value:
514 if not inSatSegment:
515 indexSwitchTrue.append(i)
516 # nbCore is the number of detected cores.
517 nbCore += 1
518 inSatSegment = True
519 elif inSatSegment:
520 indexSwitchFalse.append(i)
521 inSatSegment = False
523 # 1. we look for edge bleed in saturated cores in the footprint
524 if nbCore == 2:
525 # we now estimate the x coordinates of the edges of the subfootprint
526 # for each core
527 xEdgesCores = [0]
528 xEdgesCores.append(int((indexSwitchTrue[1] + indexSwitchFalse[0])/2))
529 xEdgesCores.append(fpCore.getSpans().asArray().shape[1])
530 # Get the X and Y footprint coordinates of the cores
531 for i in range(nbCore):
532 subfp = fpCore.getSpans().asArray()[:, xEdgesCores[i]:xEdgesCores[i+1]]
533 xCoreFP = int(xEdgesCores[i] + numpy.argmax(numpy.sum(subfp, axis=0)))
534 # turn into X coordinate in detector space
535 xCore = xCoreFP + fpCore.getBBox().getMinX()
536 # get Y footprint coordinate of the core
537 # by trimming the edges where edge bleeds are potentially dominant
538 if subfp.shape[0] <= 200: 538 ↛ 541line 538 didn't jump to line 541 because the condition on line 538 was always true
539 yCoreFP = int(numpy.argmax(numpy.sum(subfp, axis=1)))
540 else:
541 yCoreFP = int(numpy.argmax(numpy.sum(subfp[100:-100, :],
542 axis=1)))
543 yCoreFP = 100+yCoreFP
545 # Estimate the width of the saturated core
546 widthSat = numpy.sum(subfp[int(yCoreFP), :])
548 subfpArea = numpy.sum(subfp)
549 if subfpArea > itlEdgeBleedSatMinArea and subfpArea < itlEdgeBleedSatMaxArea: 549 ↛ 531line 549 didn't jump to line 531 because the condition on line 549 was always true
550 _applyMaskITLEdgeBleed(ccdExposure, xCore,
551 satLevel, widthSat,
552 itlEdgeBleedThreshold,
553 itlEdgeBleedModelConstant,
554 saturatedMaskName, log)
555 elif nbCore > 2: 555 ↛ 557line 555 didn't jump to line 557 because the condition on line 555 was never true
556 # TODO DM-49736: support N cores in saturated footprint
557 log.warning(
558 "Too many (%d) cores in saturated footprint to mask edge bleeds.",
559 nbCore,
560 )
561 else:
562 # Get centroid of saturated core
563 xCore, yCore = fpCore.getCentroid()
564 # Turn the Y detector coordinate into Y footprint coordinate
565 yCoreFP = yCore - fpCore.getBBox().getMinY()
566 # Get the number of saturated columns around the centroid
567 widthSat = numpy.sum(fpCore.getSpans().asArray()[int(yCoreFP), :])
568 _applyMaskITLEdgeBleed(ccdExposure, xCore,
569 satLevel, widthSat, itlEdgeBleedThreshold,
570 itlEdgeBleedModelConstant, saturatedMaskName, log)
573def _applyMaskITLEdgeBleed(ccdExposure, xCore,
574 satLevel, widthSat,
575 itlEdgeBleedThreshold=5000.,
576 itlEdgeBleedModelConstant=0.03,
577 saturatedMaskName="SAT", log=None):
578 """Apply ITL edge bleed masking model.
580 Parameters
581 ----------
582 ccdExposure : `lsst.afw.image.Exposure`
583 Exposure to apply masking to.
584 xCore: `int`
585 X coordinate of the saturated core.
586 satLevel: `float`
587 Minimum saturation level of the detector.
588 widthSat: `float`
589 Width of the saturated core.
590 itlEdgeBleedThreshold : `float`, optional
591 Threshold above median sky background for edge bleed detection
592 (electron units).
593 itlEdgeBleedModelConstant : `float`, optional
594 Constant in the decaying exponential in the edge bleed masking.
595 saturatedMaskName : `str`, optional
596 Mask name for saturation.
597 log : `logging.Logger`, optional
598 Logger to handle messages.
599 """
600 log = log if log else logging.getLogger(__name__)
602 maskedImage = ccdExposure.maskedImage
603 xmax = maskedImage.image.array.shape[1]
604 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
606 for amp in ccdExposure.getDetector():
607 # Select the 2 top and bottom amplifiers around the saturated
608 # core with a potential edge bleed by selecting the amplifiers
609 # that have the same X coordinate as the saturated core.
610 # As we don't care about the Y coordinate, we set it to the
611 # center of the BBox.
612 yBox = amp.getBBox().getCenter()[1]
613 if amp.getBBox().contains(xCore, yBox):
615 # Get the amp name
616 ampName = amp.getName()
618 # Because in ITLs the edge bleed happens on both edges
619 # of the detector, we make a cutout around
620 # both the top and bottom
621 # edge bleed candidates around the saturated core.
622 # We flip the cutout of the top amplifier
623 # to then work with the same coordinates for both.
624 # The way of selecting top vs bottom amp
625 # is very specific to ITL.
626 if ampName[:2] == 'C1':
627 sliceImage = maskedImage.image.array[:200, :]
628 sliceMask = maskedImage.mask.array[:200, :]
629 elif ampName[:2] == 'C0': 629 ↛ 643line 629 didn't jump to line 643 because the condition on line 629 was always true
630 sliceImage = numpy.flipud(maskedImage.image.array[-200:, :])
631 sliceMask = numpy.flipud(maskedImage.mask.array[-200:, :])
633 # The middle columns of edge bleeds often have
634 # high counts, so we check there is an edge bleed
635 # by looking at a small image up to 50 pixels from the edge
636 # and around the saturated columns
637 # of the saturated core, and checking its median is
638 # above the sky background by itlEdgeBleedThreshold
640 # If the centroid is too close to the edge of the detector
641 # (within 5 pixels), we set the limit to the mean check
642 # to the edge of the detector
643 lowerRangeSmall = int(xCore)-5
644 upperRangeSmall = int(xCore)+5
645 if lowerRangeSmall < 0:
646 lowerRangeSmall = 0
647 if upperRangeSmall > xmax: 647 ↛ 648line 647 didn't jump to line 648 because the condition on line 647 was never true
648 upperRangeSmall = xmax
649 ampImageBG = numpy.median(maskedImage[amp.getBBox()].image.array)
650 edgeMedian = numpy.median(sliceImage[:50, lowerRangeSmall:upperRangeSmall])
651 if edgeMedian > (ampImageBG + itlEdgeBleedThreshold):
653 log.info("Found ITL edge bleed in amp %s, column %d.", ampName, xCore)
655 # We need an estimate of the maximum width
656 # of the edge bleed for our masking model
657 # so we now estimate it by measuring the width of
658 # areas above 60 percent of the saturation level
659 # close to the edge,
660 # in a cutout up to 100 pixels from the edge,
661 # with a width of around the width of an amplifier.
662 subImageXMin = int(xCore)-250
663 subImageXMax = int(xCore)+250
664 if subImageXMin < 0:
665 subImageXMin = 0
666 elif subImageXMax > xmax:
667 subImageXMax = xmax
669 subImage = sliceImage[:100, subImageXMin:subImageXMax]
670 maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
671 axis=1))
673 # Mask edge bleed with a decaying exponential model
674 for y in range(200):
675 edgeBleedHalfWidth = \
676 int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
677 + widthSat)/2.)
678 lowerRange = int(xCore)-edgeBleedHalfWidth
679 upperRange = int(xCore)+edgeBleedHalfWidth
680 # If the edge bleed model goes outside the detector
681 # we set the limit for the masking
682 # to the edge of the detector
683 if lowerRange < 0:
684 lowerRange = 0
685 if upperRange > xmax:
686 upperRange = xmax
687 sliceMask[y, lowerRange:upperRange] |= saturatedBit
690def maskITLSatSag(ccdExposure, fpCore, saturatedMaskName="SAT"):
691 """Mask columns presenting saturation sag in saturated footprints in
692 ITL detectors.
694 Parameters
695 ----------
696 ccdExposure : `lsst.afw.image.Exposure`
697 Exposure to apply masking to.
698 fpCore : `lsst.afw.detection._detection.Footprint`
699 Footprint of saturated core.
700 saturatedMaskName : `str`, optional
701 Mask name for saturation.
702 """
704 # TODO DM-49736: add a flux level check to apply masking
706 maskedImage = ccdExposure.maskedImage
707 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
709 cc = numpy.sum(fpCore.getSpans().asArray(), axis=0)
710 # Mask full columns that have 20 percent of the height of the footprint
711 # saturated
712 columnsToMaskFP = numpy.where(cc > fpCore.getSpans().asArray().shape[0]/5.)
714 columnsToMask = [x + int(fpCore.getBBox().getMinX()) for x in columnsToMaskFP]
715 maskedImage.mask.array[:, columnsToMask] |= saturatedBit
718def maskITLDip(exposure, detectorConfig, maskPlaneNames=["SUSPECT", "ITL_DIP"], log=None):
719 """Add mask bits according to the ITL dip model.
721 Parameters
722 ----------
723 exposure : `lsst.afw.image.Exposure`
724 Exposure to do ITL dip masking.
725 detectorConfig : `lsst.ip.isr.overscanAmpConfig.OverscanDetectorConfig`
726 Configuration for this detector.
727 maskPlaneNames : `list [`str`], optional
728 Name of the ITL Dip mask planes.
729 log : `logging.Logger`, optional
730 If not set, a default logger will be used.
731 """
732 if detectorConfig.itlDipBackgroundFraction == 0.0:
733 # Nothing to do.
734 return
736 if log is None: 736 ↛ 739line 736 didn't jump to line 739 because the condition on line 736 was always true
737 log = logging.getLogger(__name__)
739 thresh = afwDetection.Threshold(
740 exposure.mask.getPlaneBitMask("SAT"),
741 afwDetection.Threshold.BITMASK,
742 )
743 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
745 heights = numpy.asarray([fp.getBBox().getHeight() for fp in fpList])
747 largeHeights, = numpy.where(heights >= detectorConfig.itlDipMinHeight)
749 if len(largeHeights) == 0: 749 ↛ 750line 749 didn't jump to line 750 because the condition on line 749 was never true
750 return
752 # Get the approximate image background.
753 approxBackground = numpy.median(exposure.image.array)
754 maskValue = exposure.mask.getPlaneBitMask(maskPlaneNames)
756 maskBak = exposure.mask.array.copy()
757 nMaskedCols = 0
759 for index in largeHeights:
760 fp = fpList[index]
761 center = fp.getCentroid()
763 nSat = numpy.sum(fp.getSpans().asArray(), axis=0)
764 width = numpy.sum(nSat > detectorConfig.itlDipMinHeight)
766 if width < detectorConfig.itlDipMinWidth:
767 continue
769 width = numpy.clip(width, None, detectorConfig.itlDipMaxWidth)
771 dipMax = detectorConfig.itlDipBackgroundFraction * approxBackground * width
773 # Assume sky-noise dominated; we could add in read noise here.
774 if dipMax < detectorConfig.itlDipMinBackgroundNoiseFraction * numpy.sqrt(approxBackground):
775 continue
777 minCol = int(center.getX() - (detectorConfig.itlDipWidthScale * width) / 2.)
778 maxCol = int(center.getX() + (detectorConfig.itlDipWidthScale * width) / 2.)
779 minCol = numpy.clip(minCol, 0, None)
780 maxCol = numpy.clip(maxCol, None, exposure.mask.array.shape[1] - 1)
782 log.info(
783 "Found ITL dip (width %d; bkg %.2f); masking column %d to %d.",
784 width,
785 approxBackground,
786 minCol,
787 maxCol,
788 )
790 exposure.mask.array[:, minCol: maxCol + 1] |= maskValue
792 nMaskedCols += (maxCol - minCol + 1)
794 if nMaskedCols > detectorConfig.itlDipMaxColsPerImage:
795 log.warning(
796 "Too many (%d) columns would be masked on this image from dip masking; restoring original mask.",
797 nMaskedCols,
798 )
799 exposure.mask.array[:, :] = maskBak
802def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
803 maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
804 """Interpolate over defects identified by a particular set of mask planes.
806 Parameters
807 ----------
808 maskedImage : `lsst.afw.image.MaskedImage`
809 Image to process.
810 fwhm : `float`
811 FWHM of double Gaussian smoothing kernel.
812 growSaturatedFootprints : scalar, optional
813 Number of pixels to grow footprints for saturated pixels.
814 maskNameList : `List` of `str`, optional
815 Mask plane name.
816 fallbackValue : scalar, optional
817 Value of last resort for interpolation.
819 Notes
820 -----
821 The ``fwhm`` parameter is used to create a PSF, but the underlying
822 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
823 not currently make use of this information.
824 """
825 mask = maskedImage.getMask()
827 if growSaturatedFootprints > 0 and "SAT" in maskNameList:
828 # If we are interpolating over an area larger than the original masked
829 # region, we need to expand the original mask bit to the full area to
830 # explain why we interpolated there.
831 growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
833 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
834 fpSet = afwDetection.FootprintSet(mask, thresh)
835 defectList = Defects.fromFootprintList(fpSet.getFootprints())
837 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
838 maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
840 return maskedImage
843def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
844 fallbackValue=None, useLegacyInterp=True):
845 """Mark saturated pixels and optionally interpolate over them
847 Parameters
848 ----------
849 maskedImage : `lsst.afw.image.MaskedImage`
850 Image to process.
851 saturation : scalar
852 Saturation level used as the detection threshold.
853 fwhm : `float`
854 FWHM of double Gaussian smoothing kernel.
855 growFootprints : scalar, optional
856 Number of pixels to grow footprints of detected regions.
857 interpolate : Bool, optional
858 If True, saturated pixels are interpolated over.
859 maskName : str, optional
860 Mask plane name.
861 fallbackValue : scalar, optional
862 Value of last resort for interpolation.
864 Notes
865 -----
866 The ``fwhm`` parameter is used to create a PSF, but the underlying
867 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
868 not currently make use of this information.
869 """
870 defectList = makeThresholdMask(
871 maskedImage=maskedImage,
872 threshold=saturation,
873 growFootprints=growFootprints,
874 maskName=maskName,
875 )
876 if interpolate:
877 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
878 maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
880 return maskedImage
883def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
884 """Compute number of edge trim pixels to match the calibration data.
886 Use the dimension difference between the raw exposure and the
887 calibration exposure to compute the edge trim pixels. This trim
888 is applied symmetrically, with the same number of pixels masked on
889 each side.
891 Parameters
892 ----------
893 rawMaskedImage : `lsst.afw.image.MaskedImage`
894 Image to trim.
895 calibMaskedImage : `lsst.afw.image.MaskedImage`
896 Calibration image to draw new bounding box from.
898 Returns
899 -------
900 replacementMaskedImage : `lsst.afw.image.MaskedImage`
901 ``rawMaskedImage`` trimmed to the appropriate size.
903 Raises
904 ------
905 RuntimeError
906 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
907 match ``calibMaskedImage``.
908 """
909 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
910 if nx != ny: 910 ↛ 911line 910 didn't jump to line 911 because the condition on line 910 was never true
911 raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
912 if nx % 2 != 0: 912 ↛ 913line 912 didn't jump to line 913 because the condition on line 912 was never true
913 raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
914 if nx < 0: 914 ↛ 915line 914 didn't jump to line 915 because the condition on line 914 was never true
915 raise RuntimeError("Calibration maskedImage is larger than raw data.")
917 nEdge = nx//2
918 if nEdge > 0:
919 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
920 SourceDetectionTask.setEdgeBits(
921 rawMaskedImage,
922 replacementMaskedImage.getBBox(),
923 rawMaskedImage.getMask().getPlaneBitMask("EDGE")
924 )
925 else:
926 replacementMaskedImage = rawMaskedImage
928 return replacementMaskedImage
931def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
932 """Apply bias correction in place.
934 Parameters
935 ----------
936 maskedImage : `lsst.afw.image.MaskedImage`
937 Image to process. The image is modified by this method.
938 biasMaskedImage : `lsst.afw.image.MaskedImage`
939 Bias image of the same size as ``maskedImage``
940 trimToFit : `Bool`, optional
941 If True, raw data is symmetrically trimmed to match
942 calibration size.
944 Raises
945 ------
946 RuntimeError
947 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
948 the same size.
950 """
951 if trimToFit:
952 maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
954 if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
955 raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
956 (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
957 maskedImage -= biasMaskedImage
960def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
961 """Apply dark correction in place.
963 Parameters
964 ----------
965 maskedImage : `lsst.afw.image.MaskedImage`
966 Image to process. The image is modified by this method.
967 darkMaskedImage : `lsst.afw.image.MaskedImage`
968 Dark image of the same size as ``maskedImage``.
969 expScale : scalar
970 Dark exposure time for ``maskedImage``.
971 darkScale : scalar
972 Dark exposure time for ``darkMaskedImage``.
973 invert : `Bool`, optional
974 If True, re-add the dark to an already corrected image.
975 trimToFit : `Bool`, optional
976 If True, raw data is symmetrically trimmed to match
977 calibration size.
979 Raises
980 ------
981 RuntimeError
982 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
983 the same size.
985 Notes
986 -----
987 The dark correction is applied by calculating:
988 maskedImage -= dark * expScaling / darkScaling
989 """
990 if trimToFit:
991 maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
993 if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
994 raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
995 (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
997 scale = expScale / darkScale
998 if not invert:
999 maskedImage.scaledMinus(scale, darkMaskedImage)
1000 else:
1001 maskedImage.scaledPlus(scale, darkMaskedImage)
1004def updateVariance(maskedImage, gain, readNoise, replace=True):
1005 """Set the variance plane based on the image plane.
1007 The maskedImage must have units of `adu` (if gain != 1.0) or
1008 electron (if gain == 1.0). This routine will always produce a
1009 variance plane in the same units as the image.
1011 Parameters
1012 ----------
1013 maskedImage : `lsst.afw.image.MaskedImage`
1014 Image to process. The variance plane is modified.
1015 gain : scalar
1016 The amplifier gain in electron/adu.
1017 readNoise : scalar
1018 The amplifier read noise in electron/pixel.
1019 replace : `bool`, optional
1020 Replace the current variance? If False, the image
1021 variance will be added to the current variance plane.
1022 """
1023 var = maskedImage.variance
1024 if replace:
1025 var[:, :] = maskedImage.image
1026 else:
1027 var[:, :] += maskedImage.image
1028 var /= gain
1029 var += (readNoise/gain)**2
1032def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
1033 """Apply flat correction in place.
1035 Parameters
1036 ----------
1037 maskedImage : `lsst.afw.image.MaskedImage`
1038 Image to process. The image is modified.
1039 flatMaskedImage : `lsst.afw.image.MaskedImage`
1040 Flat image of the same size as ``maskedImage``
1041 scalingType : str
1042 Flat scale computation method. Allowed values are 'MEAN',
1043 'MEDIAN', or 'USER'.
1044 userScale : scalar, optional
1045 Scale to use if ``scalingType='USER'``.
1046 invert : `Bool`, optional
1047 If True, unflatten an already flattened image.
1048 trimToFit : `Bool`, optional
1049 If True, raw data is symmetrically trimmed to match
1050 calibration size.
1052 Raises
1053 ------
1054 RuntimeError
1055 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
1056 the same size or if ``scalingType`` is not an allowed value.
1057 """
1058 if trimToFit:
1059 maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
1061 if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
1062 raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
1063 (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
1065 # Figure out scale from the data
1066 # Ideally the flats are normalized by the calibration product pipeline,
1067 # but this allows some flexibility in the case that the flat is created by
1068 # some other mechanism.
1069 if scalingType in ('MEAN', 'MEDIAN'):
1070 scalingType = afwMath.stringToStatisticsProperty(scalingType)
1071 flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
1072 elif scalingType == 'USER':
1073 flatScale = userScale
1074 else:
1075 raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
1077 if not invert:
1078 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
1079 else:
1080 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
1083def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
1084 """Apply illumination correction in place.
1086 Parameters
1087 ----------
1088 maskedImage : `lsst.afw.image.MaskedImage`
1089 Image to process. The image is modified.
1090 illumMaskedImage : `lsst.afw.image.MaskedImage`
1091 Illumination correction image of the same size as ``maskedImage``.
1092 illumScale : scalar
1093 Scale factor for the illumination correction.
1094 trimToFit : `Bool`, optional
1095 If True, raw data is symmetrically trimmed to match
1096 calibration size.
1098 Raises
1099 ------
1100 RuntimeError
1101 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
1102 the same size.
1103 """
1104 if trimToFit:
1105 maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
1107 if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
1108 raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
1109 (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
1111 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
1114@contextmanager
1115def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
1116 """Context manager that applies and removes gain.
1118 Parameters
1119 ----------
1120 exp : `lsst.afw.image.Exposure`
1121 Exposure to apply/remove gain.
1122 image : `lsst.afw.image.Image`
1123 Image to apply/remove gain.
1124 apply : `bool`
1125 If True, apply and remove the amplifier gain.
1126 gains : `dict` [`str`, `float`], optional
1127 A dictionary, keyed by amplifier name, of the gains to use.
1128 If gains is None, the nominal gains in the amplifier object are used.
1129 invert : `bool`, optional
1130 Invert the gains (e.g. convert electrons to adu temporarily)?
1131 isTrimmed : `bool`, optional
1132 Is this a trimmed exposure?
1134 Yields
1135 ------
1136 exp : `lsst.afw.image.Exposure`
1137 Exposure with the gain applied.
1138 """
1139 # check we have all of them if provided because mixing and matching would
1140 # be a real mess
1141 if gains and apply is True:
1142 ampNames = [amp.getName() for amp in exp.getDetector()]
1143 for ampName in ampNames:
1144 if ampName not in gains.keys(): 1144 ↛ 1145line 1144 didn't jump to line 1145 because the condition on line 1144 was never true
1145 raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
1147 if apply:
1148 ccd = exp.getDetector()
1149 for amp in ccd:
1150 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1151 if gains:
1152 gain = gains[amp.getName()]
1153 else:
1154 gain = amp.getGain()
1155 if invert: 1155 ↛ 1156line 1155 didn't jump to line 1156 because the condition on line 1155 was never true
1156 sim /= gain
1157 else:
1158 sim *= gain
1160 try:
1161 yield exp
1162 finally:
1163 if apply:
1164 ccd = exp.getDetector()
1165 for amp in ccd:
1166 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1167 if gains:
1168 gain = gains[amp.getName()]
1169 else:
1170 gain = amp.getGain()
1171 if invert: 1171 ↛ 1172line 1171 didn't jump to line 1172 because the condition on line 1171 was never true
1172 sim *= gain
1173 else:
1174 sim /= gain
1177def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
1178 sensorTransmission=None, atmosphereTransmission=None):
1179 """Attach a TransmissionCurve to an Exposure, given separate curves for
1180 different components.
1182 Parameters
1183 ----------
1184 exposure : `lsst.afw.image.Exposure`
1185 Exposure object to modify by attaching the product of all given
1186 ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
1187 Must have a valid ``Detector`` attached that matches the detector
1188 associated with sensorTransmission.
1189 opticsTransmission : `lsst.afw.image.TransmissionCurve`
1190 A ``TransmissionCurve`` that represents the throughput of the optics,
1191 to be evaluated in focal-plane coordinates.
1192 filterTransmission : `lsst.afw.image.TransmissionCurve`
1193 A ``TransmissionCurve`` that represents the throughput of the filter
1194 itself, to be evaluated in focal-plane coordinates.
1195 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1196 A ``TransmissionCurve`` that represents the throughput of the sensor
1197 itself, to be evaluated in post-assembly trimmed detector coordinates.
1198 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1199 A ``TransmissionCurve`` that represents the throughput of the
1200 atmosphere, assumed to be spatially constant.
1202 Returns
1203 -------
1204 combined : `lsst.afw.image.TransmissionCurve`
1205 The TransmissionCurve attached to the exposure.
1207 Notes
1208 -----
1209 All ``TransmissionCurve`` arguments are optional; if none are provided, the
1210 attached ``TransmissionCurve`` will have unit transmission everywhere.
1211 """
1212 combined = afwImage.TransmissionCurve.makeIdentity()
1213 if atmosphereTransmission is not None:
1214 combined *= atmosphereTransmission
1215 if opticsTransmission is not None:
1216 combined *= opticsTransmission
1217 if filterTransmission is not None:
1218 combined *= filterTransmission
1219 detector = exposure.getDetector()
1220 fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
1221 toSys=camGeom.PIXELS)
1222 combined = combined.transformedBy(fpToPix)
1223 if sensorTransmission is not None:
1224 combined *= sensorTransmission
1225 exposure.getInfo().setTransmissionCurve(combined)
1226 return combined
1229def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
1230 """Scale an exposure by the amplifier gains.
1232 Parameters
1233 ----------
1234 exposure : `lsst.afw.image.Exposure`
1235 Exposure to process. The image is modified.
1236 normalizeGains : `Bool`, optional
1237 If True, then amplifiers are scaled to force the median of
1238 each amplifier to equal the median of those medians.
1239 ptcGains : `dict`[`str`], optional
1240 Dictionary keyed by amp name containing the PTC gains.
1241 isTrimmed : `bool`, optional
1242 Is the input image trimmed?
1243 """
1244 ccd = exposure.getDetector()
1245 ccdImage = exposure.getMaskedImage()
1247 medians = []
1248 for amp in ccd:
1249 if isTrimmed:
1250 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1251 else:
1252 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1253 if ptcGains: 1253 ↛ 1256line 1253 didn't jump to line 1256 because the condition on line 1253 was always true
1254 sim *= ptcGains[amp.getName()]
1255 else:
1256 sim *= amp.getGain()
1258 if normalizeGains: 1258 ↛ 1259line 1258 didn't jump to line 1259 because the condition on line 1258 was never true
1259 medians.append(numpy.median(sim.getImage().getArray()))
1261 if normalizeGains: 1261 ↛ 1262line 1261 didn't jump to line 1262 because the condition on line 1261 was never true
1262 median = numpy.median(numpy.array(medians))
1263 for index, amp in enumerate(ccd):
1264 if isTrimmed:
1265 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1266 else:
1267 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1268 if medians[index] != 0.0:
1269 sim *= median/medians[index]
1272def widenSaturationTrails(mask):
1273 """Grow the saturation trails by an amount dependent on the width of the
1274 trail.
1276 Parameters
1277 ----------
1278 mask : `lsst.afw.image.Mask`
1279 Mask which will have the saturated areas grown.
1280 """
1282 extraGrowDict = {}
1283 for i in range(1, 6):
1284 extraGrowDict[i] = 0
1285 for i in range(6, 8):
1286 extraGrowDict[i] = 1
1287 for i in range(8, 10):
1288 extraGrowDict[i] = 3
1289 extraGrowMax = 4
1291 if extraGrowMax <= 0: 1291 ↛ 1292line 1291 didn't jump to line 1292 because the condition on line 1291 was never true
1292 return
1294 saturatedBit = mask.getPlaneBitMask("SAT")
1296 xmin, ymin = mask.getBBox().getMin()
1297 width = mask.getWidth()
1299 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1300 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1302 for fp in fpList: 1302 ↛ 1303line 1302 didn't jump to line 1303 because the loop on line 1302 never started
1303 for s in fp.getSpans():
1304 x0, x1 = s.getX0(), s.getX1()
1306 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1307 if extraGrow > 0:
1308 y = s.getY() - ymin
1309 x0 -= xmin + extraGrow
1310 x1 -= xmin - extraGrow
1312 if x0 < 0:
1313 x0 = 0
1314 if x1 >= width - 1:
1315 x1 = width - 1
1317 mask.array[y, x0:x1+1] |= saturatedBit
1320def setBadRegions(exposure, badStatistic="MEDIAN"):
1321 """Set all BAD areas of the chip to the average of the rest of the exposure
1323 Parameters
1324 ----------
1325 exposure : `lsst.afw.image.Exposure`
1326 Exposure to mask. The exposure mask is modified.
1327 badStatistic : `str`, optional
1328 Statistic to use to generate the replacement value from the
1329 image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
1331 Returns
1332 -------
1333 badPixelCount : scalar
1334 Number of bad pixels masked.
1335 badPixelValue : scalar
1336 Value substituted for bad pixels.
1338 Raises
1339 ------
1340 RuntimeError
1341 Raised if `badStatistic` is not an allowed value.
1342 """
1343 if badStatistic == "MEDIAN":
1344 statistic = afwMath.MEDIAN
1345 elif badStatistic == "MEANCLIP":
1346 statistic = afwMath.MEANCLIP
1347 else:
1348 raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
1350 mi = exposure.getMaskedImage()
1351 mask = mi.getMask()
1352 BAD = mask.getPlaneBitMask("BAD")
1353 INTRP = mask.getPlaneBitMask("INTRP")
1355 sctrl = afwMath.StatisticsControl()
1356 sctrl.setAndMask(BAD)
1357 value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
1359 maskArray = mask.getArray()
1360 imageArray = mi.getImage().getArray()
1361 badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
1362 imageArray[:] = numpy.where(badPixels, value, imageArray)
1364 return badPixels.sum(), value
1367def checkFilter(exposure, filterList, log):
1368 """Check to see if an exposure is in a filter specified by a list.
1370 The goal of this is to provide a unified filter checking interface
1371 for all filter dependent stages.
1373 Parameters
1374 ----------
1375 exposure : `lsst.afw.image.Exposure`
1376 Exposure to examine.
1377 filterList : `list` [`str`]
1378 List of physical_filter names to check.
1379 log : `logging.Logger`
1380 Logger to handle messages.
1382 Returns
1383 -------
1384 result : `bool`
1385 True if the exposure's filter is contained in the list.
1386 """
1387 if len(filterList) == 0:
1388 return False
1389 thisFilter = exposure.getFilter()
1390 if thisFilter is None: 1390 ↛ 1391line 1390 didn't jump to line 1391 because the condition on line 1390 was never true
1391 log.warning("No FilterLabel attached to this exposure!")
1392 return False
1394 thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
1395 if thisPhysicalFilter in filterList: 1395 ↛ 1397line 1395 didn't jump to line 1397 because the condition on line 1395 was always true
1396 return True
1397 elif thisFilter.bandLabel in filterList:
1398 if log:
1399 log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
1400 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1401 return True
1402 else:
1403 return False
1406def getPhysicalFilter(filterLabel, log):
1407 """Get the physical filter label associated with the given filterLabel.
1409 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1410 associated with the given ``filterLabel``, the returned label will be
1411 "Unknown".
1413 Parameters
1414 ----------
1415 filterLabel : `lsst.afw.image.FilterLabel`
1416 The `lsst.afw.image.FilterLabel` object from which to derive the
1417 physical filter label.
1418 log : `logging.Logger`
1419 Logger to handle messages.
1421 Returns
1422 -------
1423 physicalFilter : `str`
1424 The value returned by the physicalLabel attribute of ``filterLabel`` if
1425 it exists, otherwise set to \"Unknown\".
1426 """
1427 if filterLabel is None:
1428 physicalFilter = "Unknown"
1429 log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
1430 else:
1431 try:
1432 physicalFilter = filterLabel.physicalLabel
1433 except RuntimeError:
1434 log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
1435 physicalFilter = "Unknown"
1436 return physicalFilter
1439def countMaskedPixels(maskedIm, maskPlane):
1440 """Count the number of pixels in a given mask plane.
1442 Parameters
1443 ----------
1444 maskedIm : `~lsst.afw.image.MaskedImage`
1445 Masked image to examine.
1446 maskPlane : `str`
1447 Name of the mask plane to examine.
1449 Returns
1450 -------
1451 nPix : `int`
1452 Number of pixels in the requested mask plane.
1453 """
1454 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1455 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1456 return nPix
1459def getExposureGains(exposure):
1460 """Get the per-amplifier gains used for this exposure.
1462 Parameters
1463 ----------
1464 exposure : `lsst.afw.image.Exposure`
1465 The exposure to find gains for.
1467 Returns
1468 -------
1469 gains : `dict` [`str` `float`]
1470 Dictionary of gain values, keyed by amplifier name.
1471 Returns empty dict when detector is None.
1472 """
1473 det = exposure.getDetector()
1474 if det is None: 1474 ↛ 1475line 1474 didn't jump to line 1475 because the condition on line 1474 was never true
1475 return dict()
1477 metadata = exposure.getMetadata()
1478 gains = {}
1479 for amp in det:
1480 ampName = amp.getName()
1481 # The key may use the new LSST ISR or the old LSST prefix
1482 if (key1 := f"LSST ISR GAIN {ampName}") in metadata:
1483 gains[ampName] = metadata[key1]
1484 elif (key2 := f"LSST GAIN {ampName}") in metadata:
1485 gains[ampName] = metadata[key2]
1486 else:
1487 gains[ampName] = amp.getGain()
1488 return gains
1491def getExposureReadNoises(exposure):
1492 """Get the per-amplifier read noise used for this exposure.
1494 Parameters
1495 ----------
1496 exposure : `lsst.afw.image.Exposure`
1497 The exposure to find read noise for.
1499 Returns
1500 -------
1501 readnoises : `dict` [`str` `float`]
1502 Dictionary of read noise values, keyed by amplifier name.
1503 Returns empty dict when detector is None.
1504 """
1505 det = exposure.getDetector()
1506 if det is None: 1506 ↛ 1507line 1506 didn't jump to line 1507 because the condition on line 1506 was never true
1507 return dict()
1509 metadata = exposure.getMetadata()
1510 readnoises = {}
1511 for amp in det:
1512 ampName = amp.getName()
1513 # The key may use the new LSST ISR or the old LSST prefix
1514 if (key1 := f"LSST ISR READNOISE {ampName}") in metadata:
1515 readnoises[ampName] = metadata[key1]
1516 elif (key2 := f"LSST READNOISE {ampName}") in metadata:
1517 readnoises[ampName] = metadata[key2]
1518 else:
1519 readnoises[ampName] = amp.getReadNoise()
1520 return readnoises
1523def isTrimmedExposure(exposure):
1524 """Check if the unused pixels (pre-/over-scan pixels) have
1525 been trimmed from an exposure.
1527 Parameters
1528 ----------
1529 exposure : `lsst.afw.image.Exposure`
1530 The exposure to check.
1532 Returns
1533 -------
1534 result : `bool`
1535 True if the image is trimmed, else False.
1536 """
1537 return exposure.getDetector().getBBox() == exposure.getBBox()
1540def isTrimmedImage(image, detector):
1541 """Check if the unused pixels (pre-/over-scan pixels) have
1542 been trimmed from an image
1544 Parameters
1545 ----------
1546 image : `lsst.afw.image.Image`
1547 The image to check.
1548 detector : `lsst.afw.cameraGeom.Detector`
1549 The detector associated with the image.
1551 Returns
1552 -------
1553 result : `bool`
1554 True if the image is trimmed, else False.
1555 """
1556 return detector.getBBox() == image.getBBox()
1559def compareCameraKeywords(
1560 doRaiseOnCalibMismatch,
1561 cameraKeywordsToCompare,
1562 exposureMetadata,
1563 calib,
1564 calibName,
1565 log=None,
1566):
1567 """Compare header keywords to confirm camera states match.
1569 Parameters
1570 ----------
1571 doRaiseOnCalibMismatch : `bool`
1572 Raise on calibration mismatch? Otherwise, log a warning.
1573 cameraKeywordsToCompare : `list` [`str`]
1574 List of camera keywords to compare.
1575 exposureMetadata : `lsst.daf.base.PropertyList`
1576 Header for the exposure being processed.
1577 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1578 Calibration to be applied.
1579 calibName : `str`
1580 Calib type for log message.
1581 log : `logging.Logger`, optional
1582 Logger to handle messages.
1583 """
1584 try:
1585 calibMetadata = calib.metadata
1586 except AttributeError:
1587 return
1589 log = log if log else logging.getLogger(__name__)
1591 missingKeywords = []
1592 for keyword in cameraKeywordsToCompare:
1593 exposureValue = exposureMetadata.get(keyword, None)
1594 if exposureValue is None:
1595 log.debug("Sequencer keyword %s not found in exposure metadata.", keyword)
1596 continue
1598 calibValue = calibMetadata.get(keyword, None)
1600 # We don't log here if there is a missing keyword.
1601 if calibValue is None:
1602 missingKeywords.append(keyword)
1603 continue
1605 if exposureValue != calibValue:
1606 if doRaiseOnCalibMismatch:
1607 raise RuntimeError(
1608 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1609 calibName,
1610 keyword,
1611 exposureValue,
1612 calibValue,
1613 )
1614 else:
1615 log.warning(
1616 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1617 calibName,
1618 keyword,
1619 exposureValue,
1620 calibValue,
1621 )
1622 exposureMetadata[f"ISR {calibName.upper()} SEQUENCER MISMATCH"] = True
1624 if missingKeywords:
1625 log.info(
1626 "Calibration %s missing keywords %s, which were not checked.",
1627 calibName,
1628 ",".join(missingKeywords),
1629 )
1632def symmetrize(inputArray):
1633 """ Copy array over 4 quadrants prior to convolution.
1635 Parameters
1636 ----------
1637 inputarray : `numpy.array`
1638 Input array to symmetrize.
1640 Returns
1641 -------
1642 aSym : `numpy.array`
1643 Symmetrized array.
1644 """
1645 targetShape = list(inputArray.shape)
1646 r1, r2 = inputArray.shape[-1], inputArray.shape[-2]
1647 targetShape[-1] = 2*r1-1
1648 targetShape[-2] = 2*r2-1
1649 aSym = numpy.ndarray(tuple(targetShape))
1650 aSym[..., r2-1:, r1-1:] = inputArray
1651 aSym[..., r2-1:, r1-1::-1] = inputArray
1652 aSym[..., r2-1::-1, r1-1::-1] = inputArray
1653 aSym[..., r2-1::-1, r1-1:] = inputArray
1655 return aSym