Coverage for python/lsst/ip/isr/isrFunctions.py: 6%
579 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-24 01:38 -0800
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 "brighterFatterCorrection",
28 "checkFilter",
29 "compareCameraKeywords",
30 "countMaskedPixels",
31 "createPsf",
32 "darkCorrection",
33 "flatCorrection",
34 "fluxConservingBrighterFatterCorrection",
35 "gainContext",
36 "getPhysicalFilter",
37 "growMasks",
38 "maskE2VEdgeBleed",
39 "maskITLEdgeBleed",
40 "maskITLSatSag",
41 "maskITLDip",
42 "illuminationCorrection",
43 "interpolateDefectList",
44 "interpolateFromMask",
45 "makeThresholdMask",
46 "saturationCorrection",
47 "setBadRegions",
48 "transferFlux",
49 "transposeMaskedImage",
50 "trimToMatchCalibBBox",
51 "updateVariance",
52 "widenSaturationTrails",
53 "getExposureGains",
54 "getExposureReadNoises",
55]
57import logging
58import math
59import numpy
60import scipy.signal
62import lsst.geom
63import lsst.afw.image as afwImage
64import lsst.afw.detection as afwDetection
65import lsst.afw.math as afwMath
66import lsst.meas.algorithms as measAlg
67import lsst.afw.cameraGeom as camGeom
69from lsst.afw.geom import SpanSet, Stencil
70from lsst.meas.algorithms.detection import SourceDetectionTask
72from contextlib import contextmanager
74from .defects import Defects
77def createPsf(fwhm):
78 """Make a double Gaussian PSF.
80 Parameters
81 ----------
82 fwhm : scalar
83 FWHM of double Gaussian smoothing kernel.
85 Returns
86 -------
87 psf : `lsst.meas.algorithms.DoubleGaussianPsf`
88 The created smoothing kernel.
89 """
90 ksize = 4*int(fwhm) + 1
91 return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
94def transposeMaskedImage(maskedImage):
95 """Make a transposed copy of a masked image.
97 Parameters
98 ----------
99 maskedImage : `lsst.afw.image.MaskedImage`
100 Image to process.
102 Returns
103 -------
104 transposed : `lsst.afw.image.MaskedImage`
105 The transposed copy of the input image.
106 """
107 transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
108 transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
109 transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
110 transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
111 return transposed
114def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None,
115 maskNameList=None, useLegacyInterp=True):
116 """Interpolate over defects specified in a defect list.
118 Parameters
119 ----------
120 maskedImage : `lsst.afw.image.MaskedImage`
121 Image to process.
122 defectList : `lsst.meas.algorithms.Defects`
123 List of defects to interpolate over.
124 fwhm : `float`
125 FWHM of double Gaussian smoothing kernel.
126 fallbackValue : scalar, optional
127 Fallback value if an interpolated value cannot be determined.
128 If None, then the clipped mean of the image is used.
129 maskNameList : `list [string]`
130 List of the defects to interpolate over (used for GP interpolator).
131 useLegacyInterp : `bool`
132 Use the legacy interpolation (polynomial interpolation) if True. Use
133 Gaussian Process interpolation if False.
135 Notes
136 -----
137 The ``fwhm`` parameter is used to create a PSF, but the underlying
138 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
139 not currently make use of this information in legacy Interpolation, but use
140 if for the Gaussian Process as an estimation of the correlation lenght.
141 """
142 psf = createPsf(fwhm)
143 if fallbackValue is None:
144 fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
145 if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
146 maskedImage.getMask().addMaskPlane('INTRP')
148 # Hardcoded fwhm value. PSF estimated latter in step1,
149 # not in ISR.
150 if useLegacyInterp:
151 kwargs = {}
152 fwhm = fwhm
153 else:
154 # tested on a dozens of images and looks a good set of
155 # hyperparameters, but cannot guarrenty this is optimal,
156 # need further testing.
157 kwargs = {"bin_spacing": 20,
158 "threshold_dynamic_binning": 2000,
159 "threshold_subdivide": 20000}
160 fwhm = 15
162 measAlg.interpolateOverDefects(maskedImage, psf, defectList,
163 fallbackValue=fallbackValue,
164 useFallbackValueAtEdge=True,
165 fwhm=fwhm,
166 useLegacyInterp=useLegacyInterp,
167 maskNameList=maskNameList, **kwargs)
168 return maskedImage
171def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
172 """Mask pixels based on threshold detection.
174 Parameters
175 ----------
176 maskedImage : `lsst.afw.image.MaskedImage`
177 Image to process. Only the mask plane is updated.
178 threshold : scalar
179 Detection threshold.
180 growFootprints : scalar, optional
181 Number of pixels to grow footprints of detected regions.
182 maskName : str, optional
183 Mask plane name, or list of names to convert
185 Returns
186 -------
187 defectList : `lsst.meas.algorithms.Defects`
188 Defect list constructed from pixels above the threshold.
189 """
190 # find saturated regions
191 thresh = afwDetection.Threshold(threshold)
192 fs = afwDetection.FootprintSet(maskedImage, thresh)
194 if growFootprints > 0:
195 fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
196 fpList = fs.getFootprints()
198 # set mask
199 mask = maskedImage.getMask()
200 bitmask = mask.getPlaneBitMask(maskName)
201 afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
203 return Defects.fromFootprintList(fpList)
206def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
207 """Grow a mask by an amount and add to the requested plane.
209 Parameters
210 ----------
211 mask : `lsst.afw.image.Mask`
212 Mask image to process.
213 radius : scalar
214 Amount to grow the mask.
215 maskNameList : `str` or `list` [`str`]
216 Mask names that should be grown.
217 maskValue : `str`
218 Mask plane to assign the newly masked pixels to.
219 """
220 if radius > 0:
221 spans = SpanSet.fromMask(mask, mask.getPlaneBitMask(maskNameList))
222 # Use MANHATTAN for equivalence with 'isotropic=False` footprint grows,
223 # but CIRCLE is probably better and might be just as fast.
224 spans = spans.dilated(radius, Stencil.MANHATTAN)
225 spans = spans.clippedTo(mask.getBBox())
226 spans.setMask(mask, mask.getPlaneBitMask(maskValue))
229def maskE2VEdgeBleed(exposure, e2vEdgeBleedSatMinArea=10000,
230 e2vEdgeBleedSatMaxArea=100000,
231 e2vEdgeBleedYMax=350,
232 saturatedMaskName="SAT", log=None):
233 """Mask edge bleeds in E2V detectors.
235 Parameters
236 ----------
237 exposure : `lsst.afw.image.Exposure`
238 Exposure to apply masking to.
239 e2vEdgeBleedSatMinArea : `int`, optional
240 Minimum limit of saturated cores footprint area.
241 e2vEdgeBleedSatMaxArea : `int`, optional
242 Maximum limit of saturated cores footprint area.
243 e2vEdgeBleedYMax: `float`, optional
244 Height of edge bleed masking.
245 saturatedMaskName : `str`, optional
246 Mask name for saturation.
247 log : `logging.Logger`, optional
248 Logger to handle messages.
249 """
251 log = log if log else logging.getLogger(__name__)
253 maskedImage = exposure.maskedImage
254 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
256 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
258 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
260 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
261 largeAreas, = numpy.where((satAreas >= e2vEdgeBleedSatMinArea)
262 & (satAreas < e2vEdgeBleedSatMaxArea))
263 for largeAreasIndex in largeAreas:
264 fpCore = fpList[largeAreasIndex]
265 xCore, yCore = fpCore.getCentroid()
266 xCore = int(xCore)
267 yCore = int(yCore)
269 for amp in exposure.getDetector():
270 if amp.getBBox().contains(xCore, yCore):
271 ampName = amp.getName()
272 if ampName[:2] == 'C0':
273 # Check that the footprint reaches the bottom of the
274 # amplifier.
275 if fpCore.getBBox().getMinY() == 0:
276 # This is a large saturation footprint that hits the
277 # edge, and is thus classified as an edge bleed.
279 # TODO DM-50587: Optimize number of rows to mask by
280 # looking at the median signal level as a function of
281 # row number on the right side of the saturation trail.
283 log.info("Found E2V edge bleed in amp %s, column %d.", ampName, xCore)
284 maskedImage.mask[amp.getBBox()].array[:e2vEdgeBleedYMax, :] |= saturatedBit
287def maskITLEdgeBleed(ccdExposure, badAmpDict,
288 fpCore, itlEdgeBleedSatMinArea=10000,
289 itlEdgeBleedSatMaxArea=100000,
290 itlEdgeBleedThreshold=5000.,
291 itlEdgeBleedModelConstant=0.02,
292 saturatedMaskName="SAT", log=None):
293 """Mask edge bleeds in ITL detectors.
295 Parameters
296 ----------
297 ccdExposure : `lsst.afw.image.Exposure`
298 Exposure to apply masking to.
299 badAmpDict : `dict` [`str`, `bool`]
300 Dictionary of amplifiers, keyed by name, value is True if
301 amplifier is fully masked.
302 fpCore : `lsst.afw.detection._detection.Footprint`
303 Footprint of saturated core.
304 itlEdgeBleedThreshold : `float`, optional
305 Threshold above median sky background for edge bleed detection
306 (electron units).
307 itlEdgeBleedModelConstant : `float`, optional
308 Constant in the decaying exponential in the edge bleed masking.
309 saturatedMaskName : `str`, optional
310 Mask name for saturation.
311 log : `logging.Logger`, optional
312 Logger to handle messages.
313 """
315 log = log if log else logging.getLogger(__name__)
317 # Get median of amplifier saturation level
318 satLevel = numpy.nanmedian([ccdExposure.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
319 for amp in ccdExposure.getDetector() if not badAmpDict[amp.getName()]])
321 # 1. we check if there are several cores in the footprint:
322 # Get centroid of saturated core
323 xCore, yCore = fpCore.getCentroid()
324 # Turn the Y detector coordinate into Y footprint coordinate
325 yCoreFP = int(yCore) - fpCore.getBBox().getMinY()
326 # Now test if there is one or more cores by checking if the slice at the
327 # center is full of saturated pixels or has several segments of saturated
328 # columns (i.e. several cores with trails)
329 checkCoreNbRow = fpCore.getSpans().asArray()[yCoreFP, :]
330 nbCore = 0
331 indexSwitchTrue = []
332 indexSwitchFalse = []
333 if checkCoreNbRow[0]:
334 # If the slice starts with saturated pixels
335 inSatSegment = True
336 nbCore = 1
337 indexSwitchTrue.append(0)
338 else:
339 # If the slice starts with non saturated pixels
340 inSatSegment = False
342 for i, value in enumerate(checkCoreNbRow):
343 if value:
344 if not inSatSegment:
345 indexSwitchTrue.append(i)
346 # nbCore is the number of detected cores.
347 nbCore += 1
348 inSatSegment = True
349 elif inSatSegment:
350 indexSwitchFalse.append(i)
351 inSatSegment = False
353 # 1. we look for edge bleed in saturated cores in the footprint
354 if nbCore == 2:
355 # we now estimate the x coordinates of the edges of the subfootprint
356 # for each core
357 xEdgesCores = [0]
358 xEdgesCores.append(int((indexSwitchTrue[1] + indexSwitchFalse[0])/2))
359 xEdgesCores.append(fpCore.getSpans().asArray().shape[1])
360 # Get the X and Y footprint coordinates of the cores
361 for i in range(nbCore):
362 subfp = fpCore.getSpans().asArray()[:, xEdgesCores[i]:xEdgesCores[i+1]]
363 xCoreFP = int(xEdgesCores[i] + numpy.argmax(numpy.sum(subfp, axis=0)))
364 # turn into X coordinate in detector space
365 xCore = xCoreFP + fpCore.getBBox().getMinX()
366 # get Y footprint coordinate of the core
367 # by trimming the edges where edge bleeds are potentially dominant
368 if subfp.shape[0] <= 200:
369 yCoreFP = int(numpy.argmax(numpy.sum(subfp, axis=1)))
370 else:
371 yCoreFP = int(numpy.argmax(numpy.sum(subfp[100:-100, :],
372 axis=1)))
373 yCoreFP = 100+yCoreFP
375 # Estimate the width of the saturated core
376 widthSat = numpy.sum(subfp[int(yCoreFP), :])
378 subfpArea = numpy.sum(subfp)
379 if subfpArea > itlEdgeBleedSatMinArea and subfpArea < itlEdgeBleedSatMaxArea:
380 _applyMaskITLEdgeBleed(ccdExposure, xCore,
381 satLevel, widthSat,
382 itlEdgeBleedThreshold,
383 itlEdgeBleedModelConstant,
384 saturatedMaskName, log)
385 elif nbCore > 2:
386 # TODO DM-49736: support N cores in saturated footprint
387 log.warning(
388 "Too many (%d) cores in saturated footprint to mask edge bleeds.",
389 nbCore,
390 )
391 else:
392 # Get centroid of saturated core
393 xCore, yCore = fpCore.getCentroid()
394 # Turn the Y detector coordinate into Y footprint coordinate
395 yCoreFP = yCore - fpCore.getBBox().getMinY()
396 # Get the number of saturated columns around the centroid
397 widthSat = numpy.sum(fpCore.getSpans().asArray()[int(yCoreFP), :])
398 _applyMaskITLEdgeBleed(ccdExposure, xCore,
399 satLevel, widthSat, itlEdgeBleedThreshold,
400 itlEdgeBleedModelConstant, saturatedMaskName, log)
403def _applyMaskITLEdgeBleed(ccdExposure, xCore,
404 satLevel, widthSat,
405 itlEdgeBleedThreshold=5000.,
406 itlEdgeBleedModelConstant=0.03,
407 saturatedMaskName="SAT", log=None):
408 """Apply ITL edge bleed masking model.
410 Parameters
411 ----------
412 ccdExposure : `lsst.afw.image.Exposure`
413 Exposure to apply masking to.
414 xCore: `int`
415 X coordinate of the saturated core.
416 satLevel: `float`
417 Minimum saturation level of the detector.
418 widthSat: `float`
419 Width of the saturated core.
420 itlEdgeBleedThreshold : `float`, optional
421 Threshold above median sky background for edge bleed detection
422 (electron units).
423 itlEdgeBleedModelConstant : `float`, optional
424 Constant in the decaying exponential in the edge bleed masking.
425 saturatedMaskName : `str`, optional
426 Mask name for saturation.
427 log : `logging.Logger`, optional
428 Logger to handle messages.
429 """
430 log = log if log else logging.getLogger(__name__)
432 maskedImage = ccdExposure.maskedImage
433 xmax = maskedImage.image.array.shape[1]
434 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
436 for amp in ccdExposure.getDetector():
437 # Select the 2 top and bottom amplifiers around the saturated
438 # core with a potential edge bleed by selecting the amplifiers
439 # that have the same X coordinate as the saturated core.
440 # As we don't care about the Y coordinate, we set it to the
441 # center of the BBox.
442 yBox = amp.getBBox().getCenter()[1]
443 if amp.getBBox().contains(xCore, yBox):
445 # Get the amp name
446 ampName = amp.getName()
448 # Because in ITLs the edge bleed happens on both edges
449 # of the detector, we make a cutout around
450 # both the top and bottom
451 # edge bleed candidates around the saturated core.
452 # We flip the cutout of the top amplifier
453 # to then work with the same coordinates for both.
454 # The way of selecting top vs bottom amp
455 # is very specific to ITL.
456 if ampName[:2] == 'C1':
457 sliceImage = maskedImage.image.array[:200, :]
458 sliceMask = maskedImage.mask.array[:200, :]
459 elif ampName[:2] == 'C0':
460 sliceImage = numpy.flipud(maskedImage.image.array[-200:, :])
461 sliceMask = numpy.flipud(maskedImage.mask.array[-200:, :])
463 # The middle columns of edge bleeds often have
464 # high counts, so we check there is an edge bleed
465 # by looking at a small image up to 50 pixels from the edge
466 # and around the saturated columns
467 # of the saturated core, and checking its median is
468 # above the sky background by itlEdgeBleedThreshold
470 # If the centroid is too close to the edge of the detector
471 # (within 5 pixels), we set the limit to the mean check
472 # to the edge of the detector
473 lowerRangeSmall = int(xCore)-5
474 upperRangeSmall = int(xCore)+5
475 if lowerRangeSmall < 0:
476 lowerRangeSmall = 0
477 if upperRangeSmall > xmax:
478 upperRangeSmall = xmax
479 ampImageBG = numpy.median(maskedImage[amp.getBBox()].image.array)
480 edgeMedian = numpy.median(sliceImage[:50, lowerRangeSmall:upperRangeSmall])
481 if edgeMedian > (ampImageBG + itlEdgeBleedThreshold):
483 log.info("Found ITL edge bleed in amp %s, column %d.", ampName, xCore)
485 # We need an estimate of the maximum width
486 # of the edge bleed for our masking model
487 # so we now estimate it by measuring the width of
488 # areas above 60 percent of the saturation level
489 # close to the edge,
490 # in a cutout up to 100 pixels from the edge,
491 # with a width of around the width of an amplifier.
492 subImageXMin = int(xCore)-250
493 subImageXMax = int(xCore)+250
494 if subImageXMin < 0:
495 subImageXMin = 0
496 elif subImageXMax > xmax:
497 subImageXMax = xmax
499 subImage = sliceImage[:100, subImageXMin:subImageXMax]
500 maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
501 axis=1))
503 # Mask edge bleed with a decaying exponential model
504 for y in range(200):
505 edgeBleedHalfWidth = \
506 int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
507 + widthSat)/2.)
508 lowerRange = int(xCore)-edgeBleedHalfWidth
509 upperRange = int(xCore)+edgeBleedHalfWidth
510 # If the edge bleed model goes outside the detector
511 # we set the limit for the masking
512 # to the edge of the detector
513 if lowerRange < 0:
514 lowerRange = 0
515 if upperRange > xmax:
516 upperRange = xmax
517 sliceMask[y, lowerRange:upperRange] |= saturatedBit
520def maskITLSatSag(ccdExposure, fpCore, saturatedMaskName="SAT"):
521 """Mask columns presenting saturation sag in saturated footprints in
522 ITL detectors.
524 Parameters
525 ----------
526 ccdExposure : `lsst.afw.image.Exposure`
527 Exposure to apply masking to.
528 fpCore : `lsst.afw.detection._detection.Footprint`
529 Footprint of saturated core.
530 saturatedMaskName : `str`, optional
531 Mask name for saturation.
532 """
534 # TODO DM-49736: add a flux level check to apply masking
536 maskedImage = ccdExposure.maskedImage
537 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
539 cc = numpy.sum(fpCore.getSpans().asArray(), axis=0)
540 # Mask full columns that have 20 percent of the height of the footprint
541 # saturated
542 columnsToMaskFP = numpy.where(cc > fpCore.getSpans().asArray().shape[0]/5.)
544 columnsToMask = [x + int(fpCore.getBBox().getMinX()) for x in columnsToMaskFP]
545 maskedImage.mask.array[:, columnsToMask] |= saturatedBit
548def maskITLDip(exposure, detectorConfig, maskPlaneNames=["SUSPECT", "ITL_DIP"], log=None):
549 """Add mask bits according to the ITL dip model.
551 Parameters
552 ----------
553 exposure : `lsst.afw.image.Exposure`
554 Exposure to do ITL dip masking.
555 detectorConfig : `lsst.ip.isr.overscanAmpConfig.OverscanDetectorConfig`
556 Configuration for this detector.
557 maskPlaneNames : `list [`str`], optional
558 Name of the ITL Dip mask planes.
559 log : `logging.Logger`, optional
560 If not set, a default logger will be used.
561 """
562 if detectorConfig.itlDipBackgroundFraction == 0.0:
563 # Nothing to do.
564 return
566 if log is None:
567 log = logging.getLogger(__name__)
569 thresh = afwDetection.Threshold(
570 exposure.mask.getPlaneBitMask("SAT"),
571 afwDetection.Threshold.BITMASK,
572 )
573 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
575 heights = numpy.asarray([fp.getBBox().getHeight() for fp in fpList])
577 largeHeights, = numpy.where(heights >= detectorConfig.itlDipMinHeight)
579 if len(largeHeights) == 0:
580 return
582 # Get the approximate image background.
583 approxBackground = numpy.median(exposure.image.array)
584 maskValue = exposure.mask.getPlaneBitMask(maskPlaneNames)
586 maskBak = exposure.mask.array.copy()
587 nMaskedCols = 0
589 for index in largeHeights:
590 fp = fpList[index]
591 center = fp.getCentroid()
593 nSat = numpy.sum(fp.getSpans().asArray(), axis=0)
594 width = numpy.sum(nSat > detectorConfig.itlDipMinHeight)
596 if width < detectorConfig.itlDipMinWidth:
597 continue
599 width = numpy.clip(width, None, detectorConfig.itlDipMaxWidth)
601 dipMax = detectorConfig.itlDipBackgroundFraction * approxBackground * width
603 # Assume sky-noise dominated; we could add in read noise here.
604 if dipMax < detectorConfig.itlDipMinBackgroundNoiseFraction * numpy.sqrt(approxBackground):
605 continue
607 minCol = int(center.getX() - (detectorConfig.itlDipWidthScale * width) / 2.)
608 maxCol = int(center.getX() + (detectorConfig.itlDipWidthScale * width) / 2.)
609 minCol = numpy.clip(minCol, 0, None)
610 maxCol = numpy.clip(maxCol, None, exposure.mask.array.shape[1] - 1)
612 log.info(
613 "Found ITL dip (width %d; bkg %.2f); masking column %d to %d.",
614 width,
615 approxBackground,
616 minCol,
617 maxCol,
618 )
620 exposure.mask.array[:, minCol: maxCol + 1] |= maskValue
622 nMaskedCols += (maxCol - minCol + 1)
624 if nMaskedCols > detectorConfig.itlDipMaxColsPerImage:
625 log.warning(
626 "Too many (%d) columns would be masked on this image from dip masking; restoring original mask.",
627 nMaskedCols,
628 )
629 exposure.mask.array[:, :] = maskBak
632def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
633 maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
634 """Interpolate over defects identified by a particular set of mask planes.
636 Parameters
637 ----------
638 maskedImage : `lsst.afw.image.MaskedImage`
639 Image to process.
640 fwhm : `float`
641 FWHM of double Gaussian smoothing kernel.
642 growSaturatedFootprints : scalar, optional
643 Number of pixels to grow footprints for saturated pixels.
644 maskNameList : `List` of `str`, optional
645 Mask plane name.
646 fallbackValue : scalar, optional
647 Value of last resort for interpolation.
649 Notes
650 -----
651 The ``fwhm`` parameter is used to create a PSF, but the underlying
652 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
653 not currently make use of this information.
654 """
655 mask = maskedImage.getMask()
657 if growSaturatedFootprints > 0 and "SAT" in maskNameList:
658 # If we are interpolating over an area larger than the original masked
659 # region, we need to expand the original mask bit to the full area to
660 # explain why we interpolated there.
661 growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
663 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
664 fpSet = afwDetection.FootprintSet(mask, thresh)
665 defectList = Defects.fromFootprintList(fpSet.getFootprints())
667 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
668 maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
670 return maskedImage
673def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
674 fallbackValue=None, useLegacyInterp=True):
675 """Mark saturated pixels and optionally interpolate over them
677 Parameters
678 ----------
679 maskedImage : `lsst.afw.image.MaskedImage`
680 Image to process.
681 saturation : scalar
682 Saturation level used as the detection threshold.
683 fwhm : `float`
684 FWHM of double Gaussian smoothing kernel.
685 growFootprints : scalar, optional
686 Number of pixels to grow footprints of detected regions.
687 interpolate : Bool, optional
688 If True, saturated pixels are interpolated over.
689 maskName : str, optional
690 Mask plane name.
691 fallbackValue : scalar, optional
692 Value of last resort for interpolation.
694 Notes
695 -----
696 The ``fwhm`` parameter is used to create a PSF, but the underlying
697 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
698 not currently make use of this information.
699 """
700 defectList = makeThresholdMask(
701 maskedImage=maskedImage,
702 threshold=saturation,
703 growFootprints=growFootprints,
704 maskName=maskName,
705 )
706 if interpolate:
707 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
708 maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
710 return maskedImage
713def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
714 """Compute number of edge trim pixels to match the calibration data.
716 Use the dimension difference between the raw exposure and the
717 calibration exposure to compute the edge trim pixels. This trim
718 is applied symmetrically, with the same number of pixels masked on
719 each side.
721 Parameters
722 ----------
723 rawMaskedImage : `lsst.afw.image.MaskedImage`
724 Image to trim.
725 calibMaskedImage : `lsst.afw.image.MaskedImage`
726 Calibration image to draw new bounding box from.
728 Returns
729 -------
730 replacementMaskedImage : `lsst.afw.image.MaskedImage`
731 ``rawMaskedImage`` trimmed to the appropriate size.
733 Raises
734 ------
735 RuntimeError
736 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
737 match ``calibMaskedImage``.
738 """
739 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
740 if nx != ny:
741 raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
742 if nx % 2 != 0:
743 raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
744 if nx < 0:
745 raise RuntimeError("Calibration maskedImage is larger than raw data.")
747 nEdge = nx//2
748 if nEdge > 0:
749 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
750 SourceDetectionTask.setEdgeBits(
751 rawMaskedImage,
752 replacementMaskedImage.getBBox(),
753 rawMaskedImage.getMask().getPlaneBitMask("EDGE")
754 )
755 else:
756 replacementMaskedImage = rawMaskedImage
758 return replacementMaskedImage
761def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
762 """Apply bias correction in place.
764 Parameters
765 ----------
766 maskedImage : `lsst.afw.image.MaskedImage`
767 Image to process. The image is modified by this method.
768 biasMaskedImage : `lsst.afw.image.MaskedImage`
769 Bias image of the same size as ``maskedImage``
770 trimToFit : `Bool`, optional
771 If True, raw data is symmetrically trimmed to match
772 calibration size.
774 Raises
775 ------
776 RuntimeError
777 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
778 the same size.
780 """
781 if trimToFit:
782 maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
784 if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
785 raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
786 (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
787 maskedImage -= biasMaskedImage
790def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
791 """Apply dark correction in place.
793 Parameters
794 ----------
795 maskedImage : `lsst.afw.image.MaskedImage`
796 Image to process. The image is modified by this method.
797 darkMaskedImage : `lsst.afw.image.MaskedImage`
798 Dark image of the same size as ``maskedImage``.
799 expScale : scalar
800 Dark exposure time for ``maskedImage``.
801 darkScale : scalar
802 Dark exposure time for ``darkMaskedImage``.
803 invert : `Bool`, optional
804 If True, re-add the dark to an already corrected image.
805 trimToFit : `Bool`, optional
806 If True, raw data is symmetrically trimmed to match
807 calibration size.
809 Raises
810 ------
811 RuntimeError
812 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
813 the same size.
815 Notes
816 -----
817 The dark correction is applied by calculating:
818 maskedImage -= dark * expScaling / darkScaling
819 """
820 if trimToFit:
821 maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
823 if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
824 raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
825 (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
827 scale = expScale / darkScale
828 if not invert:
829 maskedImage.scaledMinus(scale, darkMaskedImage)
830 else:
831 maskedImage.scaledPlus(scale, darkMaskedImage)
834def updateVariance(maskedImage, gain, readNoise, replace=True):
835 """Set the variance plane based on the image plane.
837 The maskedImage must have units of `adu` (if gain != 1.0) or
838 electron (if gain == 1.0). This routine will always produce a
839 variance plane in the same units as the image.
841 Parameters
842 ----------
843 maskedImage : `lsst.afw.image.MaskedImage`
844 Image to process. The variance plane is modified.
845 gain : scalar
846 The amplifier gain in electron/adu.
847 readNoise : scalar
848 The amplifier read noise in electron/pixel.
849 replace : `bool`, optional
850 Replace the current variance? If False, the image
851 variance will be added to the current variance plane.
852 """
853 var = maskedImage.variance
854 if replace:
855 var[:, :] = maskedImage.image
856 else:
857 var[:, :] += maskedImage.image
858 var /= gain
859 var += (readNoise/gain)**2
862def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
863 """Apply flat correction in place.
865 Parameters
866 ----------
867 maskedImage : `lsst.afw.image.MaskedImage`
868 Image to process. The image is modified.
869 flatMaskedImage : `lsst.afw.image.MaskedImage`
870 Flat image of the same size as ``maskedImage``
871 scalingType : str
872 Flat scale computation method. Allowed values are 'MEAN',
873 'MEDIAN', or 'USER'.
874 userScale : scalar, optional
875 Scale to use if ``scalingType='USER'``.
876 invert : `Bool`, optional
877 If True, unflatten an already flattened image.
878 trimToFit : `Bool`, optional
879 If True, raw data is symmetrically trimmed to match
880 calibration size.
882 Raises
883 ------
884 RuntimeError
885 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
886 the same size or if ``scalingType`` is not an allowed value.
887 """
888 if trimToFit:
889 maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
891 if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
892 raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
893 (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
895 # Figure out scale from the data
896 # Ideally the flats are normalized by the calibration product pipeline,
897 # but this allows some flexibility in the case that the flat is created by
898 # some other mechanism.
899 if scalingType in ('MEAN', 'MEDIAN'):
900 scalingType = afwMath.stringToStatisticsProperty(scalingType)
901 flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
902 elif scalingType == 'USER':
903 flatScale = userScale
904 else:
905 raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
907 if not invert:
908 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
909 else:
910 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
913def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
914 """Apply illumination correction in place.
916 Parameters
917 ----------
918 maskedImage : `lsst.afw.image.MaskedImage`
919 Image to process. The image is modified.
920 illumMaskedImage : `lsst.afw.image.MaskedImage`
921 Illumination correction image of the same size as ``maskedImage``.
922 illumScale : scalar
923 Scale factor for the illumination correction.
924 trimToFit : `Bool`, optional
925 If True, raw data is symmetrically trimmed to match
926 calibration size.
928 Raises
929 ------
930 RuntimeError
931 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
932 the same size.
933 """
934 if trimToFit:
935 maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
937 if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
938 raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
939 (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
941 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
944def brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None):
945 """Apply brighter fatter correction in place for the image.
947 Parameters
948 ----------
949 exposure : `lsst.afw.image.Exposure`
950 Exposure to have brighter-fatter correction applied. Modified
951 by this method.
952 kernel : `numpy.ndarray`
953 Brighter-fatter kernel to apply.
954 maxIter : scalar
955 Number of correction iterations to run.
956 threshold : scalar
957 Convergence threshold in terms of the sum of absolute
958 deviations between an iteration and the previous one.
959 applyGain : `Bool`
960 If True, then the exposure values are scaled by the gain prior
961 to correction.
962 gains : `dict` [`str`, `float`]
963 A dictionary, keyed by amplifier name, of the gains to use.
964 If gains is None, the nominal gains in the amplifier object are used.
966 Returns
967 -------
968 diff : `float`
969 Final difference between iterations achieved in correction.
970 iteration : `int`
971 Number of iterations used to calculate correction.
973 Notes
974 -----
975 This correction takes a kernel that has been derived from flat
976 field images to redistribute the charge. The gradient of the
977 kernel is the deflection field due to the accumulated charge.
979 Given the original image I(x) and the kernel K(x) we can compute
980 the corrected image Ic(x) using the following equation:
982 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
984 To evaluate the derivative term we expand it as follows:
986 0.5 * ( d/dx(I(x))*d/dx(int(dy*K(x-y)*I(y)))
987 + I(x)*d^2/dx^2(int(dy* K(x-y)*I(y))) )
989 Because we use the measured counts instead of the incident counts
990 we apply the correction iteratively to reconstruct the original
991 counts and the correction. We stop iterating when the summed
992 difference between the current corrected image and the one from
993 the previous iteration is below the threshold. We do not require
994 convergence because the number of iterations is too large a
995 computational cost. How we define the threshold still needs to be
996 evaluated, the current default was shown to work reasonably well
997 on a small set of images. For more information on the method see
998 DocuShare Document-19407.
1000 The edges as defined by the kernel are not corrected because they
1001 have spurious values due to the convolution.
1002 """
1003 image = exposure.getMaskedImage().getImage()
1005 # The image needs to be units of electrons/holes
1006 with gainContext(exposure, image, applyGain, gains):
1008 kLx = numpy.shape(kernel)[0]
1009 kLy = numpy.shape(kernel)[1]
1010 kernelImage = afwImage.ImageD(kLx, kLy)
1011 kernelImage.getArray()[:, :] = kernel
1012 tempImage = afwImage.ImageD(image, deep=True)
1014 nanIndex = numpy.isnan(tempImage.getArray())
1015 tempImage.getArray()[nanIndex] = 0.
1017 corr = numpy.zeros(image.array.shape, dtype=numpy.float64)
1018 prev_image = numpy.zeros(image.array.shape, dtype=numpy.float64)
1020 # Define boundary by convolution region. The region that the
1021 # correction will be calculated for is one fewer in each dimension
1022 # because of the second derivative terms.
1023 # NOTE: these need to use integer math, as we're using start:end as
1024 # numpy index ranges.
1025 startX = kLx//2
1026 endX = -kLx//2
1027 startY = kLy//2
1028 endY = -kLy//2
1030 for iteration in range(maxIter):
1032 outArray = scipy.signal.convolve(
1033 tempImage.array,
1034 kernelImage.array,
1035 mode="same",
1036 method="fft",
1037 )
1038 tmpArray = tempImage.getArray()
1040 with numpy.errstate(invalid="ignore", over="ignore"):
1041 # First derivative term
1042 gradTmp = numpy.gradient(tmpArray[startY:endY, startX:endX])
1043 gradOut = numpy.gradient(outArray[startY:endY, startX:endX])
1044 first = (gradTmp[0]*gradOut[0] + gradTmp[1]*gradOut[1])[1:-1, 1:-1]
1046 # Second derivative term
1047 diffOut20 = numpy.diff(outArray, 2, 0)[startY:endY, startX + 1:endX - 1]
1048 diffOut21 = numpy.diff(outArray, 2, 1)[startY + 1:endY - 1, startX:endX]
1049 second = tmpArray[startY + 1:endY - 1, startX + 1:endX - 1]*(diffOut20 + diffOut21)
1051 corr[startY + 1:endY - 1, startX + 1:endX - 1] = 0.5*(first + second)
1053 tmpArray[:, :] = image.getArray()[:, :]
1054 tmpArray[nanIndex] = 0.
1055 tmpArray[startY:endY, startX:endX] += corr[startY:endY, startX:endX]
1057 if iteration > 0:
1058 diff = numpy.sum(numpy.abs(prev_image - tmpArray), dtype=numpy.float64)
1060 if diff < threshold:
1061 break
1062 prev_image[:, :] = tmpArray[:, :]
1064 image.getArray()[startY + 1:endY - 1, startX + 1:endX - 1] += \
1065 corr[startY + 1:endY - 1, startX + 1:endX - 1]
1067 return diff, iteration
1070def transferFlux(cFunc, fStep, correctionMode=True):
1071 """Take the input convolved deflection potential and the flux array
1072 to compute and apply the flux transfer into the correction array.
1074 Parameters
1075 ----------
1076 cFunc: `numpy.array`
1077 Deflection potential, being the convolution of the flux F with the
1078 kernel K.
1079 fStep: `numpy.array`
1080 The array of flux values which act as the source of the flux transfer.
1081 correctionMode: `bool`
1082 Defines if applying correction (True) or generating sims (False).
1084 Returns
1085 -------
1086 corr:
1087 BFE correction array
1088 """
1090 if cFunc.shape != fStep.shape:
1091 raise RuntimeError(f'transferFlux: array shapes do not match: {cFunc.shape}, {fStep.shape}')
1093 # set the sign of the correction and set its value for the
1094 # time averaged solution
1095 if correctionMode:
1096 # negative sign if applying BFE correction
1097 factor = -0.5
1098 else:
1099 # positive sign if generating BFE simulations
1100 factor = 0.5
1102 # initialise the BFE correction image to zero
1103 corr = numpy.zeros(cFunc.shape, dtype=numpy.float64)
1105 # Generate a 2D mesh of x,y coordinates
1106 yDim, xDim = cFunc.shape
1107 y = numpy.arange(yDim, dtype=int)
1108 x = numpy.arange(xDim, dtype=int)
1109 xc, yc = numpy.meshgrid(x, y)
1111 # process each axis in turn
1112 for ax in [0, 1]:
1114 # gradient of phi on right/upper edge of pixel
1115 diff = numpy.diff(cFunc, axis=ax)
1117 # expand array back to full size with zero gradient at the end
1118 gx = numpy.zeros(cFunc.shape, dtype=numpy.float64)
1119 yDiff, xDiff = diff.shape
1120 gx[:yDiff, :xDiff] += diff
1122 # select pixels with either positive gradients on the right edge,
1123 # flux flowing to the right/up
1124 # or negative gradients, flux flowing to the left/down
1125 for i, sel in enumerate([gx > 0, gx < 0]):
1126 xSelPixels = xc[sel]
1127 ySelPixels = yc[sel]
1128 # and add the flux into the pixel to the right or top
1129 # depending on which axis we are handling
1130 if ax == 0:
1131 xPix = xSelPixels
1132 yPix = ySelPixels+1
1133 else:
1134 xPix = xSelPixels+1
1135 yPix = ySelPixels
1136 # define flux as the either current pixel value or pixel
1137 # above/right
1138 # depending on whether positive or negative gradient
1139 if i == 0:
1140 # positive gradients, flux flowing to higher coordinate values
1141 flux = factor * fStep[sel]*gx[sel]
1142 else:
1143 # negative gradients, flux flowing to lower coordinate values
1144 flux = factor * fStep[yPix, xPix]*gx[sel]
1145 # change the fluxes of the donor and receiving pixels
1146 # such that flux is conserved
1147 corr[sel] -= flux
1148 corr[yPix, xPix] += flux
1150 # return correction array
1151 return corr
1154def fluxConservingBrighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain,
1155 gains=None, correctionMode=True):
1156 """Apply brighter fatter correction in place for the image.
1158 This version presents a modified version of the algorithm
1159 found in ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``
1160 which conserves the image flux, resulting in improved
1161 correction of the cores of stars. The convolution has also been
1162 modified to mitigate edge effects.
1164 Parameters
1165 ----------
1166 exposure : `lsst.afw.image.Exposure`
1167 Exposure to have brighter-fatter correction applied. Modified
1168 by this method.
1169 kernel : `numpy.ndarray`
1170 Brighter-fatter kernel to apply.
1171 maxIter : scalar
1172 Number of correction iterations to run.
1173 threshold : scalar
1174 Convergence threshold in terms of the sum of absolute
1175 deviations between an iteration and the previous one.
1176 applyGain : `Bool`
1177 If True, then the exposure values are scaled by the gain prior
1178 to correction.
1179 gains : `dict` [`str`, `float`]
1180 A dictionary, keyed by amplifier name, of the gains to use.
1181 If gains is None, the nominal gains in the amplifier object are used.
1182 correctionMode : `Bool`
1183 If True (default) the function applies correction for BFE. If False,
1184 the code can instead be used to generate a simulation of BFE (sign
1185 change in the direction of the effect)
1187 Returns
1188 -------
1189 diff : `float`
1190 Final difference between iterations achieved in correction.
1191 iteration : `int`
1192 Number of iterations used to calculate correction.
1194 Notes
1195 -----
1196 Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.
1198 This correction takes a kernel that has been derived from flat
1199 field images to redistribute the charge. The gradient of the
1200 kernel is the deflection field due to the accumulated charge.
1202 Given the original image I(x) and the kernel K(x) we can compute
1203 the corrected image Ic(x) using the following equation:
1205 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
1207 Improved algorithm at this step applies the divergence theorem to
1208 obtain a pixelised correction.
1210 Because we use the measured counts instead of the incident counts
1211 we apply the correction iteratively to reconstruct the original
1212 counts and the correction. We stop iterating when the summed
1213 difference between the current corrected image and the one from
1214 the previous iteration is below the threshold. We do not require
1215 convergence because the number of iterations is too large a
1216 computational cost. How we define the threshold still needs to be
1217 evaluated, the current default was shown to work reasonably well
1218 on a small set of images.
1220 Edges are handled in the convolution by padding. This is still not
1221 a physical model for the edge, but avoids discontinuity in the correction.
1223 Author of modified version: Lance.Miller@physics.ox.ac.uk
1224 (see DM-38555).
1225 """
1226 image = exposure.getMaskedImage().getImage()
1228 # The image needs to be units of electrons/holes
1229 with gainContext(exposure, image, applyGain, gains):
1231 # get kernel and its shape
1232 kLy, kLx = kernel.shape
1233 kernelImage = afwImage.ImageD(kLx, kLy)
1234 kernelImage.getArray()[:, :] = kernel
1235 tempImage = afwImage.ImageD(image, deep=True)
1237 nanIndex = numpy.isnan(tempImage.getArray())
1238 tempImage.getArray()[nanIndex] = 0.
1240 outImage = afwImage.ImageD(image.getDimensions())
1241 corr = numpy.zeros(image.array.shape, dtype=numpy.float64)
1242 prevImage = numpy.zeros(image.array.shape, dtype=numpy.float64)
1243 convCntrl = afwMath.ConvolutionControl(False, False, 1)
1244 fixedKernel = afwMath.FixedKernel(kernelImage)
1246 # set the padding amount
1247 # ensure we pad by an even amount larger than the kernel
1248 kLy = 2 * ((1+kLy)//2)
1249 kLx = 2 * ((1+kLx)//2)
1251 # The deflection potential only depends on the gradient of
1252 # the convolution, so we can subtract the mean, which then
1253 # allows us to pad the image with zeros and avoid wrap-around effects
1254 # (although still not handling the image edges with a physical model)
1255 # This wouldn't be great if there were a strong image gradient.
1256 imYdimension, imXdimension = tempImage.array.shape
1257 imean = numpy.mean(tempImage.getArray()[~nanIndex], dtype=numpy.float64)
1258 # subtract mean from image
1259 tempImage -= imean
1260 tempImage.array[nanIndex] = 0.0
1261 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
1262 outImage = afwImage.ImageD(numpy.pad(outImage.getArray(), ((0, kLy), (0, kLx))))
1263 # Convert array to afw image so afwMath.convolve works
1264 padImage = afwImage.ImageD(padArray.shape[1], padArray.shape[0])
1265 padImage.array[:] = padArray
1267 for iteration in range(maxIter):
1269 # create deflection potential, convolution of flux with kernel
1270 # using padded counts array
1271 afwMath.convolve(outImage, padImage, fixedKernel, convCntrl)
1272 tmpArray = tempImage.getArray()
1273 outArray = outImage.getArray()
1275 # trim convolution output back to original shape
1276 outArray = outArray[:imYdimension, :imXdimension]
1278 # generate the correction array, with correctionMode set as input
1279 corr[...] = transferFlux(outArray, tmpArray, correctionMode=correctionMode)
1281 # update the arrays for the next iteration
1282 tmpArray[:, :] = image.getArray()[:, :]
1283 tmpArray += corr
1284 tmpArray[nanIndex] = 0.
1285 # update padded array
1286 # subtract mean
1287 tmpArray -= imean
1288 tempImage.array[nanIndex] = 0.
1289 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
1291 if iteration > 0:
1292 diff = numpy.sum(numpy.abs(prevImage - tmpArray), dtype=numpy.float64)
1294 if diff < threshold:
1295 break
1296 prevImage[:, :] = tmpArray[:, :]
1298 image.getArray()[:] += corr[:]
1300 return diff, iteration
1303@contextmanager
1304def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
1305 """Context manager that applies and removes gain.
1307 Parameters
1308 ----------
1309 exp : `lsst.afw.image.Exposure`
1310 Exposure to apply/remove gain.
1311 image : `lsst.afw.image.Image`
1312 Image to apply/remove gain.
1313 apply : `bool`
1314 If True, apply and remove the amplifier gain.
1315 gains : `dict` [`str`, `float`], optional
1316 A dictionary, keyed by amplifier name, of the gains to use.
1317 If gains is None, the nominal gains in the amplifier object are used.
1318 invert : `bool`, optional
1319 Invert the gains (e.g. convert electrons to adu temporarily)?
1320 isTrimmed : `bool`, optional
1321 Is this a trimmed exposure?
1323 Yields
1324 ------
1325 exp : `lsst.afw.image.Exposure`
1326 Exposure with the gain applied.
1327 """
1328 # check we have all of them if provided because mixing and matching would
1329 # be a real mess
1330 if gains and apply is True:
1331 ampNames = [amp.getName() for amp in exp.getDetector()]
1332 for ampName in ampNames:
1333 if ampName not in gains.keys():
1334 raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
1336 if apply:
1337 ccd = exp.getDetector()
1338 for amp in ccd:
1339 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1340 if gains:
1341 gain = gains[amp.getName()]
1342 else:
1343 gain = amp.getGain()
1344 if invert:
1345 sim /= gain
1346 else:
1347 sim *= gain
1349 try:
1350 yield exp
1351 finally:
1352 if apply:
1353 ccd = exp.getDetector()
1354 for amp in ccd:
1355 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1356 if gains:
1357 gain = gains[amp.getName()]
1358 else:
1359 gain = amp.getGain()
1360 if invert:
1361 sim *= gain
1362 else:
1363 sim /= gain
1366def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
1367 sensorTransmission=None, atmosphereTransmission=None):
1368 """Attach a TransmissionCurve to an Exposure, given separate curves for
1369 different components.
1371 Parameters
1372 ----------
1373 exposure : `lsst.afw.image.Exposure`
1374 Exposure object to modify by attaching the product of all given
1375 ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
1376 Must have a valid ``Detector`` attached that matches the detector
1377 associated with sensorTransmission.
1378 opticsTransmission : `lsst.afw.image.TransmissionCurve`
1379 A ``TransmissionCurve`` that represents the throughput of the optics,
1380 to be evaluated in focal-plane coordinates.
1381 filterTransmission : `lsst.afw.image.TransmissionCurve`
1382 A ``TransmissionCurve`` that represents the throughput of the filter
1383 itself, to be evaluated in focal-plane coordinates.
1384 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1385 A ``TransmissionCurve`` that represents the throughput of the sensor
1386 itself, to be evaluated in post-assembly trimmed detector coordinates.
1387 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1388 A ``TransmissionCurve`` that represents the throughput of the
1389 atmosphere, assumed to be spatially constant.
1391 Returns
1392 -------
1393 combined : `lsst.afw.image.TransmissionCurve`
1394 The TransmissionCurve attached to the exposure.
1396 Notes
1397 -----
1398 All ``TransmissionCurve`` arguments are optional; if none are provided, the
1399 attached ``TransmissionCurve`` will have unit transmission everywhere.
1400 """
1401 combined = afwImage.TransmissionCurve.makeIdentity()
1402 if atmosphereTransmission is not None:
1403 combined *= atmosphereTransmission
1404 if opticsTransmission is not None:
1405 combined *= opticsTransmission
1406 if filterTransmission is not None:
1407 combined *= filterTransmission
1408 detector = exposure.getDetector()
1409 fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
1410 toSys=camGeom.PIXELS)
1411 combined = combined.transformedBy(fpToPix)
1412 if sensorTransmission is not None:
1413 combined *= sensorTransmission
1414 exposure.getInfo().setTransmissionCurve(combined)
1415 return combined
1418def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
1419 """Scale an exposure by the amplifier gains.
1421 Parameters
1422 ----------
1423 exposure : `lsst.afw.image.Exposure`
1424 Exposure to process. The image is modified.
1425 normalizeGains : `Bool`, optional
1426 If True, then amplifiers are scaled to force the median of
1427 each amplifier to equal the median of those medians.
1428 ptcGains : `dict`[`str`], optional
1429 Dictionary keyed by amp name containing the PTC gains.
1430 isTrimmed : `bool`, optional
1431 Is the input image trimmed?
1432 """
1433 ccd = exposure.getDetector()
1434 ccdImage = exposure.getMaskedImage()
1436 medians = []
1437 for amp in ccd:
1438 if isTrimmed:
1439 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1440 else:
1441 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1442 if ptcGains:
1443 sim *= ptcGains[amp.getName()]
1444 else:
1445 sim *= amp.getGain()
1447 if normalizeGains:
1448 medians.append(numpy.median(sim.getImage().getArray()))
1450 if normalizeGains:
1451 median = numpy.median(numpy.array(medians))
1452 for index, amp in enumerate(ccd):
1453 if isTrimmed:
1454 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1455 else:
1456 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1457 if medians[index] != 0.0:
1458 sim *= median/medians[index]
1461def widenSaturationTrails(mask):
1462 """Grow the saturation trails by an amount dependent on the width of the
1463 trail.
1465 Parameters
1466 ----------
1467 mask : `lsst.afw.image.Mask`
1468 Mask which will have the saturated areas grown.
1469 """
1471 extraGrowDict = {}
1472 for i in range(1, 6):
1473 extraGrowDict[i] = 0
1474 for i in range(6, 8):
1475 extraGrowDict[i] = 1
1476 for i in range(8, 10):
1477 extraGrowDict[i] = 3
1478 extraGrowMax = 4
1480 if extraGrowMax <= 0:
1481 return
1483 saturatedBit = mask.getPlaneBitMask("SAT")
1485 xmin, ymin = mask.getBBox().getMin()
1486 width = mask.getWidth()
1488 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1489 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1491 for fp in fpList:
1492 for s in fp.getSpans():
1493 x0, x1 = s.getX0(), s.getX1()
1495 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1496 if extraGrow > 0:
1497 y = s.getY() - ymin
1498 x0 -= xmin + extraGrow
1499 x1 -= xmin - extraGrow
1501 if x0 < 0:
1502 x0 = 0
1503 if x1 >= width - 1:
1504 x1 = width - 1
1506 mask.array[y, x0:x1+1] |= saturatedBit
1509def setBadRegions(exposure, badStatistic="MEDIAN"):
1510 """Set all BAD areas of the chip to the average of the rest of the exposure
1512 Parameters
1513 ----------
1514 exposure : `lsst.afw.image.Exposure`
1515 Exposure to mask. The exposure mask is modified.
1516 badStatistic : `str`, optional
1517 Statistic to use to generate the replacement value from the
1518 image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
1520 Returns
1521 -------
1522 badPixelCount : scalar
1523 Number of bad pixels masked.
1524 badPixelValue : scalar
1525 Value substituted for bad pixels.
1527 Raises
1528 ------
1529 RuntimeError
1530 Raised if `badStatistic` is not an allowed value.
1531 """
1532 if badStatistic == "MEDIAN":
1533 statistic = afwMath.MEDIAN
1534 elif badStatistic == "MEANCLIP":
1535 statistic = afwMath.MEANCLIP
1536 else:
1537 raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
1539 mi = exposure.getMaskedImage()
1540 mask = mi.getMask()
1541 BAD = mask.getPlaneBitMask("BAD")
1542 INTRP = mask.getPlaneBitMask("INTRP")
1544 sctrl = afwMath.StatisticsControl()
1545 sctrl.setAndMask(BAD)
1546 value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
1548 maskArray = mask.getArray()
1549 imageArray = mi.getImage().getArray()
1550 badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
1551 imageArray[:] = numpy.where(badPixels, value, imageArray)
1553 return badPixels.sum(), value
1556def checkFilter(exposure, filterList, log):
1557 """Check to see if an exposure is in a filter specified by a list.
1559 The goal of this is to provide a unified filter checking interface
1560 for all filter dependent stages.
1562 Parameters
1563 ----------
1564 exposure : `lsst.afw.image.Exposure`
1565 Exposure to examine.
1566 filterList : `list` [`str`]
1567 List of physical_filter names to check.
1568 log : `logging.Logger`
1569 Logger to handle messages.
1571 Returns
1572 -------
1573 result : `bool`
1574 True if the exposure's filter is contained in the list.
1575 """
1576 if len(filterList) == 0:
1577 return False
1578 thisFilter = exposure.getFilter()
1579 if thisFilter is None:
1580 log.warning("No FilterLabel attached to this exposure!")
1581 return False
1583 thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
1584 if thisPhysicalFilter in filterList:
1585 return True
1586 elif thisFilter.bandLabel in filterList:
1587 if log:
1588 log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
1589 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1590 return True
1591 else:
1592 return False
1595def getPhysicalFilter(filterLabel, log):
1596 """Get the physical filter label associated with the given filterLabel.
1598 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1599 associated with the given ``filterLabel``, the returned label will be
1600 "Unknown".
1602 Parameters
1603 ----------
1604 filterLabel : `lsst.afw.image.FilterLabel`
1605 The `lsst.afw.image.FilterLabel` object from which to derive the
1606 physical filter label.
1607 log : `logging.Logger`
1608 Logger to handle messages.
1610 Returns
1611 -------
1612 physicalFilter : `str`
1613 The value returned by the physicalLabel attribute of ``filterLabel`` if
1614 it exists, otherwise set to \"Unknown\".
1615 """
1616 if filterLabel is None:
1617 physicalFilter = "Unknown"
1618 log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
1619 else:
1620 try:
1621 physicalFilter = filterLabel.physicalLabel
1622 except RuntimeError:
1623 log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
1624 physicalFilter = "Unknown"
1625 return physicalFilter
1628def countMaskedPixels(maskedIm, maskPlane):
1629 """Count the number of pixels in a given mask plane.
1631 Parameters
1632 ----------
1633 maskedIm : `~lsst.afw.image.MaskedImage`
1634 Masked image to examine.
1635 maskPlane : `str`
1636 Name of the mask plane to examine.
1638 Returns
1639 -------
1640 nPix : `int`
1641 Number of pixels in the requested mask plane.
1642 """
1643 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1644 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1645 return nPix
1648def getExposureGains(exposure):
1649 """Get the per-amplifier gains used for this exposure.
1651 Parameters
1652 ----------
1653 exposure : `lsst.afw.image.Exposure`
1654 The exposure to find gains for.
1656 Returns
1657 -------
1658 gains : `dict` [`str` `float`]
1659 Dictionary of gain values, keyed by amplifier name.
1660 Returns empty dict when detector is None.
1661 """
1662 det = exposure.getDetector()
1663 if det is None:
1664 return dict()
1666 metadata = exposure.getMetadata()
1667 gains = {}
1668 for amp in det:
1669 ampName = amp.getName()
1670 # The key may use the new LSST ISR or the old LSST prefix
1671 if (key1 := f"LSST ISR GAIN {ampName}") in metadata:
1672 gains[ampName] = metadata[key1]
1673 elif (key2 := f"LSST GAIN {ampName}") in metadata:
1674 gains[ampName] = metadata[key2]
1675 else:
1676 gains[ampName] = amp.getGain()
1677 return gains
1680def getExposureReadNoises(exposure):
1681 """Get the per-amplifier read noise used for this exposure.
1683 Parameters
1684 ----------
1685 exposure : `lsst.afw.image.Exposure`
1686 The exposure to find read noise for.
1688 Returns
1689 -------
1690 readnoises : `dict` [`str` `float`]
1691 Dictionary of read noise values, keyed by amplifier name.
1692 Returns empty dict when detector is None.
1693 """
1694 det = exposure.getDetector()
1695 if det is None:
1696 return dict()
1698 metadata = exposure.getMetadata()
1699 readnoises = {}
1700 for amp in det:
1701 ampName = amp.getName()
1702 # The key may use the new LSST ISR or the old LSST prefix
1703 if (key1 := f"LSST ISR READNOISE {ampName}") in metadata:
1704 readnoises[ampName] = metadata[key1]
1705 elif (key2 := f"LSST READNOISE {ampName}") in metadata:
1706 readnoises[ampName] = metadata[key2]
1707 else:
1708 readnoises[ampName] = amp.getReadNoise()
1709 return readnoises
1712def isTrimmedExposure(exposure):
1713 """Check if the unused pixels (pre-/over-scan pixels) have
1714 been trimmed from an exposure.
1716 Parameters
1717 ----------
1718 exposure : `lsst.afw.image.Exposure`
1719 The exposure to check.
1721 Returns
1722 -------
1723 result : `bool`
1724 True if the image is trimmed, else False.
1725 """
1726 return exposure.getDetector().getBBox() == exposure.getBBox()
1729def isTrimmedImage(image, detector):
1730 """Check if the unused pixels (pre-/over-scan pixels) have
1731 been trimmed from an image
1733 Parameters
1734 ----------
1735 image : `lsst.afw.image.Image`
1736 The image to check.
1737 detector : `lsst.afw.cameraGeom.Detector`
1738 The detector associated with the image.
1740 Returns
1741 -------
1742 result : `bool`
1743 True if the image is trimmed, else False.
1744 """
1745 return detector.getBBox() == image.getBBox()
1748def compareCameraKeywords(
1749 doRaiseOnCalibMismatch,
1750 cameraKeywordsToCompare,
1751 exposureMetadata,
1752 calib,
1753 calibName,
1754 log=None,
1755):
1756 """Compare header keywords to confirm camera states match.
1758 Parameters
1759 ----------
1760 doRaiseOnCalibMismatch : `bool`
1761 Raise on calibration mismatch? Otherwise, log a warning.
1762 cameraKeywordsToCompare : `list` [`str`]
1763 List of camera keywords to compare.
1764 exposureMetadata : `lsst.daf.base.PropertyList`
1765 Header for the exposure being processed.
1766 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1767 Calibration to be applied.
1768 calibName : `str`
1769 Calib type for log message.
1770 log : `logging.Logger`, optional
1771 Logger to handle messages.
1772 """
1773 try:
1774 calibMetadata = calib.metadata
1775 except AttributeError:
1776 return
1778 log = log if log else logging.getLogger(__name__)
1780 missingKeywords = []
1781 for keyword in cameraKeywordsToCompare:
1782 exposureValue = exposureMetadata.get(keyword, None)
1783 if exposureValue is None:
1784 log.debug("Sequencer keyword %s not found in exposure metadata.", keyword)
1785 continue
1787 calibValue = calibMetadata.get(keyword, None)
1789 # We don't log here if there is a missing keyword.
1790 if calibValue is None:
1791 missingKeywords.append(keyword)
1792 continue
1794 if exposureValue != calibValue:
1795 if doRaiseOnCalibMismatch:
1796 raise RuntimeError(
1797 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1798 calibName,
1799 keyword,
1800 exposureValue,
1801 calibValue,
1802 )
1803 else:
1804 log.warning(
1805 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1806 calibName,
1807 keyword,
1808 exposureValue,
1809 calibValue,
1810 )
1811 exposureMetadata[f"ISR {calibName.upper()} SEQUENCER MISMATCH"] = True
1813 if missingKeywords:
1814 log.info(
1815 "Calibration %s missing keywords %s, which were not checked.",
1816 calibName,
1817 ",".join(missingKeywords),
1818 )