Loading [MathJax]/extensions/tex2jax.js
LSST Applications g032c94a9f9+a1301e4c20,g04a91732dc+321623803d,g07dc498a13+dc60e07e33,g0fba68d861+b7e4830700,g1409bbee79+dc60e07e33,g1a7e361dbc+dc60e07e33,g1fd858c14a+bc317df4c0,g208c678f98+64d2817f4c,g2c84ff76c0+9484f2668e,g35bb328faa+fcb1d3bbc8,g4d2262a081+fb060387ce,g4d39ba7253+0f38e7b1d1,g4e0f332c67+5d362be553,g53246c7159+fcb1d3bbc8,g60b5630c4e+0f38e7b1d1,g78460c75b0+2f9a1b4bcd,g786e29fd12+cf7ec2a62a,g7b71ed6315+fcb1d3bbc8,g8852436030+1ad2ae6bba,g89139ef638+dc60e07e33,g8d6b6b353c+0f38e7b1d1,g9125e01d80+fcb1d3bbc8,g989de1cb63+dc60e07e33,g9f33ca652e+602c5da793,ga9baa6287d+0f38e7b1d1,gaaedd4e678+dc60e07e33,gabe3b4be73+1e0a283bba,gb1101e3267+4e433ac613,gb58c049af0+f03b321e39,gb90eeb9370+6b7d01c6c0,gcf25f946ba+1ad2ae6bba,gd315a588df+cb74d54ad7,gd6cbbdb0b4+c8606af20c,gd9a9a58781+fcb1d3bbc8,gde0f65d7ad+93c67b85fe,ge278dab8ac+932305ba37,ge82c20c137+76d20ab76d,gf18def8413+4ce00804e3,w.2025.11
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
isrFunctions.py
Go to the documentation of this file.
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#
22
23__all__ = [
24 "applyGains",
25 "attachTransmissionCurve",
26 "biasCorrection",
27 "brighterFatterCorrection",
28 "checkFilter",
29 "countMaskedPixels",
30 "createPsf",
31 "darkCorrection",
32 "flatCorrection",
33 "fluxConservingBrighterFatterCorrection",
34 "gainContext",
35 "getPhysicalFilter",
36 "growMasks",
37 "maskITLEdgeBleed",
38 "illuminationCorrection",
39 "interpolateDefectList",
40 "interpolateFromMask",
41 "makeThresholdMask",
42 "saturationCorrection",
43 "setBadRegions",
44 "transferFlux",
45 "transposeMaskedImage",
46 "trimToMatchCalibBBox",
47 "updateVariance",
48 "widenSaturationTrails",
49 "getExposureGains",
50 "getExposureReadNoises",
51]
52
53import math
54import numpy
55
56import lsst.geom
57import lsst.afw.image as afwImage
58import lsst.afw.detection as afwDetection
59import lsst.afw.math as afwMath
60import lsst.meas.algorithms as measAlg
61import lsst.afw.cameraGeom as camGeom
62
63from lsst.meas.algorithms.detection import SourceDetectionTask
64
65from contextlib import contextmanager
66
67from .defects import Defects
68
69
70def createPsf(fwhm):
71 """Make a double Gaussian PSF.
72
73 Parameters
74 ----------
75 fwhm : scalar
76 FWHM of double Gaussian smoothing kernel.
77
78 Returns
79 -------
80 psf : `lsst.meas.algorithms.DoubleGaussianPsf`
81 The created smoothing kernel.
82 """
83 ksize = 4*int(fwhm) + 1
84 return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
85
86
87def transposeMaskedImage(maskedImage):
88 """Make a transposed copy of a masked image.
89
90 Parameters
91 ----------
92 maskedImage : `lsst.afw.image.MaskedImage`
93 Image to process.
94
95 Returns
96 -------
97 transposed : `lsst.afw.image.MaskedImage`
98 The transposed copy of the input image.
99 """
100 transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
101 transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
102 transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
103 transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
104 return transposed
105
106
107def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None,
108 maskNameList=None, useLegacyInterp=True):
109 """Interpolate over defects specified in a defect list.
110
111 Parameters
112 ----------
113 maskedImage : `lsst.afw.image.MaskedImage`
114 Image to process.
115 defectList : `lsst.meas.algorithms.Defects`
116 List of defects to interpolate over.
117 fwhm : `float`
118 FWHM of double Gaussian smoothing kernel.
119 fallbackValue : scalar, optional
120 Fallback value if an interpolated value cannot be determined.
121 If None, then the clipped mean of the image is used.
122 maskNameList : `list [string]`
123 List of the defects to interpolate over (used for GP interpolator).
124 useLegacyInterp : `bool`
125 Use the legacy interpolation (polynomial interpolation) if True. Use
126 Gaussian Process interpolation if False.
127
128 Notes
129 -----
130 The ``fwhm`` parameter is used to create a PSF, but the underlying
131 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
132 not currently make use of this information in legacy Interpolation, but use
133 if for the Gaussian Process as an estimation of the correlation lenght.
134 """
135 psf = createPsf(fwhm)
136 if fallbackValue is None:
137 fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
138 if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
139 maskedImage.getMask().addMaskPlane('INTRP')
140
141 # Hardcoded fwhm value. PSF estimated latter in step1,
142 # not in ISR.
143 if useLegacyInterp:
144 kwargs = {}
145 fwhm = fwhm
146 else:
147 # tested on a dozens of images and looks a good set of
148 # hyperparameters, but cannot guarrenty this is optimal,
149 # need further testing.
150 kwargs = {"bin_spacing": 20,
151 "threshold_dynamic_binning": 2000,
152 "threshold_subdivide": 20000}
153 fwhm = 15
154
155 measAlg.interpolateOverDefects(maskedImage, psf, defectList,
156 fallbackValue=fallbackValue,
157 useFallbackValueAtEdge=True,
158 fwhm=fwhm,
159 useLegacyInterp=useLegacyInterp,
160 maskNameList=maskNameList, **kwargs)
161 return maskedImage
162
163
164def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
165 """Mask pixels based on threshold detection.
166
167 Parameters
168 ----------
169 maskedImage : `lsst.afw.image.MaskedImage`
170 Image to process. Only the mask plane is updated.
171 threshold : scalar
172 Detection threshold.
173 growFootprints : scalar, optional
174 Number of pixels to grow footprints of detected regions.
175 maskName : str, optional
176 Mask plane name, or list of names to convert
177
178 Returns
179 -------
180 defectList : `lsst.meas.algorithms.Defects`
181 Defect list constructed from pixels above the threshold.
182 """
183 # find saturated regions
184 thresh = afwDetection.Threshold(threshold)
185 fs = afwDetection.FootprintSet(maskedImage, thresh)
186
187 if growFootprints > 0:
188 fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
189 fpList = fs.getFootprints()
190
191 # set mask
192 mask = maskedImage.getMask()
193 bitmask = mask.getPlaneBitMask(maskName)
194 afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
195
196 return Defects.fromFootprintList(fpList)
197
198
199def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
200 """Grow a mask by an amount and add to the requested plane.
201
202 Parameters
203 ----------
204 mask : `lsst.afw.image.Mask`
205 Mask image to process.
206 radius : scalar
207 Amount to grow the mask.
208 maskNameList : `str` or `list` [`str`]
209 Mask names that should be grown.
210 maskValue : `str`
211 Mask plane to assign the newly masked pixels to.
212 """
213 if radius > 0:
214 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
215 fpSet = afwDetection.FootprintSet(mask, thresh)
216 fpSet = afwDetection.FootprintSet(fpSet, rGrow=radius, isotropic=False)
217 fpSet.setMask(mask, maskValue)
218
219
220def maskITLEdgeBleed(ccdExposure, badAmpDict, itlEdgeBleedSatMinArea=10000,
221 itlEdgeBleedSatMaxArea=100000,
222 itlEdgeBleedThreshold=5000.,
223 itlEdgeBleedModelConstant=0.03,
224 saturatedMaskName="SAT"):
225 """Mask edge bleeds in ITL detectors.
226
227 Parameters
228 ----------
229 ccdExposure : `lsst.afw.image.Exposure`
230 Exposure to apply masking to.
231 badAmpDict : `dict` [`str`, `bool`]
232 Dictionary of amplifiers, keyed by name, value is True if
233 amplifier is fully masked.
234 itlEdgeBleedSatMinArea : `int`, optional
235 Minimal saturated footprint area where the presence of edge bleeds
236 will be checked.
237 itlEdgeBleedThreshold : `float`, optional
238 Threshold above median sky background for edge bleed detection
239 (electron units).
240 itlEdgeBleedModelConstant : `float`, optional
241 Constant in the decaying exponential in the edge bleed masking.
242 saturatedMaskName : `str`, optional
243 Mask name for saturation.
244 """
245 maskedImage = ccdExposure.maskedImage
246 xmax = maskedImage.image.array.shape[1]
247
248 # Get minimum of amplifier saturation level
249 satLevel = numpy.min([ccdExposure.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
250 for amp in ccdExposure.getDetector() if not badAmpDict[amp.getName()]])
251
252 # Get footprint list of saturated pixels
253 thresh = afwDetection.Threshold(satLevel)
254 fpList = afwDetection.FootprintSet(maskedImage, thresh).getFootprints()
255
256 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
257 largeAreas, = numpy.where((satAreas >= itlEdgeBleedSatMinArea)
258 & (satAreas < itlEdgeBleedSatMaxArea))
259
260 satMaskBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
261
262 if len(largeAreas) == 0:
263 return
264
265 for largeAreasIndex in largeAreas:
266
267 # Get centroid of saturated core
268 xCore, yCore = fpList[largeAreasIndex].getCentroid()
269 # Turn the Y detector coordinate
270 # into Y footprint coordinate
271 yCoreFP = yCore - fpList[largeAreasIndex].getBBox().getMinY()
272 # Get the number of saturated columns around the centroid
273 widthSat = numpy.sum(fpList[largeAreasIndex].getSpans().asArray()[int(yCoreFP), :])
274
275 for amp in ccdExposure.getDetector():
276 # Select the 2 top and bottom amplifiers around the saturated
277 # core with a potential edge bleed by selecting the amplifiers
278 # that have the same X coordinate as the saturated core.
279 # As we don't care about the Y coordinate, we set it to the
280 # center of the BBox.
281 yBox = amp.getBBox().getCenter()[1]
282 if amp.getBBox().contains(xCore, yBox):
283
284 # Get the amp name
285 ampName = amp.getName()
286
287 # Because in ITLs the edge bleed happens on both edges
288 # of the detector, we make a cutout around
289 # both the top and bottom
290 # edge bleed candidates around the saturated core.
291 # We flip the cutout of the top amplifier
292 # to then work with the same coordinates for both.
293 # The way of selecting top vs bottom amp
294 # is very specific to ITL.
295 if ampName[:2] == 'C1':
296 sliceImage = maskedImage.image.array[:200, :]
297 sliceMask = maskedImage.mask.array[:200, :]
298 elif ampName[:2] == 'C0':
299 sliceImage = numpy.flipud(maskedImage.image.array[-200:, :])
300 sliceMask = numpy.flipud(maskedImage.mask.array[-200:, :])
301
302 # The middle columns of edge bleeds often have
303 # high counts, so we check there is an edge bleed
304 # by looking at a small image up to 50 pixels from the edge
305 # and around the saturated columns
306 # of the saturated core, and checking its median is
307 # above the sky background by itlEdgeBleedThreshold
308
309 # If the centroid is too close to the edge of the detector
310 # (within 5 pixels), we set the limit to the mean check
311 # to the edge of the detector
312 lowerRangeSmall = int(xCore)-5
313 upperRangeSmall = int(xCore)+5
314 if lowerRangeSmall < 0:
315 lowerRangeSmall = 0
316 if upperRangeSmall > xmax:
317 upperRangeSmall = xmax
318 ampImageBG = numpy.median(maskedImage[amp.getBBox()].image.array)
319 edgeMedian = numpy.median(sliceImage[:50, lowerRangeSmall:upperRangeSmall])
320 if edgeMedian > (ampImageBG + itlEdgeBleedThreshold):
321
322 # We need an estimate of the maximum width
323 # of the edge bleed for our masking model
324 # so we now estimate it by measuring the width of
325 # areas above 60 percent of the saturation level
326 # close to the edge,
327 # in a cutout up to 100 pixels from the edge.
328 subImage = sliceImage[:100, :]
329 maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
330 axis=1))
331
332 # Mask edge bleed with a
333 # decaying exponential model
334 for y in range(200):
335 edgeBleedHalfWidth = \
336 int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
337 + widthSat)/2.)
338 lowerRange = int(xCore)-edgeBleedHalfWidth
339 upperRange = int(xCore)+edgeBleedHalfWidth
340 # If the edge bleed model goes outside the detector
341 # we set the limit for the masking
342 # to the edge of the detector
343 if lowerRange < 0:
344 lowerRange = 0
345 if upperRange > xmax:
346 upperRange = xmax
347 sliceMask[y, lowerRange:upperRange] = satMaskBit
348
349
350def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
351 maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
352 """Interpolate over defects identified by a particular set of mask planes.
353
354 Parameters
355 ----------
356 maskedImage : `lsst.afw.image.MaskedImage`
357 Image to process.
358 fwhm : `float`
359 FWHM of double Gaussian smoothing kernel.
360 growSaturatedFootprints : scalar, optional
361 Number of pixels to grow footprints for saturated pixels.
362 maskNameList : `List` of `str`, optional
363 Mask plane name.
364 fallbackValue : scalar, optional
365 Value of last resort for interpolation.
366
367 Notes
368 -----
369 The ``fwhm`` parameter is used to create a PSF, but the underlying
370 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
371 not currently make use of this information.
372 """
373 mask = maskedImage.getMask()
374
375 if growSaturatedFootprints > 0 and "SAT" in maskNameList:
376 # If we are interpolating over an area larger than the original masked
377 # region, we need to expand the original mask bit to the full area to
378 # explain why we interpolated there.
379 growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
380
381 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
382 fpSet = afwDetection.FootprintSet(mask, thresh)
383 defectList = Defects.fromFootprintList(fpSet.getFootprints())
384
385 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
386 maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
387
388 return maskedImage
389
390
391def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
392 fallbackValue=None, useLegacyInterp=True):
393 """Mark saturated pixels and optionally interpolate over them
394
395 Parameters
396 ----------
397 maskedImage : `lsst.afw.image.MaskedImage`
398 Image to process.
399 saturation : scalar
400 Saturation level used as the detection threshold.
401 fwhm : `float`
402 FWHM of double Gaussian smoothing kernel.
403 growFootprints : scalar, optional
404 Number of pixels to grow footprints of detected regions.
405 interpolate : Bool, optional
406 If True, saturated pixels are interpolated over.
407 maskName : str, optional
408 Mask plane name.
409 fallbackValue : scalar, optional
410 Value of last resort for interpolation.
411
412 Notes
413 -----
414 The ``fwhm`` parameter is used to create a PSF, but the underlying
415 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
416 not currently make use of this information.
417 """
418 defectList = makeThresholdMask(
419 maskedImage=maskedImage,
420 threshold=saturation,
421 growFootprints=growFootprints,
422 maskName=maskName,
423 )
424 if interpolate:
425 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
426 maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
427
428 return maskedImage
429
430
431def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
432 """Compute number of edge trim pixels to match the calibration data.
433
434 Use the dimension difference between the raw exposure and the
435 calibration exposure to compute the edge trim pixels. This trim
436 is applied symmetrically, with the same number of pixels masked on
437 each side.
438
439 Parameters
440 ----------
441 rawMaskedImage : `lsst.afw.image.MaskedImage`
442 Image to trim.
443 calibMaskedImage : `lsst.afw.image.MaskedImage`
444 Calibration image to draw new bounding box from.
445
446 Returns
447 -------
448 replacementMaskedImage : `lsst.afw.image.MaskedImage`
449 ``rawMaskedImage`` trimmed to the appropriate size.
450
451 Raises
452 ------
453 RuntimeError
454 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
455 match ``calibMaskedImage``.
456 """
457 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
458 if nx != ny:
459 raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
460 if nx % 2 != 0:
461 raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
462 if nx < 0:
463 raise RuntimeError("Calibration maskedImage is larger than raw data.")
464
465 nEdge = nx//2
466 if nEdge > 0:
467 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
468 SourceDetectionTask.setEdgeBits(
469 rawMaskedImage,
470 replacementMaskedImage.getBBox(),
471 rawMaskedImage.getMask().getPlaneBitMask("EDGE")
472 )
473 else:
474 replacementMaskedImage = rawMaskedImage
475
476 return replacementMaskedImage
477
478
479def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
480 """Apply bias correction in place.
481
482 Parameters
483 ----------
484 maskedImage : `lsst.afw.image.MaskedImage`
485 Image to process. The image is modified by this method.
486 biasMaskedImage : `lsst.afw.image.MaskedImage`
487 Bias image of the same size as ``maskedImage``
488 trimToFit : `Bool`, optional
489 If True, raw data is symmetrically trimmed to match
490 calibration size.
491
492 Raises
493 ------
494 RuntimeError
495 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
496 the same size.
497
498 """
499 if trimToFit:
500 maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
501
502 if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
503 raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
504 (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
505 maskedImage -= biasMaskedImage
506
507
508def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
509 """Apply dark correction in place.
510
511 Parameters
512 ----------
513 maskedImage : `lsst.afw.image.MaskedImage`
514 Image to process. The image is modified by this method.
515 darkMaskedImage : `lsst.afw.image.MaskedImage`
516 Dark image of the same size as ``maskedImage``.
517 expScale : scalar
518 Dark exposure time for ``maskedImage``.
519 darkScale : scalar
520 Dark exposure time for ``darkMaskedImage``.
521 invert : `Bool`, optional
522 If True, re-add the dark to an already corrected image.
523 trimToFit : `Bool`, optional
524 If True, raw data is symmetrically trimmed to match
525 calibration size.
526
527 Raises
528 ------
529 RuntimeError
530 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
531 the same size.
532
533 Notes
534 -----
535 The dark correction is applied by calculating:
536 maskedImage -= dark * expScaling / darkScaling
537 """
538 if trimToFit:
539 maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
540
541 if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
542 raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
543 (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
544
545 scale = expScale / darkScale
546 if not invert:
547 maskedImage.scaledMinus(scale, darkMaskedImage)
548 else:
549 maskedImage.scaledPlus(scale, darkMaskedImage)
550
551
552def updateVariance(maskedImage, gain, readNoise, replace=True):
553 """Set the variance plane based on the image plane.
554
555 The maskedImage must have units of `adu` (if gain != 1.0) or
556 electron (if gain == 1.0). This routine will always produce a
557 variance plane in the same units as the image.
558
559 Parameters
560 ----------
561 maskedImage : `lsst.afw.image.MaskedImage`
562 Image to process. The variance plane is modified.
563 gain : scalar
564 The amplifier gain in electron/adu.
565 readNoise : scalar
566 The amplifier read noise in electron/pixel.
567 replace : `bool`, optional
568 Replace the current variance? If False, the image
569 variance will be added to the current variance plane.
570 """
571 var = maskedImage.variance
572 if replace:
573 var[:, :] = maskedImage.image
574 else:
575 var[:, :] += maskedImage.image
576 var /= gain
577 var += (readNoise/gain)**2
578
579
580def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
581 """Apply flat correction in place.
582
583 Parameters
584 ----------
585 maskedImage : `lsst.afw.image.MaskedImage`
586 Image to process. The image is modified.
587 flatMaskedImage : `lsst.afw.image.MaskedImage`
588 Flat image of the same size as ``maskedImage``
589 scalingType : str
590 Flat scale computation method. Allowed values are 'MEAN',
591 'MEDIAN', or 'USER'.
592 userScale : scalar, optional
593 Scale to use if ``scalingType='USER'``.
594 invert : `Bool`, optional
595 If True, unflatten an already flattened image.
596 trimToFit : `Bool`, optional
597 If True, raw data is symmetrically trimmed to match
598 calibration size.
599
600 Raises
601 ------
602 RuntimeError
603 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
604 the same size or if ``scalingType`` is not an allowed value.
605 """
606 if trimToFit:
607 maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
608
609 if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
610 raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
611 (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
612
613 # Figure out scale from the data
614 # Ideally the flats are normalized by the calibration product pipeline,
615 # but this allows some flexibility in the case that the flat is created by
616 # some other mechanism.
617 if scalingType in ('MEAN', 'MEDIAN'):
618 scalingType = afwMath.stringToStatisticsProperty(scalingType)
619 flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
620 elif scalingType == 'USER':
621 flatScale = userScale
622 else:
623 raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
624
625 if not invert:
626 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
627 else:
628 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
629
630
631def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
632 """Apply illumination correction in place.
633
634 Parameters
635 ----------
636 maskedImage : `lsst.afw.image.MaskedImage`
637 Image to process. The image is modified.
638 illumMaskedImage : `lsst.afw.image.MaskedImage`
639 Illumination correction image of the same size as ``maskedImage``.
640 illumScale : scalar
641 Scale factor for the illumination correction.
642 trimToFit : `Bool`, optional
643 If True, raw data is symmetrically trimmed to match
644 calibration size.
645
646 Raises
647 ------
648 RuntimeError
649 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
650 the same size.
651 """
652 if trimToFit:
653 maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
654
655 if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
656 raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
657 (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
658
659 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
660
661
662def brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None):
663 """Apply brighter fatter correction in place for the image.
664
665 Parameters
666 ----------
667 exposure : `lsst.afw.image.Exposure`
668 Exposure to have brighter-fatter correction applied. Modified
669 by this method.
670 kernel : `numpy.ndarray`
671 Brighter-fatter kernel to apply.
672 maxIter : scalar
673 Number of correction iterations to run.
674 threshold : scalar
675 Convergence threshold in terms of the sum of absolute
676 deviations between an iteration and the previous one.
677 applyGain : `Bool`
678 If True, then the exposure values are scaled by the gain prior
679 to correction.
680 gains : `dict` [`str`, `float`]
681 A dictionary, keyed by amplifier name, of the gains to use.
682 If gains is None, the nominal gains in the amplifier object are used.
683
684 Returns
685 -------
686 diff : `float`
687 Final difference between iterations achieved in correction.
688 iteration : `int`
689 Number of iterations used to calculate correction.
690
691 Notes
692 -----
693 This correction takes a kernel that has been derived from flat
694 field images to redistribute the charge. The gradient of the
695 kernel is the deflection field due to the accumulated charge.
696
697 Given the original image I(x) and the kernel K(x) we can compute
698 the corrected image Ic(x) using the following equation:
699
700 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
701
702 To evaluate the derivative term we expand it as follows:
703
704 0.5 * ( d/dx(I(x))*d/dx(int(dy*K(x-y)*I(y)))
705 + I(x)*d^2/dx^2(int(dy* K(x-y)*I(y))) )
706
707 Because we use the measured counts instead of the incident counts
708 we apply the correction iteratively to reconstruct the original
709 counts and the correction. We stop iterating when the summed
710 difference between the current corrected image and the one from
711 the previous iteration is below the threshold. We do not require
712 convergence because the number of iterations is too large a
713 computational cost. How we define the threshold still needs to be
714 evaluated, the current default was shown to work reasonably well
715 on a small set of images. For more information on the method see
716 DocuShare Document-19407.
717
718 The edges as defined by the kernel are not corrected because they
719 have spurious values due to the convolution.
720 """
721 image = exposure.getMaskedImage().getImage()
722
723 # The image needs to be units of electrons/holes
724 with gainContext(exposure, image, applyGain, gains):
725
726 kLx = numpy.shape(kernel)[0]
727 kLy = numpy.shape(kernel)[1]
728 kernelImage = afwImage.ImageD(kLx, kLy)
729 kernelImage.getArray()[:, :] = kernel
730 tempImage = afwImage.ImageD(image, deep=True)
731
732 nanIndex = numpy.isnan(tempImage.getArray())
733 tempImage.getArray()[nanIndex] = 0.
734
735 outImage = afwImage.ImageD(image.getDimensions())
736 corr = numpy.zeros(image.array.shape, dtype=numpy.float64)
737 prev_image = numpy.zeros(image.array.shape, dtype=numpy.float64)
738 convCntrl = afwMath.ConvolutionControl(False, True, 1)
739 fixedKernel = afwMath.FixedKernel(kernelImage)
740
741 # Define boundary by convolution region. The region that the
742 # correction will be calculated for is one fewer in each dimension
743 # because of the second derivative terms.
744 # NOTE: these need to use integer math, as we're using start:end as
745 # numpy index ranges.
746 startX = kLx//2
747 endX = -kLx//2
748 startY = kLy//2
749 endY = -kLy//2
750
751 for iteration in range(maxIter):
752
753 afwMath.convolve(outImage, tempImage, fixedKernel, convCntrl)
754 tmpArray = tempImage.getArray()
755 outArray = outImage.getArray()
756
757 with numpy.errstate(invalid="ignore", over="ignore"):
758 # First derivative term
759 gradTmp = numpy.gradient(tmpArray[startY:endY, startX:endX])
760 gradOut = numpy.gradient(outArray[startY:endY, startX:endX])
761 first = (gradTmp[0]*gradOut[0] + gradTmp[1]*gradOut[1])[1:-1, 1:-1]
762
763 # Second derivative term
764 diffOut20 = numpy.diff(outArray, 2, 0)[startY:endY, startX + 1:endX - 1]
765 diffOut21 = numpy.diff(outArray, 2, 1)[startY + 1:endY - 1, startX:endX]
766 second = tmpArray[startY + 1:endY - 1, startX + 1:endX - 1]*(diffOut20 + diffOut21)
767
768 corr[startY + 1:endY - 1, startX + 1:endX - 1] = 0.5*(first + second)
769
770 tmpArray[:, :] = image.getArray()[:, :]
771 tmpArray[nanIndex] = 0.
772 tmpArray[startY:endY, startX:endX] += corr[startY:endY, startX:endX]
773
774 if iteration > 0:
775 diff = numpy.sum(numpy.abs(prev_image - tmpArray), dtype=numpy.float64)
776
777 if diff < threshold:
778 break
779 prev_image[:, :] = tmpArray[:, :]
780
781 image.getArray()[startY + 1:endY - 1, startX + 1:endX - 1] += \
782 corr[startY + 1:endY - 1, startX + 1:endX - 1]
783
784 return diff, iteration
785
786
787def transferFlux(cFunc, fStep, correctionMode=True):
788 """Take the input convolved deflection potential and the flux array
789 to compute and apply the flux transfer into the correction array.
790
791 Parameters
792 ----------
793 cFunc: `numpy.array`
794 Deflection potential, being the convolution of the flux F with the
795 kernel K.
796 fStep: `numpy.array`
797 The array of flux values which act as the source of the flux transfer.
798 correctionMode: `bool`
799 Defines if applying correction (True) or generating sims (False).
800
801 Returns
802 -------
803 corr:
804 BFE correction array
805 """
806
807 if cFunc.shape != fStep.shape:
808 raise RuntimeError(f'transferFlux: array shapes do not match: {cFunc.shape}, {fStep.shape}')
809
810 # set the sign of the correction and set its value for the
811 # time averaged solution
812 if correctionMode:
813 # negative sign if applying BFE correction
814 factor = -0.5
815 else:
816 # positive sign if generating BFE simulations
817 factor = 0.5
818
819 # initialise the BFE correction image to zero
820 corr = numpy.zeros(cFunc.shape, dtype=numpy.float64)
821
822 # Generate a 2D mesh of x,y coordinates
823 yDim, xDim = cFunc.shape
824 y = numpy.arange(yDim, dtype=int)
825 x = numpy.arange(xDim, dtype=int)
826 xc, yc = numpy.meshgrid(x, y)
827
828 # process each axis in turn
829 for ax in [0, 1]:
830
831 # gradient of phi on right/upper edge of pixel
832 diff = numpy.diff(cFunc, axis=ax)
833
834 # expand array back to full size with zero gradient at the end
835 gx = numpy.zeros(cFunc.shape, dtype=numpy.float64)
836 yDiff, xDiff = diff.shape
837 gx[:yDiff, :xDiff] += diff
838
839 # select pixels with either positive gradients on the right edge,
840 # flux flowing to the right/up
841 # or negative gradients, flux flowing to the left/down
842 for i, sel in enumerate([gx > 0, gx < 0]):
843 xSelPixels = xc[sel]
844 ySelPixels = yc[sel]
845 # and add the flux into the pixel to the right or top
846 # depending on which axis we are handling
847 if ax == 0:
848 xPix = xSelPixels
849 yPix = ySelPixels+1
850 else:
851 xPix = xSelPixels+1
852 yPix = ySelPixels
853 # define flux as the either current pixel value or pixel
854 # above/right
855 # depending on whether positive or negative gradient
856 if i == 0:
857 # positive gradients, flux flowing to higher coordinate values
858 flux = factor * fStep[sel]*gx[sel]
859 else:
860 # negative gradients, flux flowing to lower coordinate values
861 flux = factor * fStep[yPix, xPix]*gx[sel]
862 # change the fluxes of the donor and receiving pixels
863 # such that flux is conserved
864 corr[sel] -= flux
865 corr[yPix, xPix] += flux
866
867 # return correction array
868 return corr
869
870
871def fluxConservingBrighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain,
872 gains=None, correctionMode=True):
873 """Apply brighter fatter correction in place for the image.
874
875 This version presents a modified version of the algorithm
876 found in ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``
877 which conserves the image flux, resulting in improved
878 correction of the cores of stars. The convolution has also been
879 modified to mitigate edge effects.
880
881 Parameters
882 ----------
883 exposure : `lsst.afw.image.Exposure`
884 Exposure to have brighter-fatter correction applied. Modified
885 by this method.
886 kernel : `numpy.ndarray`
887 Brighter-fatter kernel to apply.
888 maxIter : scalar
889 Number of correction iterations to run.
890 threshold : scalar
891 Convergence threshold in terms of the sum of absolute
892 deviations between an iteration and the previous one.
893 applyGain : `Bool`
894 If True, then the exposure values are scaled by the gain prior
895 to correction.
896 gains : `dict` [`str`, `float`]
897 A dictionary, keyed by amplifier name, of the gains to use.
898 If gains is None, the nominal gains in the amplifier object are used.
899 correctionMode : `Bool`
900 If True (default) the function applies correction for BFE. If False,
901 the code can instead be used to generate a simulation of BFE (sign
902 change in the direction of the effect)
903
904 Returns
905 -------
906 diff : `float`
907 Final difference between iterations achieved in correction.
908 iteration : `int`
909 Number of iterations used to calculate correction.
910
911 Notes
912 -----
913 Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.
914
915 This correction takes a kernel that has been derived from flat
916 field images to redistribute the charge. The gradient of the
917 kernel is the deflection field due to the accumulated charge.
918
919 Given the original image I(x) and the kernel K(x) we can compute
920 the corrected image Ic(x) using the following equation:
921
922 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
923
924 Improved algorithm at this step applies the divergence theorem to
925 obtain a pixelised correction.
926
927 Because we use the measured counts instead of the incident counts
928 we apply the correction iteratively to reconstruct the original
929 counts and the correction. We stop iterating when the summed
930 difference between the current corrected image and the one from
931 the previous iteration is below the threshold. We do not require
932 convergence because the number of iterations is too large a
933 computational cost. How we define the threshold still needs to be
934 evaluated, the current default was shown to work reasonably well
935 on a small set of images.
936
937 Edges are handled in the convolution by padding. This is still not
938 a physical model for the edge, but avoids discontinuity in the correction.
939
940 Author of modified version: Lance.Miller@physics.ox.ac.uk
941 (see DM-38555).
942 """
943 image = exposure.getMaskedImage().getImage()
944
945 # The image needs to be units of electrons/holes
946 with gainContext(exposure, image, applyGain, gains):
947
948 # get kernel and its shape
949 kLy, kLx = kernel.shape
950 kernelImage = afwImage.ImageD(kLx, kLy)
951 kernelImage.getArray()[:, :] = kernel
952 tempImage = afwImage.ImageD(image, deep=True)
953
954 nanIndex = numpy.isnan(tempImage.getArray())
955 tempImage.getArray()[nanIndex] = 0.
956
957 outImage = afwImage.ImageD(image.getDimensions())
958 corr = numpy.zeros(image.array.shape, dtype=numpy.float64)
959 prevImage = numpy.zeros(image.array.shape, dtype=numpy.float64)
960 convCntrl = afwMath.ConvolutionControl(False, False, 1)
961 fixedKernel = afwMath.FixedKernel(kernelImage)
962
963 # set the padding amount
964 # ensure we pad by an even amount larger than the kernel
965 kLy = 2 * ((1+kLy)//2)
966 kLx = 2 * ((1+kLx)//2)
967
968 # The deflection potential only depends on the gradient of
969 # the convolution, so we can subtract the mean, which then
970 # allows us to pad the image with zeros and avoid wrap-around effects
971 # (although still not handling the image edges with a physical model)
972 # This wouldn't be great if there were a strong image gradient.
973 imYdimension, imXdimension = tempImage.array.shape
974 imean = numpy.mean(tempImage.getArray()[~nanIndex], dtype=numpy.float64)
975 # subtract mean from image
976 tempImage -= imean
977 tempImage.array[nanIndex] = 0.0
978 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
979 outImage = afwImage.ImageD(numpy.pad(outImage.getArray(), ((0, kLy), (0, kLx))))
980 # Convert array to afw image so afwMath.convolve works
981 padImage = afwImage.ImageD(padArray.shape[1], padArray.shape[0])
982 padImage.array[:] = padArray
983
984 for iteration in range(maxIter):
985
986 # create deflection potential, convolution of flux with kernel
987 # using padded counts array
988 afwMath.convolve(outImage, padImage, fixedKernel, convCntrl)
989 tmpArray = tempImage.getArray()
990 outArray = outImage.getArray()
991
992 # trim convolution output back to original shape
993 outArray = outArray[:imYdimension, :imXdimension]
994
995 # generate the correction array, with correctionMode set as input
996 corr[...] = transferFlux(outArray, tmpArray, correctionMode=correctionMode)
997
998 # update the arrays for the next iteration
999 tmpArray[:, :] = image.getArray()[:, :]
1000 tmpArray += corr
1001 tmpArray[nanIndex] = 0.
1002 # update padded array
1003 # subtract mean
1004 tmpArray -= imean
1005 tempImage.array[nanIndex] = 0.
1006 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
1007
1008 if iteration > 0:
1009 diff = numpy.sum(numpy.abs(prevImage - tmpArray), dtype=numpy.float64)
1010
1011 if diff < threshold:
1012 break
1013 prevImage[:, :] = tmpArray[:, :]
1014
1015 image.getArray()[:] += corr[:]
1016
1017 return diff, iteration
1018
1019
1020@contextmanager
1021def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
1022 """Context manager that applies and removes gain.
1023
1024 Parameters
1025 ----------
1026 exp : `lsst.afw.image.Exposure`
1027 Exposure to apply/remove gain.
1028 image : `lsst.afw.image.Image`
1029 Image to apply/remove gain.
1030 apply : `bool`
1031 If True, apply and remove the amplifier gain.
1032 gains : `dict` [`str`, `float`], optional
1033 A dictionary, keyed by amplifier name, of the gains to use.
1034 If gains is None, the nominal gains in the amplifier object are used.
1035 invert : `bool`, optional
1036 Invert the gains (e.g. convert electrons to adu temporarily)?
1037 isTrimmed : `bool`, optional
1038 Is this a trimmed exposure?
1039
1040 Yields
1041 ------
1042 exp : `lsst.afw.image.Exposure`
1043 Exposure with the gain applied.
1044 """
1045 # check we have all of them if provided because mixing and matching would
1046 # be a real mess
1047 if gains and apply is True:
1048 ampNames = [amp.getName() for amp in exp.getDetector()]
1049 for ampName in ampNames:
1050 if ampName not in gains.keys():
1051 raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
1052
1053 if apply:
1054 ccd = exp.getDetector()
1055 for amp in ccd:
1056 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1057 if gains:
1058 gain = gains[amp.getName()]
1059 else:
1060 gain = amp.getGain()
1061 if invert:
1062 sim /= gain
1063 else:
1064 sim *= gain
1065
1066 try:
1067 yield exp
1068 finally:
1069 if apply:
1070 ccd = exp.getDetector()
1071 for amp in ccd:
1072 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1073 if gains:
1074 gain = gains[amp.getName()]
1075 else:
1076 gain = amp.getGain()
1077 if invert:
1078 sim *= gain
1079 else:
1080 sim /= gain
1081
1082
1083def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
1084 sensorTransmission=None, atmosphereTransmission=None):
1085 """Attach a TransmissionCurve to an Exposure, given separate curves for
1086 different components.
1087
1088 Parameters
1089 ----------
1090 exposure : `lsst.afw.image.Exposure`
1091 Exposure object to modify by attaching the product of all given
1092 ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
1093 Must have a valid ``Detector`` attached that matches the detector
1094 associated with sensorTransmission.
1095 opticsTransmission : `lsst.afw.image.TransmissionCurve`
1096 A ``TransmissionCurve`` that represents the throughput of the optics,
1097 to be evaluated in focal-plane coordinates.
1098 filterTransmission : `lsst.afw.image.TransmissionCurve`
1099 A ``TransmissionCurve`` that represents the throughput of the filter
1100 itself, to be evaluated in focal-plane coordinates.
1101 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1102 A ``TransmissionCurve`` that represents the throughput of the sensor
1103 itself, to be evaluated in post-assembly trimmed detector coordinates.
1104 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1105 A ``TransmissionCurve`` that represents the throughput of the
1106 atmosphere, assumed to be spatially constant.
1107
1108 Returns
1109 -------
1110 combined : `lsst.afw.image.TransmissionCurve`
1111 The TransmissionCurve attached to the exposure.
1112
1113 Notes
1114 -----
1115 All ``TransmissionCurve`` arguments are optional; if none are provided, the
1116 attached ``TransmissionCurve`` will have unit transmission everywhere.
1117 """
1118 combined = afwImage.TransmissionCurve.makeIdentity()
1119 if atmosphereTransmission is not None:
1120 combined *= atmosphereTransmission
1121 if opticsTransmission is not None:
1122 combined *= opticsTransmission
1123 if filterTransmission is not None:
1124 combined *= filterTransmission
1125 detector = exposure.getDetector()
1126 fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
1127 toSys=camGeom.PIXELS)
1128 combined = combined.transformedBy(fpToPix)
1129 if sensorTransmission is not None:
1130 combined *= sensorTransmission
1131 exposure.getInfo().setTransmissionCurve(combined)
1132 return combined
1133
1134
1135def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
1136 """Scale an exposure by the amplifier gains.
1137
1138 Parameters
1139 ----------
1140 exposure : `lsst.afw.image.Exposure`
1141 Exposure to process. The image is modified.
1142 normalizeGains : `Bool`, optional
1143 If True, then amplifiers are scaled to force the median of
1144 each amplifier to equal the median of those medians.
1145 ptcGains : `dict`[`str`], optional
1146 Dictionary keyed by amp name containing the PTC gains.
1147 isTrimmed : `bool`, optional
1148 Is the input image trimmed?
1149 """
1150 ccd = exposure.getDetector()
1151 ccdImage = exposure.getMaskedImage()
1152
1153 medians = []
1154 for amp in ccd:
1155 if isTrimmed:
1156 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1157 else:
1158 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1159 if ptcGains:
1160 sim *= ptcGains[amp.getName()]
1161 else:
1162 sim *= amp.getGain()
1163
1164 if normalizeGains:
1165 medians.append(numpy.median(sim.getImage().getArray()))
1166
1167 if normalizeGains:
1168 median = numpy.median(numpy.array(medians))
1169 for index, amp in enumerate(ccd):
1170 if isTrimmed:
1171 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1172 else:
1173 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1174 if medians[index] != 0.0:
1175 sim *= median/medians[index]
1176
1177
1179 """Grow the saturation trails by an amount dependent on the width of the
1180 trail.
1181
1182 Parameters
1183 ----------
1184 mask : `lsst.afw.image.Mask`
1185 Mask which will have the saturated areas grown.
1186 """
1187
1188 extraGrowDict = {}
1189 for i in range(1, 6):
1190 extraGrowDict[i] = 0
1191 for i in range(6, 8):
1192 extraGrowDict[i] = 1
1193 for i in range(8, 10):
1194 extraGrowDict[i] = 3
1195 extraGrowMax = 4
1196
1197 if extraGrowMax <= 0:
1198 return
1199
1200 saturatedBit = mask.getPlaneBitMask("SAT")
1201
1202 xmin, ymin = mask.getBBox().getMin()
1203 width = mask.getWidth()
1204
1205 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1206 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1207
1208 for fp in fpList:
1209 for s in fp.getSpans():
1210 x0, x1 = s.getX0(), s.getX1()
1211
1212 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1213 if extraGrow > 0:
1214 y = s.getY() - ymin
1215 x0 -= xmin + extraGrow
1216 x1 -= xmin - extraGrow
1217
1218 if x0 < 0:
1219 x0 = 0
1220 if x1 >= width - 1:
1221 x1 = width - 1
1222
1223 mask.array[y, x0:x1+1] |= saturatedBit
1224
1225
1226def setBadRegions(exposure, badStatistic="MEDIAN"):
1227 """Set all BAD areas of the chip to the average of the rest of the exposure
1228
1229 Parameters
1230 ----------
1231 exposure : `lsst.afw.image.Exposure`
1232 Exposure to mask. The exposure mask is modified.
1233 badStatistic : `str`, optional
1234 Statistic to use to generate the replacement value from the
1235 image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
1236
1237 Returns
1238 -------
1239 badPixelCount : scalar
1240 Number of bad pixels masked.
1241 badPixelValue : scalar
1242 Value substituted for bad pixels.
1243
1244 Raises
1245 ------
1246 RuntimeError
1247 Raised if `badStatistic` is not an allowed value.
1248 """
1249 if badStatistic == "MEDIAN":
1250 statistic = afwMath.MEDIAN
1251 elif badStatistic == "MEANCLIP":
1252 statistic = afwMath.MEANCLIP
1253 else:
1254 raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
1255
1256 mi = exposure.getMaskedImage()
1257 mask = mi.getMask()
1258 BAD = mask.getPlaneBitMask("BAD")
1259 INTRP = mask.getPlaneBitMask("INTRP")
1260
1262 sctrl.setAndMask(BAD)
1263 value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
1264
1265 maskArray = mask.getArray()
1266 imageArray = mi.getImage().getArray()
1267 badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
1268 imageArray[:] = numpy.where(badPixels, value, imageArray)
1269
1270 return badPixels.sum(), value
1271
1272
1273def checkFilter(exposure, filterList, log):
1274 """Check to see if an exposure is in a filter specified by a list.
1275
1276 The goal of this is to provide a unified filter checking interface
1277 for all filter dependent stages.
1278
1279 Parameters
1280 ----------
1281 exposure : `lsst.afw.image.Exposure`
1282 Exposure to examine.
1283 filterList : `list` [`str`]
1284 List of physical_filter names to check.
1285 log : `logging.Logger`
1286 Logger to handle messages.
1287
1288 Returns
1289 -------
1290 result : `bool`
1291 True if the exposure's filter is contained in the list.
1292 """
1293 if len(filterList) == 0:
1294 return False
1295 thisFilter = exposure.getFilter()
1296 if thisFilter is None:
1297 log.warning("No FilterLabel attached to this exposure!")
1298 return False
1299
1300 thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
1301 if thisPhysicalFilter in filterList:
1302 return True
1303 elif thisFilter.bandLabel in filterList:
1304 if log:
1305 log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
1306 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1307 return True
1308 else:
1309 return False
1310
1311
1312def getPhysicalFilter(filterLabel, log):
1313 """Get the physical filter label associated with the given filterLabel.
1314
1315 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1316 associated with the given ``filterLabel``, the returned label will be
1317 "Unknown".
1318
1319 Parameters
1320 ----------
1321 filterLabel : `lsst.afw.image.FilterLabel`
1322 The `lsst.afw.image.FilterLabel` object from which to derive the
1323 physical filter label.
1324 log : `logging.Logger`
1325 Logger to handle messages.
1326
1327 Returns
1328 -------
1329 physicalFilter : `str`
1330 The value returned by the physicalLabel attribute of ``filterLabel`` if
1331 it exists, otherwise set to \"Unknown\".
1332 """
1333 if filterLabel is None:
1334 physicalFilter = "Unknown"
1335 log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
1336 else:
1337 try:
1338 physicalFilter = filterLabel.physicalLabel
1339 except RuntimeError:
1340 log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
1341 physicalFilter = "Unknown"
1342 return physicalFilter
1343
1344
1345def countMaskedPixels(maskedIm, maskPlane):
1346 """Count the number of pixels in a given mask plane.
1347
1348 Parameters
1349 ----------
1350 maskedIm : `~lsst.afw.image.MaskedImage`
1351 Masked image to examine.
1352 maskPlane : `str`
1353 Name of the mask plane to examine.
1354
1355 Returns
1356 -------
1357 nPix : `int`
1358 Number of pixels in the requested mask plane.
1359 """
1360 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1361 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1362 return nPix
1363
1364
1365def getExposureGains(exposure):
1366 """Get the per-amplifier gains used for this exposure.
1367
1368 Parameters
1369 ----------
1370 exposure : `lsst.afw.image.Exposure`
1371 The exposure to find gains for.
1372
1373 Returns
1374 -------
1375 gains : `dict` [`str` `float`]
1376 Dictionary of gain values, keyed by amplifier name.
1377 Returns empty dict when detector is None.
1378 """
1379 det = exposure.getDetector()
1380 if det is None:
1381 return dict()
1382
1383 metadata = exposure.getMetadata()
1384 gains = {}
1385 for amp in det:
1386 ampName = amp.getName()
1387 # The key may use the new LSST ISR or the old LSST prefix
1388 if (key1 := f"LSST ISR GAIN {ampName}") in metadata:
1389 gains[ampName] = metadata[key1]
1390 elif (key2 := f"LSST GAIN {ampName}") in metadata:
1391 gains[ampName] = metadata[key2]
1392 else:
1393 gains[ampName] = amp.getGain()
1394 return gains
1395
1396
1398 """Get the per-amplifier read noise used for this exposure.
1399
1400 Parameters
1401 ----------
1402 exposure : `lsst.afw.image.Exposure`
1403 The exposure to find read noise for.
1404
1405 Returns
1406 -------
1407 readnoises : `dict` [`str` `float`]
1408 Dictionary of read noise values, keyed by amplifier name.
1409 Returns empty dict when detector is None.
1410 """
1411 det = exposure.getDetector()
1412 if det is None:
1413 return dict()
1414
1415 metadata = exposure.getMetadata()
1416 readnoises = {}
1417 for amp in det:
1418 ampName = amp.getName()
1419 # The key may use the new LSST ISR or the old LSST prefix
1420 if (key1 := f"LSST ISR READNOISE {ampName}") in metadata:
1421 readnoises[ampName] = metadata[key1]
1422 elif (key2 := f"LSST READNOISE {ampName}") in metadata:
1423 readnoises[ampName] = metadata[key2]
1424 else:
1425 readnoises[ampName] = amp.getReadNoise()
1426 return readnoises
1427
1428
1429def isTrimmedExposure(exposure):
1430 """Check if the unused pixels (pre-/over-scan pixels) have
1431 been trimmed from an exposure.
1432
1433 Parameters
1434 ----------
1435 exposure : `lsst.afw.image.Exposure`
1436 The exposure to check.
1437
1438 Returns
1439 -------
1440 result : `bool`
1441 True if the image is trimmed, else False.
1442 """
1443 return exposure.getDetector().getBBox() == exposure.getBBox()
1444
1445
1446def isTrimmedImage(image, detector):
1447 """Check if the unused pixels (pre-/over-scan pixels) have
1448 been trimmed from an image
1449
1450 Parameters
1451 ----------
1452 image : `lsst.afw.image.Image`
1453 The image to check.
1454 detector : `lsst.afw.cameraGeom.Detector`
1455 The detector associated with the image.
1456
1457 Returns
1458 -------
1459 result : `bool`
1460 True if the image is trimmed, else False.
1461 """
1462 return detector.getBBox() == image.getBBox()
Parameters to control convolution.
A kernel created from an Image.
Definition Kernel.h:471
Pass parameters to a Statistics object.
Definition Statistics.h:83
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True)
interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None, maskNameList=None, useLegacyInterp=True)
setBadRegions(exposure, badStatistic="MEDIAN")
countMaskedPixels(maskedIm, maskPlane)
attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None, sensorTransmission=None, atmosphereTransmission=None)
flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False)
illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True)
growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD")
transferFlux(cFunc, fStep, correctionMode=True)
biasCorrection(maskedImage, biasMaskedImage, trimToFit=False)
checkFilter(exposure, filterList, log)
darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False)
fluxConservingBrighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None, correctionMode=True)
brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None)
makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT')
transposeMaskedImage(maskedImage)
interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1, maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True)
isTrimmedImage(image, detector)
updateVariance(maskedImage, gain, readNoise, replace=True)
applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True)
getPhysicalFilter(filterLabel, log)
saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT', fallbackValue=None, useLegacyInterp=True)
maskITLEdgeBleed(ccdExposure, badAmpDict, itlEdgeBleedSatMinArea=10000, itlEdgeBleedSatMaxArea=100000, itlEdgeBleedThreshold=5000., itlEdgeBleedModelConstant=0.03, saturatedMaskName="SAT")
trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage)