25 "attachTransmissionCurve",
27 "brighterFatterCorrection",
33 "fluxConservingBrighterFatterCorrection",
38 "illuminationCorrection",
39 "interpolateDefectList",
40 "interpolateFromMask",
42 "saturationCorrection",
45 "transposeMaskedImage",
46 "trimToMatchCalibBBox",
48 "widenSaturationTrails",
50 "getExposureReadNoises",
65from contextlib
import contextmanager
67from .defects
import Defects
71 """Make a double Gaussian PSF.
76 FWHM of double Gaussian smoothing kernel.
80 psf : `lsst.meas.algorithms.DoubleGaussianPsf`
81 The created smoothing kernel.
83 ksize = 4*int(fwhm) + 1
84 return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
88 """Make a transposed copy of a masked image.
92 maskedImage : `lsst.afw.image.MaskedImage`
97 transposed : `lsst.afw.image.MaskedImage`
98 The transposed copy of the input image.
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
108 maskNameList=None, useLegacyInterp=True):
109 """Interpolate over defects specified in a defect list.
113 maskedImage : `lsst.afw.image.MaskedImage`
115 defectList : `lsst.meas.algorithms.Defects`
116 List of defects to interpolate over.
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.
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.
136 if fallbackValue
is None:
138 if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
139 maskedImage.getMask().addMaskPlane(
'INTRP')
150 kwargs = {
"bin_spacing": 20,
151 "threshold_dynamic_binning": 2000,
152 "threshold_subdivide": 20000}
155 measAlg.interpolateOverDefects(maskedImage, psf, defectList,
156 fallbackValue=fallbackValue,
157 useFallbackValueAtEdge=
True,
159 useLegacyInterp=useLegacyInterp,
160 maskNameList=maskNameList, **kwargs)
165 """Mask pixels based on threshold detection.
169 maskedImage : `lsst.afw.image.MaskedImage`
170 Image to process. Only the mask plane is updated.
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
180 defectList : `lsst.meas.algorithms.Defects`
181 Defect list constructed from pixels above the threshold.
184 thresh = afwDetection.Threshold(threshold)
185 fs = afwDetection.FootprintSet(maskedImage, thresh)
187 if growFootprints > 0:
188 fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=
False)
189 fpList = fs.getFootprints()
192 mask = maskedImage.getMask()
193 bitmask = mask.getPlaneBitMask(maskName)
194 afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
196 return Defects.fromFootprintList(fpList)
199def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
200 """Grow a mask by an amount and add to the requested plane.
204 mask : `lsst.afw.image.Mask`
205 Mask image to process.
207 Amount to grow the mask.
208 maskNameList : `str` or `list` [`str`]
209 Mask names that should be grown.
211 Mask plane to assign the newly masked pixels to.
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)
199def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
…
221 itlEdgeBleedSatMaxArea=100000,
222 itlEdgeBleedThreshold=5000.,
223 itlEdgeBleedModelConstant=0.03,
224 saturatedMaskName="SAT"):
225 """Mask edge bleeds in ITL detectors.
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
237 itlEdgeBleedThreshold : `float`, optional
238 Threshold above median sky background for edge bleed detection
240 itlEdgeBleedModelConstant : `float`, optional
241 Constant in the decaying exponential in the edge bleed masking.
242 saturatedMaskName : `str`, optional
243 Mask name for saturation.
245 maskedImage = ccdExposure.maskedImage
246 xmax = maskedImage.image.array.shape[1]
249 satLevel = numpy.min([ccdExposure.metadata[f
"LSST ISR SATURATION LEVEL {amp.getName()}"]
250 for amp
in ccdExposure.getDetector()
if not badAmpDict[amp.getName()]])
253 thresh = afwDetection.Threshold(satLevel)
254 fpList = afwDetection.FootprintSet(maskedImage, thresh).getFootprints()
256 satAreas = numpy.asarray([fp.getArea()
for fp
in fpList])
257 largeAreas, = numpy.where((satAreas >= itlEdgeBleedSatMinArea)
258 & (satAreas < itlEdgeBleedSatMaxArea))
260 satMaskBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
262 if len(largeAreas) == 0:
265 for largeAreasIndex
in largeAreas:
268 xCore, yCore = fpList[largeAreasIndex].getCentroid()
271 yCoreFP = yCore - fpList[largeAreasIndex].getBBox().getMinY()
273 widthSat = numpy.sum(fpList[largeAreasIndex].getSpans().asArray()[int(yCoreFP), :])
275 for amp
in ccdExposure.getDetector():
281 yBox = amp.getBBox().getCenter()[1]
282 if amp.getBBox().contains(xCore, yBox):
285 ampName = amp.getName()
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:, :])
312 lowerRangeSmall = int(xCore)-5
313 upperRangeSmall = int(xCore)+5
314 if 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):
328 subImage = sliceImage[:100, :]
329 maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
335 edgeBleedHalfWidth = \
336 int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
338 lowerRange = int(xCore)-edgeBleedHalfWidth
339 upperRange = int(xCore)+edgeBleedHalfWidth
345 if upperRange > xmax:
347 sliceMask[y, lowerRange:upperRange] = satMaskBit
351 maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
352 """Interpolate over defects identified by a particular set of mask planes.
356 maskedImage : `lsst.afw.image.MaskedImage`
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
364 fallbackValue : scalar, optional
365 Value of last resort for interpolation.
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.
373 mask = maskedImage.getMask()
375 if growSaturatedFootprints > 0
and "SAT" in maskNameList:
379 growMasks(mask, radius=growSaturatedFootprints, maskNameList=[
'SAT'], maskValue=
"SAT")
381 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
382 fpSet = afwDetection.FootprintSet(mask, thresh)
383 defectList = Defects.fromFootprintList(fpSet.getFootprints())
386 maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
392 fallbackValue=None, useLegacyInterp=True):
393 """Mark saturated pixels and optionally interpolate over them
397 maskedImage : `lsst.afw.image.MaskedImage`
400 Saturation level used as the detection threshold.
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
409 fallbackValue : scalar, optional
410 Value of last resort for interpolation.
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.
419 maskedImage=maskedImage,
420 threshold=saturation,
421 growFootprints=growFootprints,
426 maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
432 """Compute number of edge trim pixels to match the calibration data.
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
441 rawMaskedImage : `lsst.afw.image.MaskedImage`
443 calibMaskedImage : `lsst.afw.image.MaskedImage`
444 Calibration image to draw new bounding box from.
448 replacementMaskedImage : `lsst.afw.image.MaskedImage`
449 ``rawMaskedImage`` trimmed to the appropriate size.
454 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
455 match ``calibMaskedImage``.
457 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
459 raise RuntimeError(
"Raw and calib maskedImages are trimmed differently in X and Y.")
461 raise RuntimeError(
"Calibration maskedImage is trimmed unevenly in X.")
463 raise RuntimeError(
"Calibration maskedImage is larger than raw data.")
467 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
468 SourceDetectionTask.setEdgeBits(
470 replacementMaskedImage.getBBox(),
471 rawMaskedImage.getMask().getPlaneBitMask(
"EDGE")
474 replacementMaskedImage = rawMaskedImage
476 return replacementMaskedImage
480 """Apply bias correction in place.
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
495 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
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
508def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
509 """Apply dark correction in place.
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``.
518 Dark exposure time for ``maskedImage``.
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
530 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
535 The dark correction is applied by calculating:
536 maskedImage -= dark * expScaling / darkScaling
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)))
545 scale = expScale / darkScale
547 maskedImage.scaledMinus(scale, darkMaskedImage)
549 maskedImage.scaledPlus(scale, darkMaskedImage)
508def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
…
553 """Set the variance plane based on the image plane.
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.
561 maskedImage : `lsst.afw.image.MaskedImage`
562 Image to process. The variance plane is modified.
564 The amplifier gain in electron/adu.
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.
571 var = maskedImage.variance
573 var[:, :] = maskedImage.image
575 var[:, :] += maskedImage.image
577 var += (readNoise/gain)**2
580def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
581 """Apply flat correction in place.
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``
590 Flat scale computation method. Allowed values are 'MEAN',
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
603 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
604 the same size or if ``scalingType`` is not an allowed value.
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)))
617 if scalingType
in (
'MEAN',
'MEDIAN'):
620 elif scalingType ==
'USER':
621 flatScale = userScale
623 raise RuntimeError(
'%s : %s not implemented' % (
"flatCorrection", scalingType))
626 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
628 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
580def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
…
632 """Apply illumination correction in place.
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``.
641 Scale factor for the illumination correction.
642 trimToFit : `Bool`, optional
643 If True, raw data is symmetrically trimmed to match
649 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
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)))
659 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
663 """Apply brighter fatter correction in place for the image.
667 exposure : `lsst.afw.image.Exposure`
668 Exposure to have brighter-fatter correction applied. Modified
670 kernel : `numpy.ndarray`
671 Brighter-fatter kernel to apply.
673 Number of correction iterations to run.
675 Convergence threshold in terms of the sum of absolute
676 deviations between an iteration and the previous one.
678 If True, then the exposure values are scaled by the gain prior
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.
687 Final difference between iterations achieved in correction.
689 Number of iterations used to calculate correction.
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.
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:
700 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
702 To evaluate the derivative term we expand it as follows:
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))) )
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.
718 The edges as defined by the kernel are not corrected because they
719 have spurious values due to the convolution.
721 image = exposure.getMaskedImage().getImage()
724 with gainContext(exposure, image, applyGain, gains):
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)
732 nanIndex = numpy.isnan(tempImage.getArray())
733 tempImage.getArray()[nanIndex] = 0.
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)
751 for iteration
in range(maxIter):
754 tmpArray = tempImage.getArray()
755 outArray = outImage.getArray()
757 with numpy.errstate(invalid=
"ignore", over=
"ignore"):
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]
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)
768 corr[startY + 1:endY - 1, startX + 1:endX - 1] = 0.5*(first + second)
770 tmpArray[:, :] = image.getArray()[:, :]
771 tmpArray[nanIndex] = 0.
772 tmpArray[startY:endY, startX:endX] += corr[startY:endY, startX:endX]
775 diff = numpy.sum(numpy.abs(prev_image - tmpArray), dtype=numpy.float64)
779 prev_image[:, :] = tmpArray[:, :]
781 image.getArray()[startY + 1:endY - 1, startX + 1:endX - 1] += \
782 corr[startY + 1:endY - 1, startX + 1:endX - 1]
784 return diff, iteration
788 """Take the input convolved deflection potential and the flux array
789 to compute and apply the flux transfer into the correction array.
794 Deflection potential, being the convolution of the flux F with the
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).
807 if cFunc.shape != fStep.shape:
808 raise RuntimeError(f
'transferFlux: array shapes do not match: {cFunc.shape}, {fStep.shape}')
820 corr = numpy.zeros(cFunc.shape, dtype=numpy.float64)
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)
832 diff = numpy.diff(cFunc, axis=ax)
835 gx = numpy.zeros(cFunc.shape, dtype=numpy.float64)
836 yDiff, xDiff = diff.shape
837 gx[:yDiff, :xDiff] += diff
842 for i, sel
in enumerate([gx > 0, gx < 0]):
858 flux = factor * fStep[sel]*gx[sel]
861 flux = factor * fStep[yPix, xPix]*gx[sel]
865 corr[yPix, xPix] += flux
872 gains=None, correctionMode=True):
873 """Apply brighter fatter correction in place for the image.
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.
883 exposure : `lsst.afw.image.Exposure`
884 Exposure to have brighter-fatter correction applied. Modified
886 kernel : `numpy.ndarray`
887 Brighter-fatter kernel to apply.
889 Number of correction iterations to run.
891 Convergence threshold in terms of the sum of absolute
892 deviations between an iteration and the previous one.
894 If True, then the exposure values are scaled by the gain prior
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)
907 Final difference between iterations achieved in correction.
909 Number of iterations used to calculate correction.
913 Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.
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.
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:
922 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
924 Improved algorithm at this step applies the divergence theorem to
925 obtain a pixelised correction.
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.
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.
940 Author of modified version: Lance.Miller@physics.ox.ac.uk
943 image = exposure.getMaskedImage().getImage()
946 with gainContext(exposure, image, applyGain, gains):
949 kLy, kLx = kernel.shape
950 kernelImage = afwImage.ImageD(kLx, kLy)
951 kernelImage.getArray()[:, :] = kernel
952 tempImage = afwImage.ImageD(image, deep=
True)
954 nanIndex = numpy.isnan(tempImage.getArray())
955 tempImage.getArray()[nanIndex] = 0.
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)
965 kLy = 2 * ((1+kLy)//2)
966 kLx = 2 * ((1+kLx)//2)
973 imYdimension, imXdimension = tempImage.array.shape
974 imean = numpy.mean(tempImage.getArray()[~nanIndex], dtype=numpy.float64)
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))))
981 padImage = afwImage.ImageD(padArray.shape[1], padArray.shape[0])
982 padImage.array[:] = padArray
984 for iteration
in range(maxIter):
989 tmpArray = tempImage.getArray()
990 outArray = outImage.getArray()
993 outArray = outArray[:imYdimension, :imXdimension]
996 corr[...] =
transferFlux(outArray, tmpArray, correctionMode=correctionMode)
999 tmpArray[:, :] = image.getArray()[:, :]
1001 tmpArray[nanIndex] = 0.
1005 tempImage.array[nanIndex] = 0.
1006 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
1009 diff = numpy.sum(numpy.abs(prevImage - tmpArray), dtype=numpy.float64)
1011 if diff < threshold:
1013 prevImage[:, :] = tmpArray[:, :]
1015 image.getArray()[:] += corr[:]
1017 return diff, iteration
1021def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
1022 """Context manager that applies and removes gain.
1026 exp : `lsst.afw.image.Exposure`
1027 Exposure to apply/remove gain.
1028 image : `lsst.afw.image.Image`
1029 Image to apply/remove gain.
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?
1042 exp : `lsst.afw.image.Exposure`
1043 Exposure with the gain applied.
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}")
1054 ccd = exp.getDetector()
1056 sim = image.Factory(image, amp.getBBox()
if isTrimmed
else amp.getRawBBox())
1058 gain = gains[amp.getName()]
1060 gain = amp.getGain()
1070 ccd = exp.getDetector()
1072 sim = image.Factory(image, amp.getBBox()
if isTrimmed
else amp.getRawBBox())
1074 gain = gains[amp.getName()]
1076 gain = amp.getGain()
1021def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
…
1084 sensorTransmission=None, atmosphereTransmission=None):
1085 """Attach a TransmissionCurve to an Exposure, given separate curves for
1086 different components.
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.
1110 combined : `lsst.afw.image.TransmissionCurve`
1111 The TransmissionCurve attached to the exposure.
1115 All ``TransmissionCurve`` arguments are optional; if none are provided, the
1116 attached ``TransmissionCurve`` will have unit transmission everywhere.
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)
1135def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
1136 """Scale an exposure by the amplifier gains.
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?
1150 ccd = exposure.getDetector()
1151 ccdImage = exposure.getMaskedImage()
1156 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1158 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1160 sim *= ptcGains[amp.getName()]
1162 sim *= amp.getGain()
1165 medians.append(numpy.median(sim.getImage().getArray()))
1168 median = numpy.median(numpy.array(medians))
1169 for index, amp
in enumerate(ccd):
1171 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1173 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1174 if medians[index] != 0.0:
1175 sim *= median/medians[index]
1135def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
…
1179 """Grow the saturation trails by an amount dependent on the width of the
1184 mask : `lsst.afw.image.Mask`
1185 Mask which will have the saturated areas grown.
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
1197 if extraGrowMax <= 0:
1200 saturatedBit = mask.getPlaneBitMask(
"SAT")
1202 xmin, ymin = mask.getBBox().getMin()
1203 width = mask.getWidth()
1205 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1206 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1209 for s
in fp.getSpans():
1210 x0, x1 = s.getX0(), s.getX1()
1212 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1215 x0 -= xmin + extraGrow
1216 x1 -= xmin - extraGrow
1223 mask.array[y, x0:x1+1] |= saturatedBit
1227 """Set all BAD areas of the chip to the average of the rest of the exposure
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'.
1239 badPixelCount : scalar
1240 Number of bad pixels masked.
1241 badPixelValue : scalar
1242 Value substituted for bad pixels.
1247 Raised if `badStatistic` is not an allowed value.
1249 if badStatistic ==
"MEDIAN":
1250 statistic = afwMath.MEDIAN
1251 elif badStatistic ==
"MEANCLIP":
1252 statistic = afwMath.MEANCLIP
1254 raise RuntimeError(
"Impossible method %s of bad region correction" % badStatistic)
1256 mi = exposure.getMaskedImage()
1258 BAD = mask.getPlaneBitMask(
"BAD")
1259 INTRP = mask.getPlaneBitMask(
"INTRP")
1262 sctrl.setAndMask(BAD)
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)
1270 return badPixels.sum(), value
1274 """Check to see if an exposure is in a filter specified by a list.
1276 The goal of this is to provide a unified filter checking interface
1277 for all filter dependent stages.
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.
1291 True if the exposure's filter is contained in the list.
1293 if len(filterList) == 0:
1295 thisFilter = exposure.getFilter()
1296 if thisFilter
is None:
1297 log.warning(
"No FilterLabel attached to this exposure!")
1301 if thisPhysicalFilter
in filterList:
1303 elif thisFilter.bandLabel
in filterList:
1305 log.warning(
"Physical filter (%s) should be used instead of band %s for filter configurations"
1306 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1313 """Get the physical filter label associated with the given filterLabel.
1315 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1316 associated with the given ``filterLabel``, the returned label will be
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.
1329 physicalFilter : `str`
1330 The value returned by the physicalLabel attribute of ``filterLabel`` if
1331 it exists, otherwise set to \"Unknown\".
1333 if filterLabel
is None:
1334 physicalFilter =
"Unknown"
1335 log.warning(
"filterLabel is None. Setting physicalFilter to \"Unknown\".")
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
1346 """Count the number of pixels in a given mask plane.
1350 maskedIm : `~lsst.afw.image.MaskedImage`
1351 Masked image to examine.
1353 Name of the mask plane to examine.
1358 Number of pixels in the requested mask plane.
1360 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1361 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1366 """Get the per-amplifier gains used for this exposure.
1370 exposure : `lsst.afw.image.Exposure`
1371 The exposure to find gains for.
1375 gains : `dict` [`str` `float`]
1376 Dictionary of gain values, keyed by amplifier name.
1377 Returns empty dict when detector is None.
1379 det = exposure.getDetector()
1383 metadata = exposure.getMetadata()
1386 ampName = amp.getName()
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]
1393 gains[ampName] = amp.getGain()
1398 """Get the per-amplifier read noise used for this exposure.
1402 exposure : `lsst.afw.image.Exposure`
1403 The exposure to find read noise for.
1407 readnoises : `dict` [`str` `float`]
1408 Dictionary of read noise values, keyed by amplifier name.
1409 Returns empty dict when detector is None.
1411 det = exposure.getDetector()
1415 metadata = exposure.getMetadata()
1418 ampName = amp.getName()
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]
1425 readnoises[ampName] = amp.getReadNoise()
1430 """Check if the unused pixels (pre-/over-scan pixels) have
1431 been trimmed from an exposure.
1435 exposure : `lsst.afw.image.Exposure`
1436 The exposure to check.
1441 True if the image is trimmed, else False.
1443 return exposure.getDetector().getBBox() == exposure.getBBox()
1447 """Check if the unused pixels (pre-/over-scan pixels) have
1448 been trimmed from an image
1452 image : `lsst.afw.image.Image`
1454 detector : `lsst.afw.cameraGeom.Detector`
1455 The detector associated with the image.
1460 True if the image is trimmed, else False.
1462 return detector.getBBox() == image.getBBox()
Parameters to control convolution.
A kernel created from an Image.
Pass parameters to a Statistics object.
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)
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)
getExposureGains(exposure)
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")
widenSaturationTrails(mask)
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')
getExposureReadNoises(exposure)
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")
isTrimmedExposure(exposure)
trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage)