LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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 "illuminationCorrection",
38 "interpolateDefectList",
39 "interpolateFromMask",
40 "makeThresholdMask",
41 "saturationCorrection",
42 "setBadRegions",
43 "transferFlux",
44 "transposeMaskedImage",
45 "trimToMatchCalibBBox",
46 "updateVariance",
47 "widenSaturationTrails",
48]
49
50import math
51import numpy
52
53import lsst.geom
54import lsst.afw.image as afwImage
55import lsst.afw.detection as afwDetection
56import lsst.afw.math as afwMath
57import lsst.meas.algorithms as measAlg
58import lsst.afw.cameraGeom as camGeom
59
60from lsst.meas.algorithms.detection import SourceDetectionTask
61
62from contextlib import contextmanager
63
64from .defects import Defects
65
66
67def createPsf(fwhm):
68 """Make a double Gaussian PSF.
69
70 Parameters
71 ----------
72 fwhm : scalar
73 FWHM of double Gaussian smoothing kernel.
74
75 Returns
76 -------
77 psf : `lsst.meas.algorithms.DoubleGaussianPsf`
78 The created smoothing kernel.
79 """
80 ksize = 4*int(fwhm) + 1
81 return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
82
83
84def transposeMaskedImage(maskedImage):
85 """Make a transposed copy of a masked image.
86
87 Parameters
88 ----------
89 maskedImage : `lsst.afw.image.MaskedImage`
90 Image to process.
91
92 Returns
93 -------
94 transposed : `lsst.afw.image.MaskedImage`
95 The transposed copy of the input image.
96 """
97 transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
98 transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
99 transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
100 transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
101 return transposed
102
103
104def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None):
105 """Interpolate over defects specified in a defect list.
106
107 Parameters
108 ----------
109 maskedImage : `lsst.afw.image.MaskedImage`
110 Image to process.
111 defectList : `lsst.meas.algorithms.Defects`
112 List of defects to interpolate over.
113 fwhm : `float`
114 FWHM of double Gaussian smoothing kernel.
115 fallbackValue : scalar, optional
116 Fallback value if an interpolated value cannot be determined.
117 If None, then the clipped mean of the image is used.
118
119 Notes
120 -----
121 The ``fwhm`` parameter is used to create a PSF, but the underlying
122 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
123 not currently make use of this information.
124 """
125 psf = createPsf(fwhm)
126 if fallbackValue is None:
127 fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
128 if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
129 maskedImage.getMask().addMaskPlane('INTRP')
130 measAlg.interpolateOverDefects(maskedImage, psf, defectList, fallbackValue, True)
131 return maskedImage
132
133
134def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
135 """Mask pixels based on threshold detection.
136
137 Parameters
138 ----------
139 maskedImage : `lsst.afw.image.MaskedImage`
140 Image to process. Only the mask plane is updated.
141 threshold : scalar
142 Detection threshold.
143 growFootprints : scalar, optional
144 Number of pixels to grow footprints of detected regions.
145 maskName : str, optional
146 Mask plane name, or list of names to convert
147
148 Returns
149 -------
150 defectList : `lsst.meas.algorithms.Defects`
151 Defect list constructed from pixels above the threshold.
152 """
153 # find saturated regions
154 thresh = afwDetection.Threshold(threshold)
155 fs = afwDetection.FootprintSet(maskedImage, thresh)
156
157 if growFootprints > 0:
158 fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
159 fpList = fs.getFootprints()
160
161 # set mask
162 mask = maskedImage.getMask()
163 bitmask = mask.getPlaneBitMask(maskName)
164 afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
165
166 return Defects.fromFootprintList(fpList)
167
168
169def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
170 """Grow a mask by an amount and add to the requested plane.
171
172 Parameters
173 ----------
174 mask : `lsst.afw.image.Mask`
175 Mask image to process.
176 radius : scalar
177 Amount to grow the mask.
178 maskNameList : `str` or `list` [`str`]
179 Mask names that should be grown.
180 maskValue : `str`
181 Mask plane to assign the newly masked pixels to.
182 """
183 if radius > 0:
184 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
185 fpSet = afwDetection.FootprintSet(mask, thresh)
186 fpSet = afwDetection.FootprintSet(fpSet, rGrow=radius, isotropic=False)
187 fpSet.setMask(mask, maskValue)
188
189
190def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
191 maskNameList=['SAT'], fallbackValue=None):
192 """Interpolate over defects identified by a particular set of mask planes.
193
194 Parameters
195 ----------
196 maskedImage : `lsst.afw.image.MaskedImage`
197 Image to process.
198 fwhm : `float`
199 FWHM of double Gaussian smoothing kernel.
200 growSaturatedFootprints : scalar, optional
201 Number of pixels to grow footprints for saturated pixels.
202 maskNameList : `List` of `str`, optional
203 Mask plane name.
204 fallbackValue : scalar, optional
205 Value of last resort for interpolation.
206
207 Notes
208 -----
209 The ``fwhm`` parameter is used to create a PSF, but the underlying
210 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
211 not currently make use of this information.
212 """
213 mask = maskedImage.getMask()
214
215 if growSaturatedFootprints > 0 and "SAT" in maskNameList:
216 # If we are interpolating over an area larger than the original masked
217 # region, we need to expand the original mask bit to the full area to
218 # explain why we interpolated there.
219 growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
220
221 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
222 fpSet = afwDetection.FootprintSet(mask, thresh)
223 defectList = Defects.fromFootprintList(fpSet.getFootprints())
224
225 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue)
226
227 return maskedImage
228
229
230def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
231 fallbackValue=None):
232 """Mark saturated pixels and optionally interpolate over them
233
234 Parameters
235 ----------
236 maskedImage : `lsst.afw.image.MaskedImage`
237 Image to process.
238 saturation : scalar
239 Saturation level used as the detection threshold.
240 fwhm : `float`
241 FWHM of double Gaussian smoothing kernel.
242 growFootprints : scalar, optional
243 Number of pixels to grow footprints of detected regions.
244 interpolate : Bool, optional
245 If True, saturated pixels are interpolated over.
246 maskName : str, optional
247 Mask plane name.
248 fallbackValue : scalar, optional
249 Value of last resort for interpolation.
250
251 Notes
252 -----
253 The ``fwhm`` parameter is used to create a PSF, but the underlying
254 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
255 not currently make use of this information.
256 """
257 defectList = makeThresholdMask(
258 maskedImage=maskedImage,
259 threshold=saturation,
260 growFootprints=growFootprints,
261 maskName=maskName,
262 )
263 if interpolate:
264 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue)
265
266 return maskedImage
267
268
269def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
270 """Compute number of edge trim pixels to match the calibration data.
271
272 Use the dimension difference between the raw exposure and the
273 calibration exposure to compute the edge trim pixels. This trim
274 is applied symmetrically, with the same number of pixels masked on
275 each side.
276
277 Parameters
278 ----------
279 rawMaskedImage : `lsst.afw.image.MaskedImage`
280 Image to trim.
281 calibMaskedImage : `lsst.afw.image.MaskedImage`
282 Calibration image to draw new bounding box from.
283
284 Returns
285 -------
286 replacementMaskedImage : `lsst.afw.image.MaskedImage`
287 ``rawMaskedImage`` trimmed to the appropriate size.
288
289 Raises
290 ------
291 RuntimeError
292 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
293 match ``calibMaskedImage``.
294 """
295 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
296 if nx != ny:
297 raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
298 if nx % 2 != 0:
299 raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
300 if nx < 0:
301 raise RuntimeError("Calibration maskedImage is larger than raw data.")
302
303 nEdge = nx//2
304 if nEdge > 0:
305 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
306 SourceDetectionTask.setEdgeBits(
307 rawMaskedImage,
308 replacementMaskedImage.getBBox(),
309 rawMaskedImage.getMask().getPlaneBitMask("EDGE")
310 )
311 else:
312 replacementMaskedImage = rawMaskedImage
313
314 return replacementMaskedImage
315
316
317def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
318 """Apply bias correction in place.
319
320 Parameters
321 ----------
322 maskedImage : `lsst.afw.image.MaskedImage`
323 Image to process. The image is modified by this method.
324 biasMaskedImage : `lsst.afw.image.MaskedImage`
325 Bias image of the same size as ``maskedImage``
326 trimToFit : `Bool`, optional
327 If True, raw data is symmetrically trimmed to match
328 calibration size.
329
330 Raises
331 ------
332 RuntimeError
333 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
334 the same size.
335
336 """
337 if trimToFit:
338 maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
339
340 if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
341 raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
342 (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
343 maskedImage -= biasMaskedImage
344
345
346def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
347 """Apply dark correction in place.
348
349 Parameters
350 ----------
351 maskedImage : `lsst.afw.image.MaskedImage`
352 Image to process. The image is modified by this method.
353 darkMaskedImage : `lsst.afw.image.MaskedImage`
354 Dark image of the same size as ``maskedImage``.
355 expScale : scalar
356 Dark exposure time for ``maskedImage``.
357 darkScale : scalar
358 Dark exposure time for ``darkMaskedImage``.
359 invert : `Bool`, optional
360 If True, re-add the dark to an already corrected image.
361 trimToFit : `Bool`, optional
362 If True, raw data is symmetrically trimmed to match
363 calibration size.
364
365 Raises
366 ------
367 RuntimeError
368 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
369 the same size.
370
371 Notes
372 -----
373 The dark correction is applied by calculating:
374 maskedImage -= dark * expScaling / darkScaling
375 """
376 if trimToFit:
377 maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
378
379 if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
380 raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
381 (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
382
383 scale = expScale / darkScale
384 if not invert:
385 maskedImage.scaledMinus(scale, darkMaskedImage)
386 else:
387 maskedImage.scaledPlus(scale, darkMaskedImage)
388
389
390def updateVariance(maskedImage, gain, readNoise):
391 """Set the variance plane based on the image plane.
392
393 Parameters
394 ----------
395 maskedImage : `lsst.afw.image.MaskedImage`
396 Image to process. The variance plane is modified.
397 gain : scalar
398 The amplifier gain in electrons/ADU.
399 readNoise : scalar
400 The amplifier read nmoise in ADU/pixel.
401 """
402 var = maskedImage.getVariance()
403 var[:] = maskedImage.getImage()
404 var /= gain
405 var += readNoise**2
406
407
408def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
409 """Apply flat correction in place.
410
411 Parameters
412 ----------
413 maskedImage : `lsst.afw.image.MaskedImage`
414 Image to process. The image is modified.
415 flatMaskedImage : `lsst.afw.image.MaskedImage`
416 Flat image of the same size as ``maskedImage``
417 scalingType : str
418 Flat scale computation method. Allowed values are 'MEAN',
419 'MEDIAN', or 'USER'.
420 userScale : scalar, optional
421 Scale to use if ``scalingType='USER'``.
422 invert : `Bool`, optional
423 If True, unflatten an already flattened image.
424 trimToFit : `Bool`, optional
425 If True, raw data is symmetrically trimmed to match
426 calibration size.
427
428 Raises
429 ------
430 RuntimeError
431 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
432 the same size or if ``scalingType`` is not an allowed value.
433 """
434 if trimToFit:
435 maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
436
437 if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
438 raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
439 (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
440
441 # Figure out scale from the data
442 # Ideally the flats are normalized by the calibration product pipeline,
443 # but this allows some flexibility in the case that the flat is created by
444 # some other mechanism.
445 if scalingType in ('MEAN', 'MEDIAN'):
446 scalingType = afwMath.stringToStatisticsProperty(scalingType)
447 flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
448 elif scalingType == 'USER':
449 flatScale = userScale
450 else:
451 raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
452
453 if not invert:
454 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
455 else:
456 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
457
458
459def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
460 """Apply illumination correction in place.
461
462 Parameters
463 ----------
464 maskedImage : `lsst.afw.image.MaskedImage`
465 Image to process. The image is modified.
466 illumMaskedImage : `lsst.afw.image.MaskedImage`
467 Illumination correction image of the same size as ``maskedImage``.
468 illumScale : scalar
469 Scale factor for the illumination correction.
470 trimToFit : `Bool`, optional
471 If True, raw data is symmetrically trimmed to match
472 calibration size.
473
474 Raises
475 ------
476 RuntimeError
477 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
478 the same size.
479 """
480 if trimToFit:
481 maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
482
483 if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
484 raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
485 (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
486
487 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
488
489
490def brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None):
491 """Apply brighter fatter correction in place for the image.
492
493 Parameters
494 ----------
495 exposure : `lsst.afw.image.Exposure`
496 Exposure to have brighter-fatter correction applied. Modified
497 by this method.
498 kernel : `numpy.ndarray`
499 Brighter-fatter kernel to apply.
500 maxIter : scalar
501 Number of correction iterations to run.
502 threshold : scalar
503 Convergence threshold in terms of the sum of absolute
504 deviations between an iteration and the previous one.
505 applyGain : `Bool`
506 If True, then the exposure values are scaled by the gain prior
507 to correction.
508 gains : `dict` [`str`, `float`]
509 A dictionary, keyed by amplifier name, of the gains to use.
510 If gains is None, the nominal gains in the amplifier object are used.
511
512 Returns
513 -------
514 diff : `float`
515 Final difference between iterations achieved in correction.
516 iteration : `int`
517 Number of iterations used to calculate correction.
518
519 Notes
520 -----
521 This correction takes a kernel that has been derived from flat
522 field images to redistribute the charge. The gradient of the
523 kernel is the deflection field due to the accumulated charge.
524
525 Given the original image I(x) and the kernel K(x) we can compute
526 the corrected image Ic(x) using the following equation:
527
528 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
529
530 To evaluate the derivative term we expand it as follows:
531
532 0.5 * ( d/dx(I(x))*d/dx(int(dy*K(x-y)*I(y)))
533 + I(x)*d^2/dx^2(int(dy* K(x-y)*I(y))) )
534
535 Because we use the measured counts instead of the incident counts
536 we apply the correction iteratively to reconstruct the original
537 counts and the correction. We stop iterating when the summed
538 difference between the current corrected image and the one from
539 the previous iteration is below the threshold. We do not require
540 convergence because the number of iterations is too large a
541 computational cost. How we define the threshold still needs to be
542 evaluated, the current default was shown to work reasonably well
543 on a small set of images. For more information on the method see
544 DocuShare Document-19407.
545
546 The edges as defined by the kernel are not corrected because they
547 have spurious values due to the convolution.
548 """
549 image = exposure.getMaskedImage().getImage()
550
551 # The image needs to be units of electrons/holes
552 with gainContext(exposure, image, applyGain, gains):
553
554 kLx = numpy.shape(kernel)[0]
555 kLy = numpy.shape(kernel)[1]
556 kernelImage = afwImage.ImageD(kLx, kLy)
557 kernelImage.getArray()[:, :] = kernel
558 tempImage = image.clone()
559
560 nanIndex = numpy.isnan(tempImage.getArray())
561 tempImage.getArray()[nanIndex] = 0.
562
563 outImage = afwImage.ImageF(image.getDimensions())
564 corr = numpy.zeros_like(image.getArray())
565 prev_image = numpy.zeros_like(image.getArray())
566 convCntrl = afwMath.ConvolutionControl(False, True, 1)
567 fixedKernel = afwMath.FixedKernel(kernelImage)
568
569 # Define boundary by convolution region. The region that the
570 # correction will be calculated for is one fewer in each dimension
571 # because of the second derivative terms.
572 # NOTE: these need to use integer math, as we're using start:end as
573 # numpy index ranges.
574 startX = kLx//2
575 endX = -kLx//2
576 startY = kLy//2
577 endY = -kLy//2
578
579 for iteration in range(maxIter):
580
581 afwMath.convolve(outImage, tempImage, fixedKernel, convCntrl)
582 tmpArray = tempImage.getArray()
583 outArray = outImage.getArray()
584
585 with numpy.errstate(invalid="ignore", over="ignore"):
586 # First derivative term
587 gradTmp = numpy.gradient(tmpArray[startY:endY, startX:endX])
588 gradOut = numpy.gradient(outArray[startY:endY, startX:endX])
589 first = (gradTmp[0]*gradOut[0] + gradTmp[1]*gradOut[1])[1:-1, 1:-1]
590
591 # Second derivative term
592 diffOut20 = numpy.diff(outArray, 2, 0)[startY:endY, startX + 1:endX - 1]
593 diffOut21 = numpy.diff(outArray, 2, 1)[startY + 1:endY - 1, startX:endX]
594 second = tmpArray[startY + 1:endY - 1, startX + 1:endX - 1]*(diffOut20 + diffOut21)
595
596 corr[startY + 1:endY - 1, startX + 1:endX - 1] = 0.5*(first + second)
597
598 tmpArray[:, :] = image.getArray()[:, :]
599 tmpArray[nanIndex] = 0.
600 tmpArray[startY:endY, startX:endX] += corr[startY:endY, startX:endX]
601
602 if iteration > 0:
603 diff = numpy.sum(numpy.abs(prev_image - tmpArray))
604
605 if diff < threshold:
606 break
607 prev_image[:, :] = tmpArray[:, :]
608
609 image.getArray()[startY + 1:endY - 1, startX + 1:endX - 1] += \
610 corr[startY + 1:endY - 1, startX + 1:endX - 1]
611
612 return diff, iteration
613
614
615def transferFlux(cFunc, fStep, correctionMode=True):
616 """Take the input convolved deflection potential and the flux array
617 to compute and apply the flux transfer into the correction array.
618
619 Parameters
620 ----------
621 cFunc: `numpy.array`
622 Deflection potential, being the convolution of the flux F with the
623 kernel K.
624 fStep: `numpy.array`
625 The array of flux values which act as the source of the flux transfer.
626 correctionMode: `bool`
627 Defines if applying correction (True) or generating sims (False).
628
629 Returns
630 -------
631 corr:
632 BFE correction array
633 """
634
635 if cFunc.shape != fStep.shape:
636 raise RuntimeError(f'transferFlux: array shapes do not match: {cFunc.shape}, {fStep.shape}')
637
638 # set the sign of the correction and set its value for the
639 # time averaged solution
640 if correctionMode:
641 # negative sign if applying BFE correction
642 factor = -0.5
643 else:
644 # positive sign if generating BFE simulations
645 factor = 0.5
646
647 # initialise the BFE correction image to zero
648 corr = numpy.zeros_like(cFunc)
649
650 # Generate a 2D mesh of x,y coordinates
651 yDim, xDim = cFunc.shape
652 y = numpy.arange(yDim, dtype=int)
653 x = numpy.arange(xDim, dtype=int)
654 xc, yc = numpy.meshgrid(x, y)
655
656 # process each axis in turn
657 for ax in [0, 1]:
658
659 # gradient of phi on right/upper edge of pixel
660 diff = numpy.diff(cFunc, axis=ax)
661
662 # expand array back to full size with zero gradient at the end
663 gx = numpy.zeros_like(cFunc)
664 yDiff, xDiff = diff.shape
665 gx[:yDiff, :xDiff] += diff
666
667 # select pixels with either positive gradients on the right edge,
668 # flux flowing to the right/up
669 # or negative gradients, flux flowing to the left/down
670 for i, sel in enumerate([gx > 0, gx < 0]):
671 xSelPixels = xc[sel]
672 ySelPixels = yc[sel]
673 # and add the flux into the pixel to the right or top
674 # depending on which axis we are handling
675 if ax == 0:
676 xPix = xSelPixels
677 yPix = ySelPixels+1
678 else:
679 xPix = xSelPixels+1
680 yPix = ySelPixels
681 # define flux as the either current pixel value or pixel
682 # above/right
683 # depending on whether positive or negative gradient
684 if i == 0:
685 # positive gradients, flux flowing to higher coordinate values
686 flux = factor * fStep[sel]*gx[sel]
687 else:
688 # negative gradients, flux flowing to lower coordinate values
689 flux = factor * fStep[yPix, xPix]*gx[sel]
690 # change the fluxes of the donor and receiving pixels
691 # such that flux is conserved
692 corr[sel] -= flux
693 corr[yPix, xPix] += flux
694
695 # return correction array
696 return corr
697
698
699def fluxConservingBrighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain,
700 gains=None, correctionMode=True):
701 """Apply brighter fatter correction in place for the image.
702
703 This version presents a modified version of the algorithm
704 found in ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``
705 which conserves the image flux, resulting in improved
706 correction of the cores of stars. The convolution has also been
707 modified to mitigate edge effects.
708
709 Parameters
710 ----------
711 exposure : `lsst.afw.image.Exposure`
712 Exposure to have brighter-fatter correction applied. Modified
713 by this method.
714 kernel : `numpy.ndarray`
715 Brighter-fatter kernel to apply.
716 maxIter : scalar
717 Number of correction iterations to run.
718 threshold : scalar
719 Convergence threshold in terms of the sum of absolute
720 deviations between an iteration and the previous one.
721 applyGain : `Bool`
722 If True, then the exposure values are scaled by the gain prior
723 to correction.
724 gains : `dict` [`str`, `float`]
725 A dictionary, keyed by amplifier name, of the gains to use.
726 If gains is None, the nominal gains in the amplifier object are used.
727 correctionMode : `Bool`
728 If True (default) the function applies correction for BFE. If False,
729 the code can instead be used to generate a simulation of BFE (sign
730 change in the direction of the effect)
731
732 Returns
733 -------
734 diff : `float`
735 Final difference between iterations achieved in correction.
736 iteration : `int`
737 Number of iterations used to calculate correction.
738
739 Notes
740 -----
741 Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.
742
743 This correction takes a kernel that has been derived from flat
744 field images to redistribute the charge. The gradient of the
745 kernel is the deflection field due to the accumulated charge.
746
747 Given the original image I(x) and the kernel K(x) we can compute
748 the corrected image Ic(x) using the following equation:
749
750 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
751
752 Improved algorithm at this step applies the divergence theorem to
753 obtain a pixelised correction.
754
755 Because we use the measured counts instead of the incident counts
756 we apply the correction iteratively to reconstruct the original
757 counts and the correction. We stop iterating when the summed
758 difference between the current corrected image and the one from
759 the previous iteration is below the threshold. We do not require
760 convergence because the number of iterations is too large a
761 computational cost. How we define the threshold still needs to be
762 evaluated, the current default was shown to work reasonably well
763 on a small set of images.
764
765 Edges are handled in the convolution by padding. This is still not
766 a physical model for the edge, but avoids discontinuity in the correction.
767
768 Author of modified version: Lance.Miller@physics.ox.ac.uk
769 (see DM-38555).
770 """
771 image = exposure.getMaskedImage().getImage()
772
773 # The image needs to be units of electrons/holes
774 with gainContext(exposure, image, applyGain, gains):
775
776 # get kernel and its shape
777 kLy, kLx = kernel.shape
778 kernelImage = afwImage.ImageD(kLx, kLy)
779 kernelImage.getArray()[:, :] = kernel
780 tempImage = image.clone()
781
782 nanIndex = numpy.isnan(tempImage.getArray())
783 tempImage.getArray()[nanIndex] = 0.
784
785 outImage = afwImage.ImageF(image.getDimensions())
786 corr = numpy.zeros_like(image.getArray())
787 prevImage = numpy.zeros_like(image.getArray())
788 convCntrl = afwMath.ConvolutionControl(False, False, 1)
789 fixedKernel = afwMath.FixedKernel(kernelImage)
790
791 # set the padding amount
792 # ensure we pad by an even amount larger than the kernel
793 kLy = 2 * ((1+kLy)//2)
794 kLx = 2 * ((1+kLx)//2)
795
796 # The deflection potential only depends on the gradient of
797 # the convolution, so we can subtract the mean, which then
798 # allows us to pad the image with zeros and avoid wrap-around effects
799 # (although still not handling the image edges with a physical model)
800 # This wouldn't be great if there were a strong image gradient.
801 imYdimension, imXdimension = tempImage.array.shape
802 imean = numpy.mean(tempImage.getArray()[~nanIndex])
803 # subtract mean from image
804 tempImage -= imean
805 tempImage.array[nanIndex] = 0.
806 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
807 outImage = afwImage.ImageF(numpy.pad(outImage.getArray(), ((0, kLy), (0, kLx))))
808 # Convert array to afw image so afwMath.convolve works
809 padImage = afwImage.ImageF(padArray.shape[1], padArray.shape[0])
810 padImage.array[:] = padArray
811
812 for iteration in range(maxIter):
813
814 # create deflection potential, convolution of flux with kernel
815 # using padded counts array
816 afwMath.convolve(outImage, padImage, fixedKernel, convCntrl)
817 tmpArray = tempImage.getArray()
818 outArray = outImage.getArray()
819
820 # trim convolution output back to original shape
821 outArray = outArray[:imYdimension, :imXdimension]
822
823 # generate the correction array, with correctionMode set as input
824 corr[...] = transferFlux(outArray, tmpArray, correctionMode=correctionMode)
825
826 # update the arrays for the next iteration
827 tmpArray[:, :] = image.getArray()[:, :]
828 tmpArray += corr
829 tmpArray[nanIndex] = 0.
830 # update padded array
831 # subtract mean
832 tmpArray -= imean
833 tempImage.array[nanIndex] = 0.
834 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
835
836 if iteration > 0:
837 diff = numpy.sum(numpy.abs(prevImage - tmpArray))
838
839 if diff < threshold:
840 break
841 prevImage[:, :] = tmpArray[:, :]
842
843 image.getArray()[:] += corr[:]
844
845 return diff, iteration
846
847
848@contextmanager
849def gainContext(exp, image, apply, gains=None):
850 """Context manager that applies and removes gain.
851
852 Parameters
853 ----------
854 exp : `lsst.afw.image.Exposure`
855 Exposure to apply/remove gain.
856 image : `lsst.afw.image.Image`
857 Image to apply/remove gain.
858 apply : `Bool`
859 If True, apply and remove the amplifier gain.
860 gains : `dict` [`str`, `float`]
861 A dictionary, keyed by amplifier name, of the gains to use.
862 If gains is None, the nominal gains in the amplifier object are used.
863
864 Yields
865 ------
866 exp : `lsst.afw.image.Exposure`
867 Exposure with the gain applied.
868 """
869 # check we have all of them if provided because mixing and matching would
870 # be a real mess
871 if gains and apply is True:
872 ampNames = [amp.getName() for amp in exp.getDetector()]
873 for ampName in ampNames:
874 if ampName not in gains.keys():
875 raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
876
877 if apply:
878 ccd = exp.getDetector()
879 for amp in ccd:
880 sim = image.Factory(image, amp.getBBox())
881 if gains:
882 gain = gains[amp.getName()]
883 else:
884 gain = amp.getGain()
885 sim *= gain
886
887 try:
888 yield exp
889 finally:
890 if apply:
891 ccd = exp.getDetector()
892 for amp in ccd:
893 sim = image.Factory(image, amp.getBBox())
894 if gains:
895 gain = gains[amp.getName()]
896 else:
897 gain = amp.getGain()
898 sim /= gain
899
900
901def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
902 sensorTransmission=None, atmosphereTransmission=None):
903 """Attach a TransmissionCurve to an Exposure, given separate curves for
904 different components.
905
906 Parameters
907 ----------
908 exposure : `lsst.afw.image.Exposure`
909 Exposure object to modify by attaching the product of all given
910 ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
911 Must have a valid ``Detector`` attached that matches the detector
912 associated with sensorTransmission.
913 opticsTransmission : `lsst.afw.image.TransmissionCurve`
914 A ``TransmissionCurve`` that represents the throughput of the optics,
915 to be evaluated in focal-plane coordinates.
916 filterTransmission : `lsst.afw.image.TransmissionCurve`
917 A ``TransmissionCurve`` that represents the throughput of the filter
918 itself, to be evaluated in focal-plane coordinates.
919 sensorTransmission : `lsst.afw.image.TransmissionCurve`
920 A ``TransmissionCurve`` that represents the throughput of the sensor
921 itself, to be evaluated in post-assembly trimmed detector coordinates.
922 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
923 A ``TransmissionCurve`` that represents the throughput of the
924 atmosphere, assumed to be spatially constant.
925
926 Returns
927 -------
928 combined : `lsst.afw.image.TransmissionCurve`
929 The TransmissionCurve attached to the exposure.
930
931 Notes
932 -----
933 All ``TransmissionCurve`` arguments are optional; if none are provided, the
934 attached ``TransmissionCurve`` will have unit transmission everywhere.
935 """
936 combined = afwImage.TransmissionCurve.makeIdentity()
937 if atmosphereTransmission is not None:
938 combined *= atmosphereTransmission
939 if opticsTransmission is not None:
940 combined *= opticsTransmission
941 if filterTransmission is not None:
942 combined *= filterTransmission
943 detector = exposure.getDetector()
944 fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
945 toSys=camGeom.PIXELS)
946 combined = combined.transformedBy(fpToPix)
947 if sensorTransmission is not None:
948 combined *= sensorTransmission
949 exposure.getInfo().setTransmissionCurve(combined)
950 return combined
951
952
953def applyGains(exposure, normalizeGains=False, ptcGains=None):
954 """Scale an exposure by the amplifier gains.
955
956 Parameters
957 ----------
958 exposure : `lsst.afw.image.Exposure`
959 Exposure to process. The image is modified.
960 normalizeGains : `Bool`, optional
961 If True, then amplifiers are scaled to force the median of
962 each amplifier to equal the median of those medians.
963 ptcGains : `dict`[`str`], optional
964 Dictionary keyed by amp name containing the PTC gains.
965 """
966 ccd = exposure.getDetector()
967 ccdImage = exposure.getMaskedImage()
968
969 medians = []
970 for amp in ccd:
971 sim = ccdImage.Factory(ccdImage, amp.getBBox())
972 if ptcGains:
973 sim *= ptcGains[amp.getName()]
974 else:
975 sim *= amp.getGain()
976
977 if normalizeGains:
978 medians.append(numpy.median(sim.getImage().getArray()))
979
980 if normalizeGains:
981 median = numpy.median(numpy.array(medians))
982 for index, amp in enumerate(ccd):
983 sim = ccdImage.Factory(ccdImage, amp.getBBox())
984 if medians[index] != 0.0:
985 sim *= median/medians[index]
986
987
989 """Grow the saturation trails by an amount dependent on the width of the
990 trail.
991
992 Parameters
993 ----------
994 mask : `lsst.afw.image.Mask`
995 Mask which will have the saturated areas grown.
996 """
997
998 extraGrowDict = {}
999 for i in range(1, 6):
1000 extraGrowDict[i] = 0
1001 for i in range(6, 8):
1002 extraGrowDict[i] = 1
1003 for i in range(8, 10):
1004 extraGrowDict[i] = 3
1005 extraGrowMax = 4
1006
1007 if extraGrowMax <= 0:
1008 return
1009
1010 saturatedBit = mask.getPlaneBitMask("SAT")
1011
1012 xmin, ymin = mask.getBBox().getMin()
1013 width = mask.getWidth()
1014
1015 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1016 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1017
1018 for fp in fpList:
1019 for s in fp.getSpans():
1020 x0, x1 = s.getX0(), s.getX1()
1021
1022 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1023 if extraGrow > 0:
1024 y = s.getY() - ymin
1025 x0 -= xmin + extraGrow
1026 x1 -= xmin - extraGrow
1027
1028 if x0 < 0:
1029 x0 = 0
1030 if x1 >= width - 1:
1031 x1 = width - 1
1032
1033 mask.array[y, x0:x1+1] |= saturatedBit
1034
1035
1036def setBadRegions(exposure, badStatistic="MEDIAN"):
1037 """Set all BAD areas of the chip to the average of the rest of the exposure
1038
1039 Parameters
1040 ----------
1041 exposure : `lsst.afw.image.Exposure`
1042 Exposure to mask. The exposure mask is modified.
1043 badStatistic : `str`, optional
1044 Statistic to use to generate the replacement value from the
1045 image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
1046
1047 Returns
1048 -------
1049 badPixelCount : scalar
1050 Number of bad pixels masked.
1051 badPixelValue : scalar
1052 Value substituted for bad pixels.
1053
1054 Raises
1055 ------
1056 RuntimeError
1057 Raised if `badStatistic` is not an allowed value.
1058 """
1059 if badStatistic == "MEDIAN":
1060 statistic = afwMath.MEDIAN
1061 elif badStatistic == "MEANCLIP":
1062 statistic = afwMath.MEANCLIP
1063 else:
1064 raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
1065
1066 mi = exposure.getMaskedImage()
1067 mask = mi.getMask()
1068 BAD = mask.getPlaneBitMask("BAD")
1069 INTRP = mask.getPlaneBitMask("INTRP")
1070
1072 sctrl.setAndMask(BAD)
1073 value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
1074
1075 maskArray = mask.getArray()
1076 imageArray = mi.getImage().getArray()
1077 badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
1078 imageArray[:] = numpy.where(badPixels, value, imageArray)
1079
1080 return badPixels.sum(), value
1081
1082
1083def checkFilter(exposure, filterList, log):
1084 """Check to see if an exposure is in a filter specified by a list.
1085
1086 The goal of this is to provide a unified filter checking interface
1087 for all filter dependent stages.
1088
1089 Parameters
1090 ----------
1091 exposure : `lsst.afw.image.Exposure`
1092 Exposure to examine.
1093 filterList : `list` [`str`]
1094 List of physical_filter names to check.
1095 log : `logging.Logger`
1096 Logger to handle messages.
1097
1098 Returns
1099 -------
1100 result : `bool`
1101 True if the exposure's filter is contained in the list.
1102 """
1103 if len(filterList) == 0:
1104 return False
1105 thisFilter = exposure.getFilter()
1106 if thisFilter is None:
1107 log.warning("No FilterLabel attached to this exposure!")
1108 return False
1109
1110 thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
1111 if thisPhysicalFilter in filterList:
1112 return True
1113 elif thisFilter.bandLabel in filterList:
1114 if log:
1115 log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
1116 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1117 return True
1118 else:
1119 return False
1120
1121
1122def getPhysicalFilter(filterLabel, log):
1123 """Get the physical filter label associated with the given filterLabel.
1124
1125 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1126 associated with the given ``filterLabel``, the returned label will be
1127 "Unknown".
1128
1129 Parameters
1130 ----------
1131 filterLabel : `lsst.afw.image.FilterLabel`
1132 The `lsst.afw.image.FilterLabel` object from which to derive the
1133 physical filter label.
1134 log : `logging.Logger`
1135 Logger to handle messages.
1136
1137 Returns
1138 -------
1139 physicalFilter : `str`
1140 The value returned by the physicalLabel attribute of ``filterLabel`` if
1141 it exists, otherwise set to \"Unknown\".
1142 """
1143 if filterLabel is None:
1144 physicalFilter = "Unknown"
1145 log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
1146 else:
1147 try:
1148 physicalFilter = filterLabel.physicalLabel
1149 except RuntimeError:
1150 log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
1151 physicalFilter = "Unknown"
1152 return physicalFilter
1153
1154
1155def countMaskedPixels(maskedIm, maskPlane):
1156 """Count the number of pixels in a given mask plane.
1157
1158 Parameters
1159 ----------
1160 maskedIm : `~lsst.afw.image.MaskedImage`
1161 Masked image to examine.
1162 maskPlane : `str`
1163 Name of the mask plane to examine.
1164
1165 Returns
1166 -------
1167 nPix : `int`
1168 Number of pixels in the requested mask plane.
1169 """
1170 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1171 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1172 return nPix
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)
setBadRegions(exposure, badStatistic="MEDIAN")
countMaskedPixels(maskedIm, maskPlane)
attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None, sensorTransmission=None, atmosphereTransmission=None)
updateVariance(maskedImage, gain, readNoise)
applyGains(exposure, normalizeGains=False, ptcGains=None)
flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False)
interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1, maskNameList=['SAT'], fallbackValue=None)
illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True)
growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD")
gainContext(exp, image, apply, gains=None)
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)
interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None)
makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT')
transposeMaskedImage(maskedImage)
getPhysicalFilter(filterLabel, log)
saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT', fallbackValue=None)
trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage)