Coverage for python / lsst / ip / isr / isrFunctions.py: 7%
479 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 08:54 +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 "maskE2VEdgeBleed",
37 "maskITLEdgeBleed",
38 "maskITLSatSag",
39 "maskITLDip",
40 "illuminationCorrection",
41 "interpolateDefectList",
42 "interpolateFromMask",
43 "makeThresholdMask",
44 "saturationCorrection",
45 "setBadRegions",
46 "transposeMaskedImage",
47 "trimToMatchCalibBBox",
48 "updateVariance",
49 "widenSaturationTrails",
50 "getExposureGains",
51 "getExposureReadNoises",
52]
54import logging
55import math
56import numpy
58import lsst.geom
59import lsst.afw.image as afwImage
60import lsst.afw.detection as afwDetection
61import lsst.afw.math as afwMath
62import lsst.meas.algorithms as measAlg
63import lsst.afw.cameraGeom as camGeom
65from lsst.afw.geom import SpanSet, Stencil
66from lsst.meas.algorithms.detection import SourceDetectionTask
68from contextlib import contextmanager
70from .defects import Defects
73def createPsf(fwhm):
74 """Make a double Gaussian PSF.
76 Parameters
77 ----------
78 fwhm : scalar
79 FWHM of double Gaussian smoothing kernel.
81 Returns
82 -------
83 psf : `lsst.meas.algorithms.DoubleGaussianPsf`
84 The created smoothing kernel.
85 """
86 ksize = 4*int(fwhm) + 1
87 return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
90def transposeMaskedImage(maskedImage):
91 """Make a transposed copy of a masked image.
93 Parameters
94 ----------
95 maskedImage : `lsst.afw.image.MaskedImage`
96 Image to process.
98 Returns
99 -------
100 transposed : `lsst.afw.image.MaskedImage`
101 The transposed copy of the input image.
102 """
103 transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
104 transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
105 transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
106 transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
107 return transposed
110def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None,
111 maskNameList=None, useLegacyInterp=True):
112 """Interpolate over defects specified in a defect list.
114 Parameters
115 ----------
116 maskedImage : `lsst.afw.image.MaskedImage`
117 Image to process.
118 defectList : `lsst.meas.algorithms.Defects`
119 List of defects to interpolate over.
120 fwhm : `float`
121 FWHM of double Gaussian smoothing kernel.
122 fallbackValue : scalar, optional
123 Fallback value if an interpolated value cannot be determined.
124 If None, then the clipped mean of the image is used.
125 maskNameList : `list [string]`
126 List of the defects to interpolate over (used for GP interpolator).
127 useLegacyInterp : `bool`
128 Use the legacy interpolation (polynomial interpolation) if True. Use
129 Gaussian Process interpolation if False.
131 Notes
132 -----
133 The ``fwhm`` parameter is used to create a PSF, but the underlying
134 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
135 not currently make use of this information in legacy Interpolation, but use
136 if for the Gaussian Process as an estimation of the correlation lenght.
137 """
138 psf = createPsf(fwhm)
139 if fallbackValue is None:
140 fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
141 if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
142 maskedImage.getMask().addMaskPlane('INTRP')
144 # Hardcoded fwhm value. PSF estimated latter in step1,
145 # not in ISR.
146 if useLegacyInterp:
147 kwargs = {}
148 fwhm = fwhm
149 else:
150 # tested on a dozens of images and looks a good set of
151 # hyperparameters, but cannot guarrenty this is optimal,
152 # need further testing.
153 kwargs = {"bin_spacing": 20,
154 "threshold_dynamic_binning": 2000,
155 "threshold_subdivide": 20000}
156 fwhm = 15
158 measAlg.interpolateOverDefects(maskedImage, psf, defectList,
159 fallbackValue=fallbackValue,
160 useFallbackValueAtEdge=True,
161 fwhm=fwhm,
162 useLegacyInterp=useLegacyInterp,
163 maskNameList=maskNameList, **kwargs)
164 return maskedImage
167def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
168 """Mask pixels based on threshold detection.
170 Parameters
171 ----------
172 maskedImage : `lsst.afw.image.MaskedImage`
173 Image to process. Only the mask plane is updated.
174 threshold : scalar
175 Detection threshold.
176 growFootprints : scalar, optional
177 Number of pixels to grow footprints of detected regions.
178 maskName : str, optional
179 Mask plane name, or list of names to convert
181 Returns
182 -------
183 defectList : `lsst.meas.algorithms.Defects`
184 Defect list constructed from pixels above the threshold.
185 """
186 # find saturated regions
187 thresh = afwDetection.Threshold(threshold)
188 fs = afwDetection.FootprintSet(maskedImage, thresh)
190 if growFootprints > 0:
191 fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
192 fpList = fs.getFootprints()
194 # set mask
195 mask = maskedImage.getMask()
196 bitmask = mask.getPlaneBitMask(maskName)
197 afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
199 return Defects.fromFootprintList(fpList)
202def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
203 """Grow a mask by an amount and add to the requested plane.
205 Parameters
206 ----------
207 mask : `lsst.afw.image.Mask`
208 Mask image to process.
209 radius : scalar
210 Amount to grow the mask.
211 maskNameList : `str` or `list` [`str`]
212 Mask names that should be grown.
213 maskValue : `str`
214 Mask plane to assign the newly masked pixels to.
215 """
216 if radius > 0:
217 spans = SpanSet.fromMask(mask, mask.getPlaneBitMask(maskNameList))
218 # Use MANHATTAN for equivalence with 'isotropic=False` footprint grows,
219 # but CIRCLE is probably better and might be just as fast.
220 spans = spans.dilated(radius, Stencil.MANHATTAN)
221 spans = spans.clippedTo(mask.getBBox())
222 spans.setMask(mask, mask.getPlaneBitMask(maskValue))
225def maskE2VEdgeBleed(exposure, e2vEdgeBleedSatMinArea=10000,
226 e2vEdgeBleedSatMaxArea=100000,
227 e2vEdgeBleedYMax=350,
228 saturatedMaskName="SAT", log=None):
229 """Mask edge bleeds in E2V detectors.
231 Parameters
232 ----------
233 exposure : `lsst.afw.image.Exposure`
234 Exposure to apply masking to.
235 e2vEdgeBleedSatMinArea : `int`, optional
236 Minimum limit of saturated cores footprint area.
237 e2vEdgeBleedSatMaxArea : `int`, optional
238 Maximum limit of saturated cores footprint area.
239 e2vEdgeBleedYMax: `float`, optional
240 Height of edge bleed masking.
241 saturatedMaskName : `str`, optional
242 Mask name for saturation.
243 log : `logging.Logger`, optional
244 Logger to handle messages.
245 """
247 log = log if log else logging.getLogger(__name__)
249 maskedImage = exposure.maskedImage
250 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
252 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
254 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
256 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
257 largeAreas, = numpy.where((satAreas >= e2vEdgeBleedSatMinArea)
258 & (satAreas < e2vEdgeBleedSatMaxArea))
259 for largeAreasIndex in largeAreas:
260 fpCore = fpList[largeAreasIndex]
261 xCore, yCore = fpCore.getCentroid()
262 xCore = int(xCore)
263 yCore = int(yCore)
265 for amp in exposure.getDetector():
266 if amp.getBBox().contains(xCore, yCore):
267 ampName = amp.getName()
268 if ampName[:2] == 'C0':
269 # Check that the footprint reaches the bottom of the
270 # amplifier.
271 if fpCore.getBBox().getMinY() == 0:
272 # This is a large saturation footprint that hits the
273 # edge, and is thus classified as an edge bleed.
275 # TODO DM-50587: Optimize number of rows to mask by
276 # looking at the median signal level as a function of
277 # row number on the right side of the saturation trail.
279 log.info("Found E2V edge bleed in amp %s, column %d.", ampName, xCore)
280 maskedImage.mask[amp.getBBox()].array[:e2vEdgeBleedYMax, :] |= saturatedBit
283def maskITLEdgeBleed(ccdExposure, badAmpDict,
284 fpCore, itlEdgeBleedSatMinArea=10000,
285 itlEdgeBleedSatMaxArea=100000,
286 itlEdgeBleedThreshold=5000.,
287 itlEdgeBleedModelConstant=0.02,
288 saturatedMaskName="SAT", log=None):
289 """Mask edge bleeds in ITL detectors.
291 Parameters
292 ----------
293 ccdExposure : `lsst.afw.image.Exposure`
294 Exposure to apply masking to.
295 badAmpDict : `dict` [`str`, `bool`]
296 Dictionary of amplifiers, keyed by name, value is True if
297 amplifier is fully masked.
298 fpCore : `lsst.afw.detection._detection.Footprint`
299 Footprint of saturated core.
300 itlEdgeBleedThreshold : `float`, optional
301 Threshold above median sky background for edge bleed detection
302 (electron units).
303 itlEdgeBleedModelConstant : `float`, optional
304 Constant in the decaying exponential in the edge bleed masking.
305 saturatedMaskName : `str`, optional
306 Mask name for saturation.
307 log : `logging.Logger`, optional
308 Logger to handle messages.
309 """
311 log = log if log else logging.getLogger(__name__)
313 # Get median of amplifier saturation level
314 satLevel = numpy.nanmedian([ccdExposure.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
315 for amp in ccdExposure.getDetector() if not badAmpDict[amp.getName()]])
317 # 1. we check if there are several cores in the footprint:
318 # Get centroid of saturated core
319 xCore, yCore = fpCore.getCentroid()
320 # Turn the Y detector coordinate into Y footprint coordinate
321 yCoreFP = int(yCore) - fpCore.getBBox().getMinY()
322 # Now test if there is one or more cores by checking if the slice at the
323 # center is full of saturated pixels or has several segments of saturated
324 # columns (i.e. several cores with trails)
325 checkCoreNbRow = fpCore.getSpans().asArray()[yCoreFP, :]
326 nbCore = 0
327 indexSwitchTrue = []
328 indexSwitchFalse = []
329 if checkCoreNbRow[0]:
330 # If the slice starts with saturated pixels
331 inSatSegment = True
332 nbCore = 1
333 indexSwitchTrue.append(0)
334 else:
335 # If the slice starts with non saturated pixels
336 inSatSegment = False
338 for i, value in enumerate(checkCoreNbRow):
339 if value:
340 if not inSatSegment:
341 indexSwitchTrue.append(i)
342 # nbCore is the number of detected cores.
343 nbCore += 1
344 inSatSegment = True
345 elif inSatSegment:
346 indexSwitchFalse.append(i)
347 inSatSegment = False
349 # 1. we look for edge bleed in saturated cores in the footprint
350 if nbCore == 2:
351 # we now estimate the x coordinates of the edges of the subfootprint
352 # for each core
353 xEdgesCores = [0]
354 xEdgesCores.append(int((indexSwitchTrue[1] + indexSwitchFalse[0])/2))
355 xEdgesCores.append(fpCore.getSpans().asArray().shape[1])
356 # Get the X and Y footprint coordinates of the cores
357 for i in range(nbCore):
358 subfp = fpCore.getSpans().asArray()[:, xEdgesCores[i]:xEdgesCores[i+1]]
359 xCoreFP = int(xEdgesCores[i] + numpy.argmax(numpy.sum(subfp, axis=0)))
360 # turn into X coordinate in detector space
361 xCore = xCoreFP + fpCore.getBBox().getMinX()
362 # get Y footprint coordinate of the core
363 # by trimming the edges where edge bleeds are potentially dominant
364 if subfp.shape[0] <= 200:
365 yCoreFP = int(numpy.argmax(numpy.sum(subfp, axis=1)))
366 else:
367 yCoreFP = int(numpy.argmax(numpy.sum(subfp[100:-100, :],
368 axis=1)))
369 yCoreFP = 100+yCoreFP
371 # Estimate the width of the saturated core
372 widthSat = numpy.sum(subfp[int(yCoreFP), :])
374 subfpArea = numpy.sum(subfp)
375 if subfpArea > itlEdgeBleedSatMinArea and subfpArea < itlEdgeBleedSatMaxArea:
376 _applyMaskITLEdgeBleed(ccdExposure, xCore,
377 satLevel, widthSat,
378 itlEdgeBleedThreshold,
379 itlEdgeBleedModelConstant,
380 saturatedMaskName, log)
381 elif nbCore > 2:
382 # TODO DM-49736: support N cores in saturated footprint
383 log.warning(
384 "Too many (%d) cores in saturated footprint to mask edge bleeds.",
385 nbCore,
386 )
387 else:
388 # Get centroid of saturated core
389 xCore, yCore = fpCore.getCentroid()
390 # Turn the Y detector coordinate into Y footprint coordinate
391 yCoreFP = yCore - fpCore.getBBox().getMinY()
392 # Get the number of saturated columns around the centroid
393 widthSat = numpy.sum(fpCore.getSpans().asArray()[int(yCoreFP), :])
394 _applyMaskITLEdgeBleed(ccdExposure, xCore,
395 satLevel, widthSat, itlEdgeBleedThreshold,
396 itlEdgeBleedModelConstant, saturatedMaskName, log)
399def _applyMaskITLEdgeBleed(ccdExposure, xCore,
400 satLevel, widthSat,
401 itlEdgeBleedThreshold=5000.,
402 itlEdgeBleedModelConstant=0.03,
403 saturatedMaskName="SAT", log=None):
404 """Apply ITL edge bleed masking model.
406 Parameters
407 ----------
408 ccdExposure : `lsst.afw.image.Exposure`
409 Exposure to apply masking to.
410 xCore: `int`
411 X coordinate of the saturated core.
412 satLevel: `float`
413 Minimum saturation level of the detector.
414 widthSat: `float`
415 Width of the saturated core.
416 itlEdgeBleedThreshold : `float`, optional
417 Threshold above median sky background for edge bleed detection
418 (electron units).
419 itlEdgeBleedModelConstant : `float`, optional
420 Constant in the decaying exponential in the edge bleed masking.
421 saturatedMaskName : `str`, optional
422 Mask name for saturation.
423 log : `logging.Logger`, optional
424 Logger to handle messages.
425 """
426 log = log if log else logging.getLogger(__name__)
428 maskedImage = ccdExposure.maskedImage
429 xmax = maskedImage.image.array.shape[1]
430 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
432 for amp in ccdExposure.getDetector():
433 # Select the 2 top and bottom amplifiers around the saturated
434 # core with a potential edge bleed by selecting the amplifiers
435 # that have the same X coordinate as the saturated core.
436 # As we don't care about the Y coordinate, we set it to the
437 # center of the BBox.
438 yBox = amp.getBBox().getCenter()[1]
439 if amp.getBBox().contains(xCore, yBox):
441 # Get the amp name
442 ampName = amp.getName()
444 # Because in ITLs the edge bleed happens on both edges
445 # of the detector, we make a cutout around
446 # both the top and bottom
447 # edge bleed candidates around the saturated core.
448 # We flip the cutout of the top amplifier
449 # to then work with the same coordinates for both.
450 # The way of selecting top vs bottom amp
451 # is very specific to ITL.
452 if ampName[:2] == 'C1':
453 sliceImage = maskedImage.image.array[:200, :]
454 sliceMask = maskedImage.mask.array[:200, :]
455 elif ampName[:2] == 'C0':
456 sliceImage = numpy.flipud(maskedImage.image.array[-200:, :])
457 sliceMask = numpy.flipud(maskedImage.mask.array[-200:, :])
459 # The middle columns of edge bleeds often have
460 # high counts, so we check there is an edge bleed
461 # by looking at a small image up to 50 pixels from the edge
462 # and around the saturated columns
463 # of the saturated core, and checking its median is
464 # above the sky background by itlEdgeBleedThreshold
466 # If the centroid is too close to the edge of the detector
467 # (within 5 pixels), we set the limit to the mean check
468 # to the edge of the detector
469 lowerRangeSmall = int(xCore)-5
470 upperRangeSmall = int(xCore)+5
471 if lowerRangeSmall < 0:
472 lowerRangeSmall = 0
473 if upperRangeSmall > xmax:
474 upperRangeSmall = xmax
475 ampImageBG = numpy.median(maskedImage[amp.getBBox()].image.array)
476 edgeMedian = numpy.median(sliceImage[:50, lowerRangeSmall:upperRangeSmall])
477 if edgeMedian > (ampImageBG + itlEdgeBleedThreshold):
479 log.info("Found ITL edge bleed in amp %s, column %d.", ampName, xCore)
481 # We need an estimate of the maximum width
482 # of the edge bleed for our masking model
483 # so we now estimate it by measuring the width of
484 # areas above 60 percent of the saturation level
485 # close to the edge,
486 # in a cutout up to 100 pixels from the edge,
487 # with a width of around the width of an amplifier.
488 subImageXMin = int(xCore)-250
489 subImageXMax = int(xCore)+250
490 if subImageXMin < 0:
491 subImageXMin = 0
492 elif subImageXMax > xmax:
493 subImageXMax = xmax
495 subImage = sliceImage[:100, subImageXMin:subImageXMax]
496 maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
497 axis=1))
499 # Mask edge bleed with a decaying exponential model
500 for y in range(200):
501 edgeBleedHalfWidth = \
502 int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
503 + widthSat)/2.)
504 lowerRange = int(xCore)-edgeBleedHalfWidth
505 upperRange = int(xCore)+edgeBleedHalfWidth
506 # If the edge bleed model goes outside the detector
507 # we set the limit for the masking
508 # to the edge of the detector
509 if lowerRange < 0:
510 lowerRange = 0
511 if upperRange > xmax:
512 upperRange = xmax
513 sliceMask[y, lowerRange:upperRange] |= saturatedBit
516def maskITLSatSag(ccdExposure, fpCore, saturatedMaskName="SAT"):
517 """Mask columns presenting saturation sag in saturated footprints in
518 ITL detectors.
520 Parameters
521 ----------
522 ccdExposure : `lsst.afw.image.Exposure`
523 Exposure to apply masking to.
524 fpCore : `lsst.afw.detection._detection.Footprint`
525 Footprint of saturated core.
526 saturatedMaskName : `str`, optional
527 Mask name for saturation.
528 """
530 # TODO DM-49736: add a flux level check to apply masking
532 maskedImage = ccdExposure.maskedImage
533 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
535 cc = numpy.sum(fpCore.getSpans().asArray(), axis=0)
536 # Mask full columns that have 20 percent of the height of the footprint
537 # saturated
538 columnsToMaskFP = numpy.where(cc > fpCore.getSpans().asArray().shape[0]/5.)
540 columnsToMask = [x + int(fpCore.getBBox().getMinX()) for x in columnsToMaskFP]
541 maskedImage.mask.array[:, columnsToMask] |= saturatedBit
544def maskITLDip(exposure, detectorConfig, maskPlaneNames=["SUSPECT", "ITL_DIP"], log=None):
545 """Add mask bits according to the ITL dip model.
547 Parameters
548 ----------
549 exposure : `lsst.afw.image.Exposure`
550 Exposure to do ITL dip masking.
551 detectorConfig : `lsst.ip.isr.overscanAmpConfig.OverscanDetectorConfig`
552 Configuration for this detector.
553 maskPlaneNames : `list [`str`], optional
554 Name of the ITL Dip mask planes.
555 log : `logging.Logger`, optional
556 If not set, a default logger will be used.
557 """
558 if detectorConfig.itlDipBackgroundFraction == 0.0:
559 # Nothing to do.
560 return
562 if log is None:
563 log = logging.getLogger(__name__)
565 thresh = afwDetection.Threshold(
566 exposure.mask.getPlaneBitMask("SAT"),
567 afwDetection.Threshold.BITMASK,
568 )
569 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
571 heights = numpy.asarray([fp.getBBox().getHeight() for fp in fpList])
573 largeHeights, = numpy.where(heights >= detectorConfig.itlDipMinHeight)
575 if len(largeHeights) == 0:
576 return
578 # Get the approximate image background.
579 approxBackground = numpy.median(exposure.image.array)
580 maskValue = exposure.mask.getPlaneBitMask(maskPlaneNames)
582 maskBak = exposure.mask.array.copy()
583 nMaskedCols = 0
585 for index in largeHeights:
586 fp = fpList[index]
587 center = fp.getCentroid()
589 nSat = numpy.sum(fp.getSpans().asArray(), axis=0)
590 width = numpy.sum(nSat > detectorConfig.itlDipMinHeight)
592 if width < detectorConfig.itlDipMinWidth:
593 continue
595 width = numpy.clip(width, None, detectorConfig.itlDipMaxWidth)
597 dipMax = detectorConfig.itlDipBackgroundFraction * approxBackground * width
599 # Assume sky-noise dominated; we could add in read noise here.
600 if dipMax < detectorConfig.itlDipMinBackgroundNoiseFraction * numpy.sqrt(approxBackground):
601 continue
603 minCol = int(center.getX() - (detectorConfig.itlDipWidthScale * width) / 2.)
604 maxCol = int(center.getX() + (detectorConfig.itlDipWidthScale * width) / 2.)
605 minCol = numpy.clip(minCol, 0, None)
606 maxCol = numpy.clip(maxCol, None, exposure.mask.array.shape[1] - 1)
608 log.info(
609 "Found ITL dip (width %d; bkg %.2f); masking column %d to %d.",
610 width,
611 approxBackground,
612 minCol,
613 maxCol,
614 )
616 exposure.mask.array[:, minCol: maxCol + 1] |= maskValue
618 nMaskedCols += (maxCol - minCol + 1)
620 if nMaskedCols > detectorConfig.itlDipMaxColsPerImage:
621 log.warning(
622 "Too many (%d) columns would be masked on this image from dip masking; restoring original mask.",
623 nMaskedCols,
624 )
625 exposure.mask.array[:, :] = maskBak
628def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
629 maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
630 """Interpolate over defects identified by a particular set of mask planes.
632 Parameters
633 ----------
634 maskedImage : `lsst.afw.image.MaskedImage`
635 Image to process.
636 fwhm : `float`
637 FWHM of double Gaussian smoothing kernel.
638 growSaturatedFootprints : scalar, optional
639 Number of pixels to grow footprints for saturated pixels.
640 maskNameList : `List` of `str`, optional
641 Mask plane name.
642 fallbackValue : scalar, optional
643 Value of last resort for interpolation.
645 Notes
646 -----
647 The ``fwhm`` parameter is used to create a PSF, but the underlying
648 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
649 not currently make use of this information.
650 """
651 mask = maskedImage.getMask()
653 if growSaturatedFootprints > 0 and "SAT" in maskNameList:
654 # If we are interpolating over an area larger than the original masked
655 # region, we need to expand the original mask bit to the full area to
656 # explain why we interpolated there.
657 growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
659 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
660 fpSet = afwDetection.FootprintSet(mask, thresh)
661 defectList = Defects.fromFootprintList(fpSet.getFootprints())
663 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
664 maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
666 return maskedImage
669def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
670 fallbackValue=None, useLegacyInterp=True):
671 """Mark saturated pixels and optionally interpolate over them
673 Parameters
674 ----------
675 maskedImage : `lsst.afw.image.MaskedImage`
676 Image to process.
677 saturation : scalar
678 Saturation level used as the detection threshold.
679 fwhm : `float`
680 FWHM of double Gaussian smoothing kernel.
681 growFootprints : scalar, optional
682 Number of pixels to grow footprints of detected regions.
683 interpolate : Bool, optional
684 If True, saturated pixels are interpolated over.
685 maskName : str, optional
686 Mask plane name.
687 fallbackValue : scalar, optional
688 Value of last resort for interpolation.
690 Notes
691 -----
692 The ``fwhm`` parameter is used to create a PSF, but the underlying
693 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
694 not currently make use of this information.
695 """
696 defectList = makeThresholdMask(
697 maskedImage=maskedImage,
698 threshold=saturation,
699 growFootprints=growFootprints,
700 maskName=maskName,
701 )
702 if interpolate:
703 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
704 maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
706 return maskedImage
709def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
710 """Compute number of edge trim pixels to match the calibration data.
712 Use the dimension difference between the raw exposure and the
713 calibration exposure to compute the edge trim pixels. This trim
714 is applied symmetrically, with the same number of pixels masked on
715 each side.
717 Parameters
718 ----------
719 rawMaskedImage : `lsst.afw.image.MaskedImage`
720 Image to trim.
721 calibMaskedImage : `lsst.afw.image.MaskedImage`
722 Calibration image to draw new bounding box from.
724 Returns
725 -------
726 replacementMaskedImage : `lsst.afw.image.MaskedImage`
727 ``rawMaskedImage`` trimmed to the appropriate size.
729 Raises
730 ------
731 RuntimeError
732 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
733 match ``calibMaskedImage``.
734 """
735 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
736 if nx != ny:
737 raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
738 if nx % 2 != 0:
739 raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
740 if nx < 0:
741 raise RuntimeError("Calibration maskedImage is larger than raw data.")
743 nEdge = nx//2
744 if nEdge > 0:
745 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
746 SourceDetectionTask.setEdgeBits(
747 rawMaskedImage,
748 replacementMaskedImage.getBBox(),
749 rawMaskedImage.getMask().getPlaneBitMask("EDGE")
750 )
751 else:
752 replacementMaskedImage = rawMaskedImage
754 return replacementMaskedImage
757def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
758 """Apply bias correction in place.
760 Parameters
761 ----------
762 maskedImage : `lsst.afw.image.MaskedImage`
763 Image to process. The image is modified by this method.
764 biasMaskedImage : `lsst.afw.image.MaskedImage`
765 Bias image of the same size as ``maskedImage``
766 trimToFit : `Bool`, optional
767 If True, raw data is symmetrically trimmed to match
768 calibration size.
770 Raises
771 ------
772 RuntimeError
773 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
774 the same size.
776 """
777 if trimToFit:
778 maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
780 if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
781 raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
782 (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
783 maskedImage -= biasMaskedImage
786def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
787 """Apply dark correction in place.
789 Parameters
790 ----------
791 maskedImage : `lsst.afw.image.MaskedImage`
792 Image to process. The image is modified by this method.
793 darkMaskedImage : `lsst.afw.image.MaskedImage`
794 Dark image of the same size as ``maskedImage``.
795 expScale : scalar
796 Dark exposure time for ``maskedImage``.
797 darkScale : scalar
798 Dark exposure time for ``darkMaskedImage``.
799 invert : `Bool`, optional
800 If True, re-add the dark to an already corrected image.
801 trimToFit : `Bool`, optional
802 If True, raw data is symmetrically trimmed to match
803 calibration size.
805 Raises
806 ------
807 RuntimeError
808 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
809 the same size.
811 Notes
812 -----
813 The dark correction is applied by calculating:
814 maskedImage -= dark * expScaling / darkScaling
815 """
816 if trimToFit:
817 maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
819 if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
820 raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
821 (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
823 scale = expScale / darkScale
824 if not invert:
825 maskedImage.scaledMinus(scale, darkMaskedImage)
826 else:
827 maskedImage.scaledPlus(scale, darkMaskedImage)
830def updateVariance(maskedImage, gain, readNoise, replace=True):
831 """Set the variance plane based on the image plane.
833 The maskedImage must have units of `adu` (if gain != 1.0) or
834 electron (if gain == 1.0). This routine will always produce a
835 variance plane in the same units as the image.
837 Parameters
838 ----------
839 maskedImage : `lsst.afw.image.MaskedImage`
840 Image to process. The variance plane is modified.
841 gain : scalar
842 The amplifier gain in electron/adu.
843 readNoise : scalar
844 The amplifier read noise in electron/pixel.
845 replace : `bool`, optional
846 Replace the current variance? If False, the image
847 variance will be added to the current variance plane.
848 """
849 var = maskedImage.variance
850 if replace:
851 var[:, :] = maskedImage.image
852 else:
853 var[:, :] += maskedImage.image
854 var /= gain
855 var += (readNoise/gain)**2
858def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
859 """Apply flat correction in place.
861 Parameters
862 ----------
863 maskedImage : `lsst.afw.image.MaskedImage`
864 Image to process. The image is modified.
865 flatMaskedImage : `lsst.afw.image.MaskedImage`
866 Flat image of the same size as ``maskedImage``
867 scalingType : str
868 Flat scale computation method. Allowed values are 'MEAN',
869 'MEDIAN', or 'USER'.
870 userScale : scalar, optional
871 Scale to use if ``scalingType='USER'``.
872 invert : `Bool`, optional
873 If True, unflatten an already flattened image.
874 trimToFit : `Bool`, optional
875 If True, raw data is symmetrically trimmed to match
876 calibration size.
878 Raises
879 ------
880 RuntimeError
881 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
882 the same size or if ``scalingType`` is not an allowed value.
883 """
884 if trimToFit:
885 maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
887 if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
888 raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
889 (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
891 # Figure out scale from the data
892 # Ideally the flats are normalized by the calibration product pipeline,
893 # but this allows some flexibility in the case that the flat is created by
894 # some other mechanism.
895 if scalingType in ('MEAN', 'MEDIAN'):
896 scalingType = afwMath.stringToStatisticsProperty(scalingType)
897 flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
898 elif scalingType == 'USER':
899 flatScale = userScale
900 else:
901 raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
903 if not invert:
904 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
905 else:
906 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
909def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
910 """Apply illumination correction in place.
912 Parameters
913 ----------
914 maskedImage : `lsst.afw.image.MaskedImage`
915 Image to process. The image is modified.
916 illumMaskedImage : `lsst.afw.image.MaskedImage`
917 Illumination correction image of the same size as ``maskedImage``.
918 illumScale : scalar
919 Scale factor for the illumination correction.
920 trimToFit : `Bool`, optional
921 If True, raw data is symmetrically trimmed to match
922 calibration size.
924 Raises
925 ------
926 RuntimeError
927 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
928 the same size.
929 """
930 if trimToFit:
931 maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
933 if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
934 raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
935 (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
937 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
940@contextmanager
941def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
942 """Context manager that applies and removes gain.
944 Parameters
945 ----------
946 exp : `lsst.afw.image.Exposure`
947 Exposure to apply/remove gain.
948 image : `lsst.afw.image.Image`
949 Image to apply/remove gain.
950 apply : `bool`
951 If True, apply and remove the amplifier gain.
952 gains : `dict` [`str`, `float`], optional
953 A dictionary, keyed by amplifier name, of the gains to use.
954 If gains is None, the nominal gains in the amplifier object are used.
955 invert : `bool`, optional
956 Invert the gains (e.g. convert electrons to adu temporarily)?
957 isTrimmed : `bool`, optional
958 Is this a trimmed exposure?
960 Yields
961 ------
962 exp : `lsst.afw.image.Exposure`
963 Exposure with the gain applied.
964 """
965 # check we have all of them if provided because mixing and matching would
966 # be a real mess
967 if gains and apply is True:
968 ampNames = [amp.getName() for amp in exp.getDetector()]
969 for ampName in ampNames:
970 if ampName not in gains.keys():
971 raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
973 if apply:
974 ccd = exp.getDetector()
975 for amp in ccd:
976 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
977 if gains:
978 gain = gains[amp.getName()]
979 else:
980 gain = amp.getGain()
981 if invert:
982 sim /= gain
983 else:
984 sim *= gain
986 try:
987 yield exp
988 finally:
989 if apply:
990 ccd = exp.getDetector()
991 for amp in ccd:
992 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
993 if gains:
994 gain = gains[amp.getName()]
995 else:
996 gain = amp.getGain()
997 if invert:
998 sim *= gain
999 else:
1000 sim /= gain
1003def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
1004 sensorTransmission=None, atmosphereTransmission=None):
1005 """Attach a TransmissionCurve to an Exposure, given separate curves for
1006 different components.
1008 Parameters
1009 ----------
1010 exposure : `lsst.afw.image.Exposure`
1011 Exposure object to modify by attaching the product of all given
1012 ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
1013 Must have a valid ``Detector`` attached that matches the detector
1014 associated with sensorTransmission.
1015 opticsTransmission : `lsst.afw.image.TransmissionCurve`
1016 A ``TransmissionCurve`` that represents the throughput of the optics,
1017 to be evaluated in focal-plane coordinates.
1018 filterTransmission : `lsst.afw.image.TransmissionCurve`
1019 A ``TransmissionCurve`` that represents the throughput of the filter
1020 itself, to be evaluated in focal-plane coordinates.
1021 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1022 A ``TransmissionCurve`` that represents the throughput of the sensor
1023 itself, to be evaluated in post-assembly trimmed detector coordinates.
1024 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1025 A ``TransmissionCurve`` that represents the throughput of the
1026 atmosphere, assumed to be spatially constant.
1028 Returns
1029 -------
1030 combined : `lsst.afw.image.TransmissionCurve`
1031 The TransmissionCurve attached to the exposure.
1033 Notes
1034 -----
1035 All ``TransmissionCurve`` arguments are optional; if none are provided, the
1036 attached ``TransmissionCurve`` will have unit transmission everywhere.
1037 """
1038 combined = afwImage.TransmissionCurve.makeIdentity()
1039 if atmosphereTransmission is not None:
1040 combined *= atmosphereTransmission
1041 if opticsTransmission is not None:
1042 combined *= opticsTransmission
1043 if filterTransmission is not None:
1044 combined *= filterTransmission
1045 detector = exposure.getDetector()
1046 fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
1047 toSys=camGeom.PIXELS)
1048 combined = combined.transformedBy(fpToPix)
1049 if sensorTransmission is not None:
1050 combined *= sensorTransmission
1051 exposure.getInfo().setTransmissionCurve(combined)
1052 return combined
1055def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
1056 """Scale an exposure by the amplifier gains.
1058 Parameters
1059 ----------
1060 exposure : `lsst.afw.image.Exposure`
1061 Exposure to process. The image is modified.
1062 normalizeGains : `Bool`, optional
1063 If True, then amplifiers are scaled to force the median of
1064 each amplifier to equal the median of those medians.
1065 ptcGains : `dict`[`str`], optional
1066 Dictionary keyed by amp name containing the PTC gains.
1067 isTrimmed : `bool`, optional
1068 Is the input image trimmed?
1069 """
1070 ccd = exposure.getDetector()
1071 ccdImage = exposure.getMaskedImage()
1073 medians = []
1074 for amp in ccd:
1075 if isTrimmed:
1076 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1077 else:
1078 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1079 if ptcGains:
1080 sim *= ptcGains[amp.getName()]
1081 else:
1082 sim *= amp.getGain()
1084 if normalizeGains:
1085 medians.append(numpy.median(sim.getImage().getArray()))
1087 if normalizeGains:
1088 median = numpy.median(numpy.array(medians))
1089 for index, amp in enumerate(ccd):
1090 if isTrimmed:
1091 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1092 else:
1093 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1094 if medians[index] != 0.0:
1095 sim *= median/medians[index]
1098def widenSaturationTrails(mask):
1099 """Grow the saturation trails by an amount dependent on the width of the
1100 trail.
1102 Parameters
1103 ----------
1104 mask : `lsst.afw.image.Mask`
1105 Mask which will have the saturated areas grown.
1106 """
1108 extraGrowDict = {}
1109 for i in range(1, 6):
1110 extraGrowDict[i] = 0
1111 for i in range(6, 8):
1112 extraGrowDict[i] = 1
1113 for i in range(8, 10):
1114 extraGrowDict[i] = 3
1115 extraGrowMax = 4
1117 if extraGrowMax <= 0:
1118 return
1120 saturatedBit = mask.getPlaneBitMask("SAT")
1122 xmin, ymin = mask.getBBox().getMin()
1123 width = mask.getWidth()
1125 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1126 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1128 for fp in fpList:
1129 for s in fp.getSpans():
1130 x0, x1 = s.getX0(), s.getX1()
1132 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1133 if extraGrow > 0:
1134 y = s.getY() - ymin
1135 x0 -= xmin + extraGrow
1136 x1 -= xmin - extraGrow
1138 if x0 < 0:
1139 x0 = 0
1140 if x1 >= width - 1:
1141 x1 = width - 1
1143 mask.array[y, x0:x1+1] |= saturatedBit
1146def setBadRegions(exposure, badStatistic="MEDIAN"):
1147 """Set all BAD areas of the chip to the average of the rest of the exposure
1149 Parameters
1150 ----------
1151 exposure : `lsst.afw.image.Exposure`
1152 Exposure to mask. The exposure mask is modified.
1153 badStatistic : `str`, optional
1154 Statistic to use to generate the replacement value from the
1155 image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
1157 Returns
1158 -------
1159 badPixelCount : scalar
1160 Number of bad pixels masked.
1161 badPixelValue : scalar
1162 Value substituted for bad pixels.
1164 Raises
1165 ------
1166 RuntimeError
1167 Raised if `badStatistic` is not an allowed value.
1168 """
1169 if badStatistic == "MEDIAN":
1170 statistic = afwMath.MEDIAN
1171 elif badStatistic == "MEANCLIP":
1172 statistic = afwMath.MEANCLIP
1173 else:
1174 raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
1176 mi = exposure.getMaskedImage()
1177 mask = mi.getMask()
1178 BAD = mask.getPlaneBitMask("BAD")
1179 INTRP = mask.getPlaneBitMask("INTRP")
1181 sctrl = afwMath.StatisticsControl()
1182 sctrl.setAndMask(BAD)
1183 value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
1185 maskArray = mask.getArray()
1186 imageArray = mi.getImage().getArray()
1187 badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
1188 imageArray[:] = numpy.where(badPixels, value, imageArray)
1190 return badPixels.sum(), value
1193def checkFilter(exposure, filterList, log):
1194 """Check to see if an exposure is in a filter specified by a list.
1196 The goal of this is to provide a unified filter checking interface
1197 for all filter dependent stages.
1199 Parameters
1200 ----------
1201 exposure : `lsst.afw.image.Exposure`
1202 Exposure to examine.
1203 filterList : `list` [`str`]
1204 List of physical_filter names to check.
1205 log : `logging.Logger`
1206 Logger to handle messages.
1208 Returns
1209 -------
1210 result : `bool`
1211 True if the exposure's filter is contained in the list.
1212 """
1213 if len(filterList) == 0:
1214 return False
1215 thisFilter = exposure.getFilter()
1216 if thisFilter is None:
1217 log.warning("No FilterLabel attached to this exposure!")
1218 return False
1220 thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
1221 if thisPhysicalFilter in filterList:
1222 return True
1223 elif thisFilter.bandLabel in filterList:
1224 if log:
1225 log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
1226 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1227 return True
1228 else:
1229 return False
1232def getPhysicalFilter(filterLabel, log):
1233 """Get the physical filter label associated with the given filterLabel.
1235 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1236 associated with the given ``filterLabel``, the returned label will be
1237 "Unknown".
1239 Parameters
1240 ----------
1241 filterLabel : `lsst.afw.image.FilterLabel`
1242 The `lsst.afw.image.FilterLabel` object from which to derive the
1243 physical filter label.
1244 log : `logging.Logger`
1245 Logger to handle messages.
1247 Returns
1248 -------
1249 physicalFilter : `str`
1250 The value returned by the physicalLabel attribute of ``filterLabel`` if
1251 it exists, otherwise set to \"Unknown\".
1252 """
1253 if filterLabel is None:
1254 physicalFilter = "Unknown"
1255 log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
1256 else:
1257 try:
1258 physicalFilter = filterLabel.physicalLabel
1259 except RuntimeError:
1260 log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
1261 physicalFilter = "Unknown"
1262 return physicalFilter
1265def countMaskedPixels(maskedIm, maskPlane):
1266 """Count the number of pixels in a given mask plane.
1268 Parameters
1269 ----------
1270 maskedIm : `~lsst.afw.image.MaskedImage`
1271 Masked image to examine.
1272 maskPlane : `str`
1273 Name of the mask plane to examine.
1275 Returns
1276 -------
1277 nPix : `int`
1278 Number of pixels in the requested mask plane.
1279 """
1280 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1281 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1282 return nPix
1285def getExposureGains(exposure):
1286 """Get the per-amplifier gains used for this exposure.
1288 Parameters
1289 ----------
1290 exposure : `lsst.afw.image.Exposure`
1291 The exposure to find gains for.
1293 Returns
1294 -------
1295 gains : `dict` [`str` `float`]
1296 Dictionary of gain values, keyed by amplifier name.
1297 Returns empty dict when detector is None.
1298 """
1299 det = exposure.getDetector()
1300 if det is None:
1301 return dict()
1303 metadata = exposure.getMetadata()
1304 gains = {}
1305 for amp in det:
1306 ampName = amp.getName()
1307 # The key may use the new LSST ISR or the old LSST prefix
1308 if (key1 := f"LSST ISR GAIN {ampName}") in metadata:
1309 gains[ampName] = metadata[key1]
1310 elif (key2 := f"LSST GAIN {ampName}") in metadata:
1311 gains[ampName] = metadata[key2]
1312 else:
1313 gains[ampName] = amp.getGain()
1314 return gains
1317def getExposureReadNoises(exposure):
1318 """Get the per-amplifier read noise used for this exposure.
1320 Parameters
1321 ----------
1322 exposure : `lsst.afw.image.Exposure`
1323 The exposure to find read noise for.
1325 Returns
1326 -------
1327 readnoises : `dict` [`str` `float`]
1328 Dictionary of read noise values, keyed by amplifier name.
1329 Returns empty dict when detector is None.
1330 """
1331 det = exposure.getDetector()
1332 if det is None:
1333 return dict()
1335 metadata = exposure.getMetadata()
1336 readnoises = {}
1337 for amp in det:
1338 ampName = amp.getName()
1339 # The key may use the new LSST ISR or the old LSST prefix
1340 if (key1 := f"LSST ISR READNOISE {ampName}") in metadata:
1341 readnoises[ampName] = metadata[key1]
1342 elif (key2 := f"LSST READNOISE {ampName}") in metadata:
1343 readnoises[ampName] = metadata[key2]
1344 else:
1345 readnoises[ampName] = amp.getReadNoise()
1346 return readnoises
1349def isTrimmedExposure(exposure):
1350 """Check if the unused pixels (pre-/over-scan pixels) have
1351 been trimmed from an exposure.
1353 Parameters
1354 ----------
1355 exposure : `lsst.afw.image.Exposure`
1356 The exposure to check.
1358 Returns
1359 -------
1360 result : `bool`
1361 True if the image is trimmed, else False.
1362 """
1363 return exposure.getDetector().getBBox() == exposure.getBBox()
1366def isTrimmedImage(image, detector):
1367 """Check if the unused pixels (pre-/over-scan pixels) have
1368 been trimmed from an image
1370 Parameters
1371 ----------
1372 image : `lsst.afw.image.Image`
1373 The image to check.
1374 detector : `lsst.afw.cameraGeom.Detector`
1375 The detector associated with the image.
1377 Returns
1378 -------
1379 result : `bool`
1380 True if the image is trimmed, else False.
1381 """
1382 return detector.getBBox() == image.getBBox()
1385def compareCameraKeywords(
1386 doRaiseOnCalibMismatch,
1387 cameraKeywordsToCompare,
1388 exposureMetadata,
1389 calib,
1390 calibName,
1391 log=None,
1392):
1393 """Compare header keywords to confirm camera states match.
1395 Parameters
1396 ----------
1397 doRaiseOnCalibMismatch : `bool`
1398 Raise on calibration mismatch? Otherwise, log a warning.
1399 cameraKeywordsToCompare : `list` [`str`]
1400 List of camera keywords to compare.
1401 exposureMetadata : `lsst.daf.base.PropertyList`
1402 Header for the exposure being processed.
1403 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1404 Calibration to be applied.
1405 calibName : `str`
1406 Calib type for log message.
1407 log : `logging.Logger`, optional
1408 Logger to handle messages.
1409 """
1410 try:
1411 calibMetadata = calib.metadata
1412 except AttributeError:
1413 return
1415 log = log if log else logging.getLogger(__name__)
1417 missingKeywords = []
1418 for keyword in cameraKeywordsToCompare:
1419 exposureValue = exposureMetadata.get(keyword, None)
1420 if exposureValue is None:
1421 log.debug("Sequencer keyword %s not found in exposure metadata.", keyword)
1422 continue
1424 calibValue = calibMetadata.get(keyword, None)
1426 # We don't log here if there is a missing keyword.
1427 if calibValue is None:
1428 missingKeywords.append(keyword)
1429 continue
1431 if exposureValue != calibValue:
1432 if doRaiseOnCalibMismatch:
1433 raise RuntimeError(
1434 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1435 calibName,
1436 keyword,
1437 exposureValue,
1438 calibValue,
1439 )
1440 else:
1441 log.warning(
1442 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1443 calibName,
1444 keyword,
1445 exposureValue,
1446 calibValue,
1447 )
1448 exposureMetadata[f"ISR {calibName.upper()} SEQUENCER MISMATCH"] = True
1450 if missingKeywords:
1451 log.info(
1452 "Calibration %s missing keywords %s, which were not checked.",
1453 calibName,
1454 ",".join(missingKeywords),
1455 )
1458def symmetrize(inputArray):
1459 """ Copy array over 4 quadrants prior to convolution.
1461 Parameters
1462 ----------
1463 inputarray : `numpy.array`
1464 Input array to symmetrize.
1466 Returns
1467 -------
1468 aSym : `numpy.array`
1469 Symmetrized array.
1470 """
1471 targetShape = list(inputArray.shape)
1472 r1, r2 = inputArray.shape[-1], inputArray.shape[-2]
1473 targetShape[-1] = 2*r1-1
1474 targetShape[-2] = 2*r2-1
1475 aSym = numpy.ndarray(tuple(targetShape))
1476 aSym[..., r2-1:, r1-1:] = inputArray
1477 aSym[..., r2-1:, r1-1::-1] = inputArray
1478 aSym[..., r2-1::-1, r1-1::-1] = inputArray
1479 aSym[..., r2-1::-1, r1-1:] = inputArray
1481 return aSym