LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
isrFunctions.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import math
23 import numpy
24 
25 import lsst.geom
26 import lsst.afw.image as afwImage
27 import lsst.afw.detection as afwDetection
28 import lsst.afw.math as afwMath
29 import lsst.meas.algorithms as measAlg
30 import lsst.afw.cameraGeom as camGeom
31 
32 from lsst.meas.algorithms.detection import SourceDetectionTask
33 
34 from contextlib import contextmanager
35 
36 from .overscan import OverscanCorrectionTask, OverscanCorrectionTaskConfig
37 from .defects import Defects
38 
39 
40 def createPsf(fwhm):
41  """Make a double Gaussian PSF.
42 
43  Parameters
44  ----------
45  fwhm : scalar
46  FWHM of double Gaussian smoothing kernel.
47 
48  Returns
49  -------
50  psf : `lsst.meas.algorithms.DoubleGaussianPsf`
51  The created smoothing kernel.
52  """
53  ksize = 4*int(fwhm) + 1
54  return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
55 
56 
57 def transposeMaskedImage(maskedImage):
58  """Make a transposed copy of a masked image.
59 
60  Parameters
61  ----------
62  maskedImage : `lsst.afw.image.MaskedImage`
63  Image to process.
64 
65  Returns
66  -------
67  transposed : `lsst.afw.image.MaskedImage`
68  The transposed copy of the input image.
69  """
70  transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
71  transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
72  transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
73  transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
74  return transposed
75 
76 
77 def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None):
78  """Interpolate over defects specified in a defect list.
79 
80  Parameters
81  ----------
82  maskedImage : `lsst.afw.image.MaskedImage`
83  Image to process.
84  defectList : `lsst.meas.algorithms.Defects`
85  List of defects to interpolate over.
86  fwhm : scalar
87  FWHM of double Gaussian smoothing kernel.
88  fallbackValue : scalar, optional
89  Fallback value if an interpolated value cannot be determined.
90  If None, then the clipped mean of the image is used.
91  """
92  psf = createPsf(fwhm)
93  if fallbackValue is None:
94  fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
95  if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
96  maskedImage.getMask().addMaskPlane('INTRP')
97  measAlg.interpolateOverDefects(maskedImage, psf, defectList, fallbackValue, True)
98  return maskedImage
99 
100 
101 def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
102  """Mask pixels based on threshold detection.
103 
104  Parameters
105  ----------
106  maskedImage : `lsst.afw.image.MaskedImage`
107  Image to process. Only the mask plane is updated.
108  threshold : scalar
109  Detection threshold.
110  growFootprints : scalar, optional
111  Number of pixels to grow footprints of detected regions.
112  maskName : str, optional
113  Mask plane name, or list of names to convert
114 
115  Returns
116  -------
117  defectList : `lsst.meas.algorithms.Defects`
118  Defect list constructed from pixels above the threshold.
119  """
120  # find saturated regions
121  thresh = afwDetection.Threshold(threshold)
122  fs = afwDetection.FootprintSet(maskedImage, thresh)
123 
124  if growFootprints > 0:
125  fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
126  fpList = fs.getFootprints()
127 
128  # set mask
129  mask = maskedImage.getMask()
130  bitmask = mask.getPlaneBitMask(maskName)
131  afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
132 
133  return Defects.fromFootprintList(fpList)
134 
135 
136 def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
137  """Grow a mask by an amount and add to the requested plane.
138 
139  Parameters
140  ----------
141  mask : `lsst.afw.image.Mask`
142  Mask image to process.
143  radius : scalar
144  Amount to grow the mask.
145  maskNameList : `str` or `list` [`str`]
146  Mask names that should be grown.
147  maskValue : `str`
148  Mask plane to assign the newly masked pixels to.
149  """
150  if radius > 0:
151  thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
152  fpSet = afwDetection.FootprintSet(mask, thresh)
153  fpSet = afwDetection.FootprintSet(fpSet, rGrow=radius, isotropic=False)
154  fpSet.setMask(mask, maskValue)
155 
156 
157 def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1,
158  maskNameList=['SAT'], fallbackValue=None):
159  """Interpolate over defects identified by a particular set of mask planes.
160 
161  Parameters
162  ----------
163  maskedImage : `lsst.afw.image.MaskedImage`
164  Image to process.
165  fwhm : scalar
166  FWHM of double Gaussian smoothing kernel.
167  growSaturatedFootprints : scalar, optional
168  Number of pixels to grow footprints for saturated pixels.
169  maskNameList : `List` of `str`, optional
170  Mask plane name.
171  fallbackValue : scalar, optional
172  Value of last resort for interpolation.
173  """
174  mask = maskedImage.getMask()
175 
176  if growSaturatedFootprints > 0 and "SAT" in maskNameList:
177  # If we are interpolating over an area larger than the original masked region, we need
178  # to expand the original mask bit to the full area to explain why we interpolated there.
179  growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
180 
181  thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
182  fpSet = afwDetection.FootprintSet(mask, thresh)
183  defectList = Defects.fromFootprintList(fpSet.getFootprints())
184 
185  interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue)
186 
187  return maskedImage
188 
189 
190 def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
191  fallbackValue=None):
192  """Mark saturated pixels and optionally interpolate over them
193 
194  Parameters
195  ----------
196  maskedImage : `lsst.afw.image.MaskedImage`
197  Image to process.
198  saturation : scalar
199  Saturation level used as the detection threshold.
200  fwhm : scalar
201  FWHM of double Gaussian smoothing kernel.
202  growFootprints : scalar, optional
203  Number of pixels to grow footprints of detected regions.
204  interpolate : Bool, optional
205  If True, saturated pixels are interpolated over.
206  maskName : str, optional
207  Mask plane name.
208  fallbackValue : scalar, optional
209  Value of last resort for interpolation.
210  """
211  defectList = makeThresholdMask(
212  maskedImage=maskedImage,
213  threshold=saturation,
214  growFootprints=growFootprints,
215  maskName=maskName,
216  )
217  if interpolate:
218  interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue)
219 
220  return maskedImage
221 
222 
223 def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
224  """Compute number of edge trim pixels to match the calibration data.
225 
226  Use the dimension difference between the raw exposure and the
227  calibration exposure to compute the edge trim pixels. This trim
228  is applied symmetrically, with the same number of pixels masked on
229  each side.
230 
231  Parameters
232  ----------
233  rawMaskedImage : `lsst.afw.image.MaskedImage`
234  Image to trim.
235  calibMaskedImage : `lsst.afw.image.MaskedImage`
236  Calibration image to draw new bounding box from.
237 
238  Returns
239  -------
240  replacementMaskedImage : `lsst.afw.image.MaskedImage`
241  ``rawMaskedImage`` trimmed to the appropriate size
242  Raises
243  ------
244  RuntimeError
245  Rasied if ``rawMaskedImage`` cannot be symmetrically trimmed to
246  match ``calibMaskedImage``.
247  """
248  nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
249  if nx != ny:
250  raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
251  if nx % 2 != 0:
252  raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
253  if nx < 0:
254  raise RuntimeError("Calibration maskedImage is larger than raw data.")
255 
256  nEdge = nx//2
257  if nEdge > 0:
258  replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
259  SourceDetectionTask.setEdgeBits(
260  rawMaskedImage,
261  replacementMaskedImage.getBBox(),
262  rawMaskedImage.getMask().getPlaneBitMask("EDGE")
263  )
264  else:
265  replacementMaskedImage = rawMaskedImage
266 
267  return replacementMaskedImage
268 
269 
270 def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
271  """Apply bias correction in place.
272 
273  Parameters
274  ----------
275  maskedImage : `lsst.afw.image.MaskedImage`
276  Image to process. The image is modified by this method.
277  biasMaskedImage : `lsst.afw.image.MaskedImage`
278  Bias image of the same size as ``maskedImage``
279  trimToFit : `Bool`, optional
280  If True, raw data is symmetrically trimmed to match
281  calibration size.
282 
283  Raises
284  ------
285  RuntimeError
286  Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
287  the same size.
288 
289  """
290  if trimToFit:
291  maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
292 
293  if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
294  raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
295  (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
296  maskedImage -= biasMaskedImage
297 
298 
299 def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
300  """Apply dark correction in place.
301 
302  Parameters
303  ----------
304  maskedImage : `lsst.afw.image.MaskedImage`
305  Image to process. The image is modified by this method.
306  darkMaskedImage : `lsst.afw.image.MaskedImage`
307  Dark image of the same size as ``maskedImage``.
308  expScale : scalar
309  Dark exposure time for ``maskedImage``.
310  darkScale : scalar
311  Dark exposure time for ``darkMaskedImage``.
312  invert : `Bool`, optional
313  If True, re-add the dark to an already corrected image.
314  trimToFit : `Bool`, optional
315  If True, raw data is symmetrically trimmed to match
316  calibration size.
317 
318  Raises
319  ------
320  RuntimeError
321  Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
322  the same size.
323 
324  Notes
325  -----
326  The dark correction is applied by calculating:
327  maskedImage -= dark * expScaling / darkScaling
328  """
329  if trimToFit:
330  maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
331 
332  if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
333  raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
334  (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
335 
336  scale = expScale / darkScale
337  if not invert:
338  maskedImage.scaledMinus(scale, darkMaskedImage)
339  else:
340  maskedImage.scaledPlus(scale, darkMaskedImage)
341 
342 
343 def updateVariance(maskedImage, gain, readNoise):
344  """Set the variance plane based on the image plane.
345 
346  Parameters
347  ----------
348  maskedImage : `lsst.afw.image.MaskedImage`
349  Image to process. The variance plane is modified.
350  gain : scalar
351  The amplifier gain in electrons/ADU.
352  readNoise : scalar
353  The amplifier read nmoise in ADU/pixel.
354  """
355  var = maskedImage.getVariance()
356  var[:] = maskedImage.getImage()
357  var /= gain
358  var += readNoise**2
359 
360 
361 def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
362  """Apply flat correction in place.
363 
364  Parameters
365  ----------
366  maskedImage : `lsst.afw.image.MaskedImage`
367  Image to process. The image is modified.
368  flatMaskedImage : `lsst.afw.image.MaskedImage`
369  Flat image of the same size as ``maskedImage``
370  scalingType : str
371  Flat scale computation method. Allowed values are 'MEAN',
372  'MEDIAN', or 'USER'.
373  userScale : scalar, optional
374  Scale to use if ``scalingType``='USER'.
375  invert : `Bool`, optional
376  If True, unflatten an already flattened image.
377  trimToFit : `Bool`, optional
378  If True, raw data is symmetrically trimmed to match
379  calibration size.
380 
381  Raises
382  ------
383  RuntimeError
384  Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
385  the same size or if ``scalingType`` is not an allowed value.
386  """
387  if trimToFit:
388  maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
389 
390  if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
391  raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
392  (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
393 
394  # Figure out scale from the data
395  # Ideally the flats are normalized by the calibration product pipeline, but this allows some flexibility
396  # in the case that the flat is created by some other mechanism.
397  if scalingType in ('MEAN', 'MEDIAN'):
398  scalingType = afwMath.stringToStatisticsProperty(scalingType)
399  flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
400  elif scalingType == 'USER':
401  flatScale = userScale
402  else:
403  raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
404 
405  if not invert:
406  maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
407  else:
408  maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
409 
410 
411 def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
412  """Apply illumination correction in place.
413 
414  Parameters
415  ----------
416  maskedImage : `lsst.afw.image.MaskedImage`
417  Image to process. The image is modified.
418  illumMaskedImage : `lsst.afw.image.MaskedImage`
419  Illumination correction image of the same size as ``maskedImage``.
420  illumScale : scalar
421  Scale factor for the illumination correction.
422  trimToFit : `Bool`, optional
423  If True, raw data is symmetrically trimmed to match
424  calibration size.
425 
426  Raises
427  ------
428  RuntimeError
429  Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
430  the same size.
431  """
432  if trimToFit:
433  maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
434 
435  if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
436  raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
437  (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
438 
439  maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
440 
441 
442 def overscanCorrection(ampMaskedImage, overscanImage, fitType='MEDIAN', order=1, collapseRej=3.0,
443  statControl=None, overscanIsInt=True):
444  """Apply overscan correction in place.
445 
446  Parameters
447  ----------
448  ampMaskedImage : `lsst.afw.image.MaskedImage`
449  Image of amplifier to correct; modified.
450  overscanImage : `lsst.afw.image.Image` or `lsst.afw.image.MaskedImage`
451  Image of overscan; modified.
452  fitType : `str`
453  Type of fit for overscan correction. May be one of:
454 
455  - ``MEAN``: use mean of overscan.
456  - ``MEANCLIP``: use clipped mean of overscan.
457  - ``MEDIAN``: use median of overscan.
458  - ``MEDIAN_PER_ROW``: use median per row of overscan.
459  - ``POLY``: fit with ordinary polynomial.
460  - ``CHEB``: fit with Chebyshev polynomial.
461  - ``LEG``: fit with Legendre polynomial.
462  - ``NATURAL_SPLINE``: fit with natural spline.
463  - ``CUBIC_SPLINE``: fit with cubic spline.
464  - ``AKIMA_SPLINE``: fit with Akima spline.
465 
466  order : `int`
467  Polynomial order or number of spline knots; ignored unless
468  ``fitType`` indicates a polynomial or spline.
469  statControl : `lsst.afw.math.StatisticsControl`
470  Statistics control object. In particular, we pay attention to numSigmaClip
471  overscanIsInt : `bool`
472  Treat the overscan region as consisting of integers, even if it's been
473  converted to float. E.g. handle ties properly.
474 
475  Returns
476  -------
477  result : `lsst.pipe.base.Struct`
478  Result struct with components:
479 
480  - ``imageFit``: Value(s) removed from image (scalar or
481  `lsst.afw.image.Image`)
482  - ``overscanFit``: Value(s) removed from overscan (scalar or
483  `lsst.afw.image.Image`)
484  - ``overscanImage``: Overscan corrected overscan region
485  (`lsst.afw.image.Image`)
486  Raises
487  ------
488  RuntimeError
489  Raised if ``fitType`` is not an allowed value.
490 
491  Notes
492  -----
493  The ``ampMaskedImage`` and ``overscanImage`` are modified, with the fit
494  subtracted. Note that the ``overscanImage`` should not be a subimage of
495  the ``ampMaskedImage``, to avoid being subtracted twice.
496 
497  Debug plots are available for the SPLINE fitTypes by setting the
498  `debug.display` for `name` == "lsst.ip.isr.isrFunctions". These
499  plots show the scatter plot of the overscan data (collapsed along
500  the perpendicular dimension) as a function of position on the CCD
501  (normalized between +/-1).
502  """
503  ampImage = ampMaskedImage.getImage()
504 
506  if fitType:
507  config.fitType = fitType
508  if order:
509  config.order = order
510  if collapseRej:
511  config.numSigmaClip = collapseRej
512  if overscanIsInt:
513  config.overscanIsInt = True
514 
515  overscanTask = OverscanCorrectionTask(config=config)
516  return overscanTask.run(ampImage, overscanImage)
517 
518 
519 def brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None):
520  """Apply brighter fatter correction in place for the image.
521 
522  Parameters
523  ----------
524  exposure : `lsst.afw.image.Exposure`
525  Exposure to have brighter-fatter correction applied. Modified
526  by this method.
527  kernel : `numpy.ndarray`
528  Brighter-fatter kernel to apply.
529  maxIter : scalar
530  Number of correction iterations to run.
531  threshold : scalar
532  Convergence threshold in terms of the sum of absolute
533  deviations between an iteration and the previous one.
534  applyGain : `Bool`
535  If True, then the exposure values are scaled by the gain prior
536  to correction.
537  gains : `dict` [`str`, `float`]
538  A dictionary, keyed by amplifier name, of the gains to use.
539  If gains is None, the nominal gains in the amplifier object are used.
540 
541  Returns
542  -------
543  diff : `float`
544  Final difference between iterations achieved in correction.
545  iteration : `int`
546  Number of iterations used to calculate correction.
547 
548  Notes
549  -----
550  This correction takes a kernel that has been derived from flat
551  field images to redistribute the charge. The gradient of the
552  kernel is the deflection field due to the accumulated charge.
553 
554  Given the original image I(x) and the kernel K(x) we can compute
555  the corrected image Ic(x) using the following equation:
556 
557  Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
558 
559  To evaluate the derivative term we expand it as follows:
560 
561  0.5 * ( d/dx(I(x))*d/dx(int(dy*K(x-y)*I(y))) + I(x)*d^2/dx^2(int(dy* K(x-y)*I(y))) )
562 
563  Because we use the measured counts instead of the incident counts
564  we apply the correction iteratively to reconstruct the original
565  counts and the correction. We stop iterating when the summed
566  difference between the current corrected image and the one from
567  the previous iteration is below the threshold. We do not require
568  convergence because the number of iterations is too large a
569  computational cost. How we define the threshold still needs to be
570  evaluated, the current default was shown to work reasonably well
571  on a small set of images. For more information on the method see
572  DocuShare Document-19407.
573 
574  The edges as defined by the kernel are not corrected because they
575  have spurious values due to the convolution.
576  """
577  image = exposure.getMaskedImage().getImage()
578 
579  # The image needs to be units of electrons/holes
580  with gainContext(exposure, image, applyGain, gains):
581 
582  kLx = numpy.shape(kernel)[0]
583  kLy = numpy.shape(kernel)[1]
584  kernelImage = afwImage.ImageD(kLx, kLy)
585  kernelImage.getArray()[:, :] = kernel
586  tempImage = image.clone()
587 
588  nanIndex = numpy.isnan(tempImage.getArray())
589  tempImage.getArray()[nanIndex] = 0.
590 
591  outImage = afwImage.ImageF(image.getDimensions())
592  corr = numpy.zeros_like(image.getArray())
593  prev_image = numpy.zeros_like(image.getArray())
594  convCntrl = afwMath.ConvolutionControl(False, True, 1)
595  fixedKernel = afwMath.FixedKernel(kernelImage)
596 
597  # Define boundary by convolution region. The region that the correction will be
598  # calculated for is one fewer in each dimension because of the second derivative terms.
599  # NOTE: these need to use integer math, as we're using start:end as numpy index ranges.
600  startX = kLx//2
601  endX = -kLx//2
602  startY = kLy//2
603  endY = -kLy//2
604 
605  for iteration in range(maxIter):
606 
607  afwMath.convolve(outImage, tempImage, fixedKernel, convCntrl)
608  tmpArray = tempImage.getArray()
609  outArray = outImage.getArray()
610 
611  with numpy.errstate(invalid="ignore", over="ignore"):
612  # First derivative term
613  gradTmp = numpy.gradient(tmpArray[startY:endY, startX:endX])
614  gradOut = numpy.gradient(outArray[startY:endY, startX:endX])
615  first = (gradTmp[0]*gradOut[0] + gradTmp[1]*gradOut[1])[1:-1, 1:-1]
616 
617  # Second derivative term
618  diffOut20 = numpy.diff(outArray, 2, 0)[startY:endY, startX + 1:endX - 1]
619  diffOut21 = numpy.diff(outArray, 2, 1)[startY + 1:endY - 1, startX:endX]
620  second = tmpArray[startY + 1:endY - 1, startX + 1:endX - 1]*(diffOut20 + diffOut21)
621 
622  corr[startY + 1:endY - 1, startX + 1:endX - 1] = 0.5*(first + second)
623 
624  tmpArray[:, :] = image.getArray()[:, :]
625  tmpArray[nanIndex] = 0.
626  tmpArray[startY:endY, startX:endX] += corr[startY:endY, startX:endX]
627 
628  if iteration > 0:
629  diff = numpy.sum(numpy.abs(prev_image - tmpArray))
630 
631  if diff < threshold:
632  break
633  prev_image[:, :] = tmpArray[:, :]
634 
635  image.getArray()[startY + 1:endY - 1, startX + 1:endX - 1] += \
636  corr[startY + 1:endY - 1, startX + 1:endX - 1]
637 
638  return diff, iteration
639 
640 
641 @contextmanager
642 def gainContext(exp, image, apply, gains=None):
643  """Context manager that applies and removes gain.
644 
645  Parameters
646  ----------
647  exp : `lsst.afw.image.Exposure`
648  Exposure to apply/remove gain.
649  image : `lsst.afw.image.Image`
650  Image to apply/remove gain.
651  apply : `Bool`
652  If True, apply and remove the amplifier gain.
653  gains : `dict` [`str`, `float`]
654  A dictionary, keyed by amplifier name, of the gains to use.
655  If gains is None, the nominal gains in the amplifier object are used.
656 
657  Yields
658  ------
659  exp : `lsst.afw.image.Exposure`
660  Exposure with the gain applied.
661  """
662  # check we have all of them if provided because mixing and matching would
663  # be a real mess
664  if gains and apply is True:
665  ampNames = [amp.getName() for amp in exp.getDetector()]
666  for ampName in ampNames:
667  if ampName not in gains.keys():
668  raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
669 
670  if apply:
671  ccd = exp.getDetector()
672  for amp in ccd:
673  sim = image.Factory(image, amp.getBBox())
674  if gains:
675  gain = gains[amp.getName()]
676  else:
677  gain = amp.getGain()
678  sim *= gain
679 
680  try:
681  yield exp
682  finally:
683  if apply:
684  ccd = exp.getDetector()
685  for amp in ccd:
686  sim = image.Factory(image, amp.getBBox())
687  if gains:
688  gain = gains[amp.getName()]
689  else:
690  gain = amp.getGain()
691  sim /= gain
692 
693 
694 def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None,
695  sensorTransmission=None, atmosphereTransmission=None):
696  """Attach a TransmissionCurve to an Exposure, given separate curves for
697  different components.
698 
699  Parameters
700  ----------
701  exposure : `lsst.afw.image.Exposure`
702  Exposure object to modify by attaching the product of all given
703  ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
704  Must have a valid ``Detector`` attached that matches the detector
705  associated with sensorTransmission.
706  opticsTransmission : `lsst.afw.image.TransmissionCurve`
707  A ``TransmissionCurve`` that represents the throughput of the optics,
708  to be evaluated in focal-plane coordinates.
709  filterTransmission : `lsst.afw.image.TransmissionCurve`
710  A ``TransmissionCurve`` that represents the throughput of the filter
711  itself, to be evaluated in focal-plane coordinates.
712  sensorTransmission : `lsst.afw.image.TransmissionCurve`
713  A ``TransmissionCurve`` that represents the throughput of the sensor
714  itself, to be evaluated in post-assembly trimmed detector coordinates.
715  atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
716  A ``TransmissionCurve`` that represents the throughput of the
717  atmosphere, assumed to be spatially constant.
718 
719  Returns
720  -------
721  combined : `lsst.afw.image.TransmissionCurve`
722  The TransmissionCurve attached to the exposure.
723 
724  Notes
725  -----
726  All ``TransmissionCurve`` arguments are optional; if none are provided, the
727  attached ``TransmissionCurve`` will have unit transmission everywhere.
728  """
729  combined = afwImage.TransmissionCurve.makeIdentity()
730  if atmosphereTransmission is not None:
731  combined *= atmosphereTransmission
732  if opticsTransmission is not None:
733  combined *= opticsTransmission
734  if filterTransmission is not None:
735  combined *= filterTransmission
736  detector = exposure.getDetector()
737  fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
738  toSys=camGeom.PIXELS)
739  combined = combined.transformedBy(fpToPix)
740  if sensorTransmission is not None:
741  combined *= sensorTransmission
742  exposure.getInfo().setTransmissionCurve(combined)
743  return combined
744 
745 
746 def applyGains(exposure, normalizeGains=False, ptcGains=None):
747  """Scale an exposure by the amplifier gains.
748 
749  Parameters
750  ----------
751  exposure : `lsst.afw.image.Exposure`
752  Exposure to process. The image is modified.
753  normalizeGains : `Bool`, optional
754  If True, then amplifiers are scaled to force the median of
755  each amplifier to equal the median of those medians.
756  ptcGains : `dict`[`str`], optional
757  Dictionary keyed by amp name containing the PTC gains.
758  """
759  ccd = exposure.getDetector()
760  ccdImage = exposure.getMaskedImage()
761 
762  medians = []
763  for amp in ccd:
764  sim = ccdImage.Factory(ccdImage, amp.getBBox())
765  if ptcGains:
766  sim *= ptcGains[amp.getName()]
767  else:
768  sim *= amp.getGain()
769 
770  if normalizeGains:
771  medians.append(numpy.median(sim.getImage().getArray()))
772 
773  if normalizeGains:
774  median = numpy.median(numpy.array(medians))
775  for index, amp in enumerate(ccd):
776  sim = ccdImage.Factory(ccdImage, amp.getBBox())
777  if medians[index] != 0.0:
778  sim *= median/medians[index]
779 
780 
782  """Grow the saturation trails by an amount dependent on the width of the trail.
783 
784  Parameters
785  ----------
786  mask : `lsst.afw.image.Mask`
787  Mask which will have the saturated areas grown.
788  """
789 
790  extraGrowDict = {}
791  for i in range(1, 6):
792  extraGrowDict[i] = 0
793  for i in range(6, 8):
794  extraGrowDict[i] = 1
795  for i in range(8, 10):
796  extraGrowDict[i] = 3
797  extraGrowMax = 4
798 
799  if extraGrowMax <= 0:
800  return
801 
802  saturatedBit = mask.getPlaneBitMask("SAT")
803 
804  xmin, ymin = mask.getBBox().getMin()
805  width = mask.getWidth()
806 
807  thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
808  fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
809 
810  for fp in fpList:
811  for s in fp.getSpans():
812  x0, x1 = s.getX0(), s.getX1()
813 
814  extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
815  if extraGrow > 0:
816  y = s.getY() - ymin
817  x0 -= xmin + extraGrow
818  x1 -= xmin - extraGrow
819 
820  if x0 < 0:
821  x0 = 0
822  if x1 >= width - 1:
823  x1 = width - 1
824 
825  mask.array[y, x0:x1+1] |= saturatedBit
826 
827 
828 def setBadRegions(exposure, badStatistic="MEDIAN"):
829  """Set all BAD areas of the chip to the average of the rest of the exposure
830 
831  Parameters
832  ----------
833  exposure : `lsst.afw.image.Exposure`
834  Exposure to mask. The exposure mask is modified.
835  badStatistic : `str`, optional
836  Statistic to use to generate the replacement value from the
837  image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
838 
839  Returns
840  -------
841  badPixelCount : scalar
842  Number of bad pixels masked.
843  badPixelValue : scalar
844  Value substituted for bad pixels.
845 
846  Raises
847  ------
848  RuntimeError
849  Raised if `badStatistic` is not an allowed value.
850  """
851  if badStatistic == "MEDIAN":
852  statistic = afwMath.MEDIAN
853  elif badStatistic == "MEANCLIP":
854  statistic = afwMath.MEANCLIP
855  else:
856  raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
857 
858  mi = exposure.getMaskedImage()
859  mask = mi.getMask()
860  BAD = mask.getPlaneBitMask("BAD")
861  INTRP = mask.getPlaneBitMask("INTRP")
862 
863  sctrl = afwMath.StatisticsControl()
864  sctrl.setAndMask(BAD)
865  value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
866 
867  maskArray = mask.getArray()
868  imageArray = mi.getImage().getArray()
869  badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
870  imageArray[:] = numpy.where(badPixels, value, imageArray)
871 
872  return badPixels.sum(), value
873 
874 
875 def checkFilter(exposure, filterList, log):
876  """Check to see if an exposure is in a filter specified by a list.
877 
878  The goal of this is to provide a unified filter checking interface
879  for all filter dependent stages.
880 
881  Parameters
882  ----------
883  exposure : `lsst.afw.image.Exposure`
884  Exposure to examine.
885  filterList : `list` [`str`]
886  List of physical_filter names to check.
887  log : `logging.Logger`
888  Logger to handle messages.
889 
890  Returns
891  -------
892  result : `bool`
893  True if the exposure's filter is contained in the list.
894  """
895  thisFilter = exposure.getFilterLabel()
896  if thisFilter is None:
897  log.warning("No FilterLabel attached to this exposure!")
898  return False
899 
900  thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
901  if thisPhysicalFilter in filterList:
902  return True
903  elif thisFilter.bandLabel in filterList:
904  if log:
905  log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
906  " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
907  return True
908  else:
909  return False
910 
911 
912 def getPhysicalFilter(filterLabel, log):
913  """Get the physical filter label associated with the given filterLabel.
914 
915  If ``filterLabel`` is `None` or there is no physicalLabel attribute
916  associated with the given ``filterLabel``, the returned label will be
917  "Unknown".
918 
919  Parameters
920  ----------
921  filterLabel : `lsst.afw.image.FilterLabel`
922  The `lsst.afw.image.FilterLabel` object from which to derive the
923  physical filter label.
924  log : `logging.Logger`
925  Logger to handle messages.
926 
927  Returns
928  -------
929  physicalFilter : `str`
930  The value returned by the physicalLabel attribute of ``filterLabel`` if
931  it exists, otherwise set to \"Unknown\".
932  """
933  if filterLabel is None:
934  physicalFilter = "Unknown"
935  log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
936  else:
937  try:
938  physicalFilter = filterLabel.physicalLabel
939  except RuntimeError:
940  log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
941  physicalFilter = "Unknown"
942  return physicalFilter
Parameters to control convolution.
Definition: ConvolveImage.h:50
A kernel created from an Image.
Definition: Kernel.h:471
Pass parameters to a Statistics object.
Definition: Statistics.h:92
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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:359
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)
Definition: Statistics.cc:738
def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False)
def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD")
def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False)
def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT', fallbackValue=None)
def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None)
Definition: isrFunctions.py:77
def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True)
def checkFilter(exposure, filterList, log)
def gainContext(exp, image, apply, gains=None)
def brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None)
def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT')
def setBadRegions(exposure, badStatistic="MEDIAN")
def applyGains(exposure, normalizeGains=False, ptcGains=None)
def getPhysicalFilter(filterLabel, log)
def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage)
def updateVariance(maskedImage, gain, readNoise)
def overscanCorrection(ampMaskedImage, overscanImage, fitType='MEDIAN', order=1, collapseRej=3.0, statControl=None, overscanIsInt=True)
def interpolateFromMask(maskedImage, fwhm, growSaturatedFootprints=1, maskNameList=['SAT'], fallbackValue=None)
def transposeMaskedImage(maskedImage)
Definition: isrFunctions.py:57
def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False)
def attachTransmissionCurve(exposure, opticsTransmission=None, filterTransmission=None, sensorTransmission=None, atmosphereTransmission=None)
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.