LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
lsst.ip.isr.isrFunctions Namespace Reference

Functions

 createPsf (fwhm)
 
 transposeMaskedImage (maskedImage)
 
 interpolateDefectList (maskedImage, defectList, fwhm, fallbackValue=None, maskNameList=None, useLegacyInterp=True)
 
 makeThresholdMask (maskedImage, threshold, growFootprints=1, maskName='SAT')
 
 growMasks (mask, radius=0, maskNameList=['BAD'], maskValue="BAD")
 
 maskE2VEdgeBleed (exposure, e2vEdgeBleedSatMinArea=10000, e2vEdgeBleedSatMaxArea=100000, e2vEdgeBleedYMax=350, saturatedMaskName="SAT", log=None)
 
 maskITLEdgeBleed (ccdExposure, badAmpDict, fpCore, itlEdgeBleedSatMinArea=10000, itlEdgeBleedSatMaxArea=100000, itlEdgeBleedThreshold=5000., itlEdgeBleedModelConstant=0.02, saturatedMaskName="SAT", log=None)
 
 _applyMaskITLEdgeBleed (ccdExposure, xCore, satLevel, widthSat, itlEdgeBleedThreshold=5000., itlEdgeBleedModelConstant=0.03, saturatedMaskName="SAT", log=None)
 
 maskITLSatSag (ccdExposure, fpCore, saturatedMaskName="SAT")
 
 maskITLDip (exposure, detectorConfig, maskPlaneNames=["SUSPECT", "ITL_DIP"], log=None)
 
 interpolateFromMask (maskedImage, fwhm, growSaturatedFootprints=1, maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True)
 
 saturationCorrection (maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT', fallbackValue=None, useLegacyInterp=True)
 
 trimToMatchCalibBBox (rawMaskedImage, calibMaskedImage)
 
 biasCorrection (maskedImage, biasMaskedImage, trimToFit=False)
 
 darkCorrection (maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False)
 
 updateVariance (maskedImage, gain, readNoise, replace=True)
 
 flatCorrection (maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False)
 
 illuminationCorrection (maskedImage, illumMaskedImage, illumScale, trimToFit=True)
 
 brighterFatterCorrection (exposure, kernel, maxIter, threshold, applyGain, gains=None)
 
 transferFlux (cFunc, fStep, correctionMode=True)
 
 fluxConservingBrighterFatterCorrection (exposure, kernel, maxIter, threshold, applyGain, gains=None, correctionMode=True)
 
 gainContext (exp, image, apply, gains=None, invert=False, isTrimmed=True)
 
 attachTransmissionCurve (exposure, opticsTransmission=None, filterTransmission=None, sensorTransmission=None, atmosphereTransmission=None)
 
 applyGains (exposure, normalizeGains=False, ptcGains=None, isTrimmed=True)
 
 widenSaturationTrails (mask)
 
 setBadRegions (exposure, badStatistic="MEDIAN")
 
 checkFilter (exposure, filterList, log)
 
 getPhysicalFilter (filterLabel, log)
 
 countMaskedPixels (maskedIm, maskPlane)
 
 getExposureGains (exposure)
 
 getExposureReadNoises (exposure)
 
 isTrimmedExposure (exposure)
 
 isTrimmedImage (image, detector)
 
 compareCameraKeywords (doRaiseOnCalibMismatch, cameraKeywordsToCompare, exposureMetadata, calib, calibName, log=None)
 

Function Documentation

◆ _applyMaskITLEdgeBleed()

lsst.ip.isr.isrFunctions._applyMaskITLEdgeBleed ( ccdExposure,
xCore,
satLevel,
widthSat,
itlEdgeBleedThreshold = 5000.,
itlEdgeBleedModelConstant = 0.03,
saturatedMaskName = "SAT",
log = None )
protected
Apply ITL edge bleed masking model.

Parameters
----------
ccdExposure : `lsst.afw.image.Exposure`
    Exposure to apply masking to.
xCore: `int`
    X coordinate of the saturated core.
satLevel: `float`
    Minimum saturation level of the detector.
widthSat: `float`
    Width of the saturated core.
itlEdgeBleedThreshold : `float`, optional
    Threshold above median sky background for edge bleed detection
    (electron units).
itlEdgeBleedModelConstant : `float`, optional
    Constant in the decaying exponential in the edge bleed masking.
saturatedMaskName : `str`, optional
    Mask name for saturation.
log : `logging.Logger`, optional
    Logger to handle messages.

Definition at line 403 of file isrFunctions.py.

407 saturatedMaskName="SAT", log=None):
408 """Apply ITL edge bleed masking model.
409
410 Parameters
411 ----------
412 ccdExposure : `lsst.afw.image.Exposure`
413 Exposure to apply masking to.
414 xCore: `int`
415 X coordinate of the saturated core.
416 satLevel: `float`
417 Minimum saturation level of the detector.
418 widthSat: `float`
419 Width of the saturated core.
420 itlEdgeBleedThreshold : `float`, optional
421 Threshold above median sky background for edge bleed detection
422 (electron units).
423 itlEdgeBleedModelConstant : `float`, optional
424 Constant in the decaying exponential in the edge bleed masking.
425 saturatedMaskName : `str`, optional
426 Mask name for saturation.
427 log : `logging.Logger`, optional
428 Logger to handle messages.
429 """
430 log = log if log else logging.getLogger(__name__)
431
432 maskedImage = ccdExposure.maskedImage
433 xmax = maskedImage.image.array.shape[1]
434 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
435
436 for amp in ccdExposure.getDetector():
437 # Select the 2 top and bottom amplifiers around the saturated
438 # core with a potential edge bleed by selecting the amplifiers
439 # that have the same X coordinate as the saturated core.
440 # As we don't care about the Y coordinate, we set it to the
441 # center of the BBox.
442 yBox = amp.getBBox().getCenter()[1]
443 if amp.getBBox().contains(xCore, yBox):
444
445 # Get the amp name
446 ampName = amp.getName()
447
448 # Because in ITLs the edge bleed happens on both edges
449 # of the detector, we make a cutout around
450 # both the top and bottom
451 # edge bleed candidates around the saturated core.
452 # We flip the cutout of the top amplifier
453 # to then work with the same coordinates for both.
454 # The way of selecting top vs bottom amp
455 # is very specific to ITL.
456 if ampName[:2] == 'C1':
457 sliceImage = maskedImage.image.array[:200, :]
458 sliceMask = maskedImage.mask.array[:200, :]
459 elif ampName[:2] == 'C0':
460 sliceImage = numpy.flipud(maskedImage.image.array[-200:, :])
461 sliceMask = numpy.flipud(maskedImage.mask.array[-200:, :])
462
463 # The middle columns of edge bleeds often have
464 # high counts, so we check there is an edge bleed
465 # by looking at a small image up to 50 pixels from the edge
466 # and around the saturated columns
467 # of the saturated core, and checking its median is
468 # above the sky background by itlEdgeBleedThreshold
469
470 # If the centroid is too close to the edge of the detector
471 # (within 5 pixels), we set the limit to the mean check
472 # to the edge of the detector
473 lowerRangeSmall = int(xCore)-5
474 upperRangeSmall = int(xCore)+5
475 if lowerRangeSmall < 0:
476 lowerRangeSmall = 0
477 if upperRangeSmall > xmax:
478 upperRangeSmall = xmax
479 ampImageBG = numpy.median(maskedImage[amp.getBBox()].image.array)
480 edgeMedian = numpy.median(sliceImage[:50, lowerRangeSmall:upperRangeSmall])
481 if edgeMedian > (ampImageBG + itlEdgeBleedThreshold):
482
483 log.info("Found ITL edge bleed in amp %s, column %d.", ampName, xCore)
484
485 # We need an estimate of the maximum width
486 # of the edge bleed for our masking model
487 # so we now estimate it by measuring the width of
488 # areas above 60 percent of the saturation level
489 # close to the edge,
490 # in a cutout up to 100 pixels from the edge,
491 # with a width of around the width of an amplifier.
492 subImageXMin = int(xCore)-250
493 subImageXMax = int(xCore)+250
494 if subImageXMin < 0:
495 subImageXMin = 0
496 elif subImageXMax > xmax:
497 subImageXMax = xmax
498
499 subImage = sliceImage[:100, subImageXMin:subImageXMax]
500 maxWidthEdgeBleed = numpy.max(numpy.sum(subImage > 0.45*satLevel,
501 axis=1))
502
503 # Mask edge bleed with a decaying exponential model
504 for y in range(200):
505 edgeBleedHalfWidth = \
506 int(((maxWidthEdgeBleed)*numpy.exp(-itlEdgeBleedModelConstant*y)
507 + widthSat)/2.)
508 lowerRange = int(xCore)-edgeBleedHalfWidth
509 upperRange = int(xCore)+edgeBleedHalfWidth
510 # If the edge bleed model goes outside the detector
511 # we set the limit for the masking
512 # to the edge of the detector
513 if lowerRange < 0:
514 lowerRange = 0
515 if upperRange > xmax:
516 upperRange = xmax
517 sliceMask[y, lowerRange:upperRange] |= saturatedBit
518
519

◆ applyGains()

lsst.ip.isr.isrFunctions.applyGains ( exposure,
normalizeGains = False,
ptcGains = None,
isTrimmed = True )
Scale an exposure by the amplifier gains.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to process.  The image is modified.
normalizeGains : `Bool`, optional
    If True, then amplifiers are scaled to force the median of
    each amplifier to equal the median of those medians.
ptcGains : `dict`[`str`], optional
    Dictionary keyed by amp name containing the PTC gains.
isTrimmed : `bool`, optional
    Is the input image trimmed?

Definition at line 1418 of file isrFunctions.py.

1418def applyGains(exposure, normalizeGains=False, ptcGains=None, isTrimmed=True):
1419 """Scale an exposure by the amplifier gains.
1420
1421 Parameters
1422 ----------
1423 exposure : `lsst.afw.image.Exposure`
1424 Exposure to process. The image is modified.
1425 normalizeGains : `Bool`, optional
1426 If True, then amplifiers are scaled to force the median of
1427 each amplifier to equal the median of those medians.
1428 ptcGains : `dict`[`str`], optional
1429 Dictionary keyed by amp name containing the PTC gains.
1430 isTrimmed : `bool`, optional
1431 Is the input image trimmed?
1432 """
1433 ccd = exposure.getDetector()
1434 ccdImage = exposure.getMaskedImage()
1435
1436 medians = []
1437 for amp in ccd:
1438 if isTrimmed:
1439 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1440 else:
1441 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1442 if ptcGains:
1443 sim *= ptcGains[amp.getName()]
1444 else:
1445 sim *= amp.getGain()
1446
1447 if normalizeGains:
1448 medians.append(numpy.median(sim.getImage().getArray()))
1449
1450 if normalizeGains:
1451 median = numpy.median(numpy.array(medians))
1452 for index, amp in enumerate(ccd):
1453 if isTrimmed:
1454 sim = ccdImage.Factory(ccdImage, amp.getBBox())
1455 else:
1456 sim = ccdImage.Factory(ccdImage, amp.getRawBBox())
1457 if medians[index] != 0.0:
1458 sim *= median/medians[index]
1459
1460

◆ attachTransmissionCurve()

lsst.ip.isr.isrFunctions.attachTransmissionCurve ( exposure,
opticsTransmission = None,
filterTransmission = None,
sensorTransmission = None,
atmosphereTransmission = None )
Attach a TransmissionCurve to an Exposure, given separate curves for
different components.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure object to modify by attaching the product of all given
    ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
    Must have a valid ``Detector`` attached that matches the detector
    associated with sensorTransmission.
opticsTransmission : `lsst.afw.image.TransmissionCurve`
    A ``TransmissionCurve`` that represents the throughput of the optics,
    to be evaluated in focal-plane coordinates.
filterTransmission : `lsst.afw.image.TransmissionCurve`
    A ``TransmissionCurve`` that represents the throughput of the filter
    itself, to be evaluated in focal-plane coordinates.
sensorTransmission : `lsst.afw.image.TransmissionCurve`
    A ``TransmissionCurve`` that represents the throughput of the sensor
    itself, to be evaluated in post-assembly trimmed detector coordinates.
atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
    A ``TransmissionCurve`` that represents the throughput of the
    atmosphere, assumed to be spatially constant.

Returns
-------
combined : `lsst.afw.image.TransmissionCurve`
    The TransmissionCurve attached to the exposure.

Notes
-----
All ``TransmissionCurve`` arguments are optional; if none are provided, the
attached ``TransmissionCurve`` will have unit transmission everywhere.

Definition at line 1366 of file isrFunctions.py.

1367 sensorTransmission=None, atmosphereTransmission=None):
1368 """Attach a TransmissionCurve to an Exposure, given separate curves for
1369 different components.
1370
1371 Parameters
1372 ----------
1373 exposure : `lsst.afw.image.Exposure`
1374 Exposure object to modify by attaching the product of all given
1375 ``TransmissionCurves`` in post-assembly trimmed detector coordinates.
1376 Must have a valid ``Detector`` attached that matches the detector
1377 associated with sensorTransmission.
1378 opticsTransmission : `lsst.afw.image.TransmissionCurve`
1379 A ``TransmissionCurve`` that represents the throughput of the optics,
1380 to be evaluated in focal-plane coordinates.
1381 filterTransmission : `lsst.afw.image.TransmissionCurve`
1382 A ``TransmissionCurve`` that represents the throughput of the filter
1383 itself, to be evaluated in focal-plane coordinates.
1384 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1385 A ``TransmissionCurve`` that represents the throughput of the sensor
1386 itself, to be evaluated in post-assembly trimmed detector coordinates.
1387 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1388 A ``TransmissionCurve`` that represents the throughput of the
1389 atmosphere, assumed to be spatially constant.
1390
1391 Returns
1392 -------
1393 combined : `lsst.afw.image.TransmissionCurve`
1394 The TransmissionCurve attached to the exposure.
1395
1396 Notes
1397 -----
1398 All ``TransmissionCurve`` arguments are optional; if none are provided, the
1399 attached ``TransmissionCurve`` will have unit transmission everywhere.
1400 """
1401 combined = afwImage.TransmissionCurve.makeIdentity()
1402 if atmosphereTransmission is not None:
1403 combined *= atmosphereTransmission
1404 if opticsTransmission is not None:
1405 combined *= opticsTransmission
1406 if filterTransmission is not None:
1407 combined *= filterTransmission
1408 detector = exposure.getDetector()
1409 fpToPix = detector.getTransform(fromSys=camGeom.FOCAL_PLANE,
1410 toSys=camGeom.PIXELS)
1411 combined = combined.transformedBy(fpToPix)
1412 if sensorTransmission is not None:
1413 combined *= sensorTransmission
1414 exposure.getInfo().setTransmissionCurve(combined)
1415 return combined
1416
1417

◆ biasCorrection()

lsst.ip.isr.isrFunctions.biasCorrection ( maskedImage,
biasMaskedImage,
trimToFit = False )
Apply bias correction in place.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
   Image to process.  The image is modified by this method.
biasMaskedImage : `lsst.afw.image.MaskedImage`
    Bias image of the same size as ``maskedImage``
trimToFit : `Bool`, optional
    If True, raw data is symmetrically trimmed to match
    calibration size.

Raises
------
RuntimeError
    Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
    the same size.

Definition at line 761 of file isrFunctions.py.

761def biasCorrection(maskedImage, biasMaskedImage, trimToFit=False):
762 """Apply bias correction in place.
763
764 Parameters
765 ----------
766 maskedImage : `lsst.afw.image.MaskedImage`
767 Image to process. The image is modified by this method.
768 biasMaskedImage : `lsst.afw.image.MaskedImage`
769 Bias image of the same size as ``maskedImage``
770 trimToFit : `Bool`, optional
771 If True, raw data is symmetrically trimmed to match
772 calibration size.
773
774 Raises
775 ------
776 RuntimeError
777 Raised if ``maskedImage`` and ``biasMaskedImage`` do not have
778 the same size.
779
780 """
781 if trimToFit:
782 maskedImage = trimToMatchCalibBBox(maskedImage, biasMaskedImage)
783
784 if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
785 raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" %
786 (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
787 maskedImage -= biasMaskedImage
788
789

◆ brighterFatterCorrection()

lsst.ip.isr.isrFunctions.brighterFatterCorrection ( exposure,
kernel,
maxIter,
threshold,
applyGain,
gains = None )
Apply brighter fatter correction in place for the image.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to have brighter-fatter correction applied.  Modified
    by this method.
kernel : `numpy.ndarray`
    Brighter-fatter kernel to apply.
maxIter : scalar
    Number of correction iterations to run.
threshold : scalar
    Convergence threshold in terms of the sum of absolute
    deviations between an iteration and the previous one.
applyGain : `Bool`
    If True, then the exposure values are scaled by the gain prior
    to correction.
gains : `dict` [`str`, `float`]
    A dictionary, keyed by amplifier name, of the gains to use.
    If gains is None, the nominal gains in the amplifier object are used.

Returns
-------
diff : `float`
    Final difference between iterations achieved in correction.
iteration : `int`
    Number of iterations used to calculate correction.

Notes
-----
This correction takes a kernel that has been derived from flat
field images to redistribute the charge.  The gradient of the
kernel is the deflection field due to the accumulated charge.

Given the original image I(x) and the kernel K(x) we can compute
the corrected image Ic(x) using the following equation:

Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))

To evaluate the derivative term we expand it as follows:

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))) )

Because we use the measured counts instead of the incident counts
we apply the correction iteratively to reconstruct the original
counts and the correction.  We stop iterating when the summed
difference between the current corrected image and the one from
the previous iteration is below the threshold.  We do not require
convergence because the number of iterations is too large a
computational cost.  How we define the threshold still needs to be
evaluated, the current default was shown to work reasonably well
on a small set of images.  For more information on the method see
DocuShare Document-19407.

The edges as defined by the kernel are not corrected because they
have spurious values due to the convolution.

Definition at line 944 of file isrFunctions.py.

944def brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None):
945 """Apply brighter fatter correction in place for the image.
946
947 Parameters
948 ----------
949 exposure : `lsst.afw.image.Exposure`
950 Exposure to have brighter-fatter correction applied. Modified
951 by this method.
952 kernel : `numpy.ndarray`
953 Brighter-fatter kernel to apply.
954 maxIter : scalar
955 Number of correction iterations to run.
956 threshold : scalar
957 Convergence threshold in terms of the sum of absolute
958 deviations between an iteration and the previous one.
959 applyGain : `Bool`
960 If True, then the exposure values are scaled by the gain prior
961 to correction.
962 gains : `dict` [`str`, `float`]
963 A dictionary, keyed by amplifier name, of the gains to use.
964 If gains is None, the nominal gains in the amplifier object are used.
965
966 Returns
967 -------
968 diff : `float`
969 Final difference between iterations achieved in correction.
970 iteration : `int`
971 Number of iterations used to calculate correction.
972
973 Notes
974 -----
975 This correction takes a kernel that has been derived from flat
976 field images to redistribute the charge. The gradient of the
977 kernel is the deflection field due to the accumulated charge.
978
979 Given the original image I(x) and the kernel K(x) we can compute
980 the corrected image Ic(x) using the following equation:
981
982 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
983
984 To evaluate the derivative term we expand it as follows:
985
986 0.5 * ( d/dx(I(x))*d/dx(int(dy*K(x-y)*I(y)))
987 + I(x)*d^2/dx^2(int(dy* K(x-y)*I(y))) )
988
989 Because we use the measured counts instead of the incident counts
990 we apply the correction iteratively to reconstruct the original
991 counts and the correction. We stop iterating when the summed
992 difference between the current corrected image and the one from
993 the previous iteration is below the threshold. We do not require
994 convergence because the number of iterations is too large a
995 computational cost. How we define the threshold still needs to be
996 evaluated, the current default was shown to work reasonably well
997 on a small set of images. For more information on the method see
998 DocuShare Document-19407.
999
1000 The edges as defined by the kernel are not corrected because they
1001 have spurious values due to the convolution.
1002 """
1003 image = exposure.getMaskedImage().getImage()
1004
1005 # The image needs to be units of electrons/holes
1006 with gainContext(exposure, image, applyGain, gains):
1007
1008 kLx = numpy.shape(kernel)[0]
1009 kLy = numpy.shape(kernel)[1]
1010 kernelImage = afwImage.ImageD(kLx, kLy)
1011 kernelImage.getArray()[:, :] = kernel
1012 tempImage = afwImage.ImageD(image, deep=True)
1013
1014 nanIndex = numpy.isnan(tempImage.getArray())
1015 tempImage.getArray()[nanIndex] = 0.
1016
1017 corr = numpy.zeros(image.array.shape, dtype=numpy.float64)
1018 prev_image = numpy.zeros(image.array.shape, dtype=numpy.float64)
1019
1020 # Define boundary by convolution region. The region that the
1021 # correction will be calculated for is one fewer in each dimension
1022 # because of the second derivative terms.
1023 # NOTE: these need to use integer math, as we're using start:end as
1024 # numpy index ranges.
1025 startX = kLx//2
1026 endX = -kLx//2
1027 startY = kLy//2
1028 endY = -kLy//2
1029
1030 for iteration in range(maxIter):
1031
1032 outArray = scipy.signal.convolve(
1033 tempImage.array,
1034 kernelImage.array,
1035 mode="same",
1036 method="fft",
1037 )
1038 tmpArray = tempImage.getArray()
1039
1040 with numpy.errstate(invalid="ignore", over="ignore"):
1041 # First derivative term
1042 gradTmp = numpy.gradient(tmpArray[startY:endY, startX:endX])
1043 gradOut = numpy.gradient(outArray[startY:endY, startX:endX])
1044 first = (gradTmp[0]*gradOut[0] + gradTmp[1]*gradOut[1])[1:-1, 1:-1]
1045
1046 # Second derivative term
1047 diffOut20 = numpy.diff(outArray, 2, 0)[startY:endY, startX + 1:endX - 1]
1048 diffOut21 = numpy.diff(outArray, 2, 1)[startY + 1:endY - 1, startX:endX]
1049 second = tmpArray[startY + 1:endY - 1, startX + 1:endX - 1]*(diffOut20 + diffOut21)
1050
1051 corr[startY + 1:endY - 1, startX + 1:endX - 1] = 0.5*(first + second)
1052
1053 tmpArray[:, :] = image.getArray()[:, :]
1054 tmpArray[nanIndex] = 0.
1055 tmpArray[startY:endY, startX:endX] += corr[startY:endY, startX:endX]
1056
1057 if iteration > 0:
1058 diff = numpy.sum(numpy.abs(prev_image - tmpArray), dtype=numpy.float64)
1059
1060 if diff < threshold:
1061 break
1062 prev_image[:, :] = tmpArray[:, :]
1063
1064 image.getArray()[startY + 1:endY - 1, startX + 1:endX - 1] += \
1065 corr[startY + 1:endY - 1, startX + 1:endX - 1]
1066
1067 return diff, iteration
1068
1069

◆ checkFilter()

lsst.ip.isr.isrFunctions.checkFilter ( exposure,
filterList,
log )
Check to see if an exposure is in a filter specified by a list.

The goal of this is to provide a unified filter checking interface
for all filter dependent stages.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to examine.
filterList : `list` [`str`]
    List of physical_filter names to check.
log : `logging.Logger`
    Logger to handle messages.

Returns
-------
result : `bool`
    True if the exposure's filter is contained in the list.

Definition at line 1556 of file isrFunctions.py.

1556def checkFilter(exposure, filterList, log):
1557 """Check to see if an exposure is in a filter specified by a list.
1558
1559 The goal of this is to provide a unified filter checking interface
1560 for all filter dependent stages.
1561
1562 Parameters
1563 ----------
1564 exposure : `lsst.afw.image.Exposure`
1565 Exposure to examine.
1566 filterList : `list` [`str`]
1567 List of physical_filter names to check.
1568 log : `logging.Logger`
1569 Logger to handle messages.
1570
1571 Returns
1572 -------
1573 result : `bool`
1574 True if the exposure's filter is contained in the list.
1575 """
1576 if len(filterList) == 0:
1577 return False
1578 thisFilter = exposure.getFilter()
1579 if thisFilter is None:
1580 log.warning("No FilterLabel attached to this exposure!")
1581 return False
1582
1583 thisPhysicalFilter = getPhysicalFilter(thisFilter, log)
1584 if thisPhysicalFilter in filterList:
1585 return True
1586 elif thisFilter.bandLabel in filterList:
1587 if log:
1588 log.warning("Physical filter (%s) should be used instead of band %s for filter configurations"
1589 " (%s)", thisPhysicalFilter, thisFilter.bandLabel, filterList)
1590 return True
1591 else:
1592 return False
1593
1594

◆ compareCameraKeywords()

lsst.ip.isr.isrFunctions.compareCameraKeywords ( doRaiseOnCalibMismatch,
cameraKeywordsToCompare,
exposureMetadata,
calib,
calibName,
log = None )
Compare header keywords to confirm camera states match.

Parameters
----------
doRaiseOnCalibMismatch : `bool`
    Raise on calibration mismatch?  Otherwise, log a warning.
cameraKeywordsToCompare : `list` [`str`]
    List of camera keywords to compare.
exposureMetadata : `lsst.daf.base.PropertyList`
    Header for the exposure being processed.
calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
    Calibration to be applied.
calibName : `str`
    Calib type for log message.
log : `logging.Logger`, optional
    Logger to handle messages.

Definition at line 1748 of file isrFunctions.py.

1755):
1756 """Compare header keywords to confirm camera states match.
1757
1758 Parameters
1759 ----------
1760 doRaiseOnCalibMismatch : `bool`
1761 Raise on calibration mismatch? Otherwise, log a warning.
1762 cameraKeywordsToCompare : `list` [`str`]
1763 List of camera keywords to compare.
1764 exposureMetadata : `lsst.daf.base.PropertyList`
1765 Header for the exposure being processed.
1766 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1767 Calibration to be applied.
1768 calibName : `str`
1769 Calib type for log message.
1770 log : `logging.Logger`, optional
1771 Logger to handle messages.
1772 """
1773 try:
1774 calibMetadata = calib.metadata
1775 except AttributeError:
1776 return
1777
1778 log = log if log else logging.getLogger(__name__)
1779
1780 missingKeywords = []
1781 for keyword in cameraKeywordsToCompare:
1782 exposureValue = exposureMetadata.get(keyword, None)
1783 if exposureValue is None:
1784 log.debug("Sequencer keyword %s not found in exposure metadata.", keyword)
1785 continue
1786
1787 calibValue = calibMetadata.get(keyword, None)
1788
1789 # We don't log here if there is a missing keyword.
1790 if calibValue is None:
1791 missingKeywords.append(keyword)
1792 continue
1793
1794 if exposureValue != calibValue:
1795 if doRaiseOnCalibMismatch:
1796 raise RuntimeError(
1797 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1798 calibName,
1799 keyword,
1800 exposureValue,
1801 calibValue,
1802 )
1803 else:
1804 log.warning(
1805 "Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
1806 calibName,
1807 keyword,
1808 exposureValue,
1809 calibValue,
1810 )
1811 exposureMetadata[f"ISR {calibName.upper()} SEQUENCER MISMATCH"] = True
1812
1813 if missingKeywords:
1814 log.info(
1815 "Calibration %s missing keywords %s, which were not checked.",
1816 calibName,
1817 ",".join(missingKeywords),
1818 )

◆ countMaskedPixels()

lsst.ip.isr.isrFunctions.countMaskedPixels ( maskedIm,
maskPlane )
Count the number of pixels in a given mask plane.

Parameters
----------
maskedIm : `~lsst.afw.image.MaskedImage`
    Masked image to examine.
maskPlane : `str`
    Name of the mask plane to examine.

Returns
-------
nPix : `int`
    Number of pixels in the requested mask plane.

Definition at line 1628 of file isrFunctions.py.

1628def countMaskedPixels(maskedIm, maskPlane):
1629 """Count the number of pixels in a given mask plane.
1630
1631 Parameters
1632 ----------
1633 maskedIm : `~lsst.afw.image.MaskedImage`
1634 Masked image to examine.
1635 maskPlane : `str`
1636 Name of the mask plane to examine.
1637
1638 Returns
1639 -------
1640 nPix : `int`
1641 Number of pixels in the requested mask plane.
1642 """
1643 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)
1644 nPix = numpy.where(numpy.bitwise_and(maskedIm.mask.array, maskBit))[0].flatten().size
1645 return nPix
1646
1647

◆ createPsf()

lsst.ip.isr.isrFunctions.createPsf ( fwhm)
Make a double Gaussian PSF.

Parameters
----------
fwhm : scalar
    FWHM of double Gaussian smoothing kernel.

Returns
-------
psf : `lsst.meas.algorithms.DoubleGaussianPsf`
    The created smoothing kernel.

Definition at line 77 of file isrFunctions.py.

77def createPsf(fwhm):
78 """Make a double Gaussian PSF.
79
80 Parameters
81 ----------
82 fwhm : scalar
83 FWHM of double Gaussian smoothing kernel.
84
85 Returns
86 -------
87 psf : `lsst.meas.algorithms.DoubleGaussianPsf`
88 The created smoothing kernel.
89 """
90 ksize = 4*int(fwhm) + 1
91 return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
92
93

◆ darkCorrection()

lsst.ip.isr.isrFunctions.darkCorrection ( maskedImage,
darkMaskedImage,
expScale,
darkScale,
invert = False,
trimToFit = False )
Apply dark correction in place.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
   Image to process.  The image is modified by this method.
darkMaskedImage : `lsst.afw.image.MaskedImage`
    Dark image of the same size as ``maskedImage``.
expScale : scalar
    Dark exposure time for ``maskedImage``.
darkScale : scalar
    Dark exposure time for ``darkMaskedImage``.
invert : `Bool`, optional
    If True, re-add the dark to an already corrected image.
trimToFit : `Bool`, optional
    If True, raw data is symmetrically trimmed to match
    calibration size.

Raises
------
RuntimeError
    Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
    the same size.

Notes
-----
The dark correction is applied by calculating:
    maskedImage -= dark * expScaling / darkScaling

Definition at line 790 of file isrFunctions.py.

790def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale, invert=False, trimToFit=False):
791 """Apply dark correction in place.
792
793 Parameters
794 ----------
795 maskedImage : `lsst.afw.image.MaskedImage`
796 Image to process. The image is modified by this method.
797 darkMaskedImage : `lsst.afw.image.MaskedImage`
798 Dark image of the same size as ``maskedImage``.
799 expScale : scalar
800 Dark exposure time for ``maskedImage``.
801 darkScale : scalar
802 Dark exposure time for ``darkMaskedImage``.
803 invert : `Bool`, optional
804 If True, re-add the dark to an already corrected image.
805 trimToFit : `Bool`, optional
806 If True, raw data is symmetrically trimmed to match
807 calibration size.
808
809 Raises
810 ------
811 RuntimeError
812 Raised if ``maskedImage`` and ``darkMaskedImage`` do not have
813 the same size.
814
815 Notes
816 -----
817 The dark correction is applied by calculating:
818 maskedImage -= dark * expScaling / darkScaling
819 """
820 if trimToFit:
821 maskedImage = trimToMatchCalibBBox(maskedImage, darkMaskedImage)
822
823 if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
824 raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" %
825 (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
826
827 scale = expScale / darkScale
828 if not invert:
829 maskedImage.scaledMinus(scale, darkMaskedImage)
830 else:
831 maskedImage.scaledPlus(scale, darkMaskedImage)
832
833

◆ flatCorrection()

lsst.ip.isr.isrFunctions.flatCorrection ( maskedImage,
flatMaskedImage,
scalingType,
userScale = 1.0,
invert = False,
trimToFit = False )
Apply flat correction in place.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.  The image is modified.
flatMaskedImage : `lsst.afw.image.MaskedImage`
    Flat image of the same size as ``maskedImage``
scalingType : str
    Flat scale computation method.  Allowed values are 'MEAN',
    'MEDIAN', or 'USER'.
userScale : scalar, optional
    Scale to use if ``scalingType='USER'``.
invert : `Bool`, optional
    If True, unflatten an already flattened image.
trimToFit : `Bool`, optional
    If True, raw data is symmetrically trimmed to match
    calibration size.

Raises
------
RuntimeError
    Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
    the same size or if ``scalingType`` is not an allowed value.

Definition at line 862 of file isrFunctions.py.

862def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0, invert=False, trimToFit=False):
863 """Apply flat correction in place.
864
865 Parameters
866 ----------
867 maskedImage : `lsst.afw.image.MaskedImage`
868 Image to process. The image is modified.
869 flatMaskedImage : `lsst.afw.image.MaskedImage`
870 Flat image of the same size as ``maskedImage``
871 scalingType : str
872 Flat scale computation method. Allowed values are 'MEAN',
873 'MEDIAN', or 'USER'.
874 userScale : scalar, optional
875 Scale to use if ``scalingType='USER'``.
876 invert : `Bool`, optional
877 If True, unflatten an already flattened image.
878 trimToFit : `Bool`, optional
879 If True, raw data is symmetrically trimmed to match
880 calibration size.
881
882 Raises
883 ------
884 RuntimeError
885 Raised if ``maskedImage`` and ``flatMaskedImage`` do not have
886 the same size or if ``scalingType`` is not an allowed value.
887 """
888 if trimToFit:
889 maskedImage = trimToMatchCalibBBox(maskedImage, flatMaskedImage)
890
891 if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
892 raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" %
893 (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
894
895 # Figure out scale from the data
896 # Ideally the flats are normalized by the calibration product pipeline,
897 # but this allows some flexibility in the case that the flat is created by
898 # some other mechanism.
899 if scalingType in ('MEAN', 'MEDIAN'):
900 scalingType = afwMath.stringToStatisticsProperty(scalingType)
901 flatScale = afwMath.makeStatistics(flatMaskedImage.image, scalingType).getValue()
902 elif scalingType == 'USER':
903 flatScale = userScale
904 else:
905 raise RuntimeError('%s : %s not implemented' % ("flatCorrection", scalingType))
906
907 if not invert:
908 maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
909 else:
910 maskedImage.scaledMultiplies(1.0/flatScale, flatMaskedImage)
911
912
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
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)

◆ fluxConservingBrighterFatterCorrection()

lsst.ip.isr.isrFunctions.fluxConservingBrighterFatterCorrection ( exposure,
kernel,
maxIter,
threshold,
applyGain,
gains = None,
correctionMode = True )
Apply brighter fatter correction in place for the image.

This version presents a modified version of the algorithm
found in ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``
which conserves the image flux, resulting in improved
correction of the cores of stars. The convolution has also been
modified to mitigate edge effects.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to have brighter-fatter correction applied.  Modified
    by this method.
kernel : `numpy.ndarray`
    Brighter-fatter kernel to apply.
maxIter : scalar
    Number of correction iterations to run.
threshold : scalar
    Convergence threshold in terms of the sum of absolute
    deviations between an iteration and the previous one.
applyGain : `Bool`
    If True, then the exposure values are scaled by the gain prior
    to correction.
gains : `dict` [`str`, `float`]
    A dictionary, keyed by amplifier name, of the gains to use.
    If gains is None, the nominal gains in the amplifier object are used.
correctionMode : `Bool`
    If True (default) the function applies correction for BFE.  If False,
    the code can instead be used to generate a simulation of BFE (sign
    change in the direction of the effect)

Returns
-------
diff : `float`
    Final difference between iterations achieved in correction.
iteration : `int`
    Number of iterations used to calculate correction.

Notes
-----
Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.

This correction takes a kernel that has been derived from flat
field images to redistribute the charge.  The gradient of the
kernel is the deflection field due to the accumulated charge.

Given the original image I(x) and the kernel K(x) we can compute
the corrected image Ic(x) using the following equation:

Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))

Improved algorithm at this step applies the divergence theorem to
obtain a pixelised correction.

Because we use the measured counts instead of the incident counts
we apply the correction iteratively to reconstruct the original
counts and the correction.  We stop iterating when the summed
difference between the current corrected image and the one from
the previous iteration is below the threshold.  We do not require
convergence because the number of iterations is too large a
computational cost.  How we define the threshold still needs to be
evaluated, the current default was shown to work reasonably well
on a small set of images.

Edges are handled in the convolution by padding.  This is still not
a physical model for the edge, but avoids discontinuity in the correction.

Author of modified version: Lance.Miller@physics.ox.ac.uk
(see DM-38555).

Definition at line 1154 of file isrFunctions.py.

1155 gains=None, correctionMode=True):
1156 """Apply brighter fatter correction in place for the image.
1157
1158 This version presents a modified version of the algorithm
1159 found in ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``
1160 which conserves the image flux, resulting in improved
1161 correction of the cores of stars. The convolution has also been
1162 modified to mitigate edge effects.
1163
1164 Parameters
1165 ----------
1166 exposure : `lsst.afw.image.Exposure`
1167 Exposure to have brighter-fatter correction applied. Modified
1168 by this method.
1169 kernel : `numpy.ndarray`
1170 Brighter-fatter kernel to apply.
1171 maxIter : scalar
1172 Number of correction iterations to run.
1173 threshold : scalar
1174 Convergence threshold in terms of the sum of absolute
1175 deviations between an iteration and the previous one.
1176 applyGain : `Bool`
1177 If True, then the exposure values are scaled by the gain prior
1178 to correction.
1179 gains : `dict` [`str`, `float`]
1180 A dictionary, keyed by amplifier name, of the gains to use.
1181 If gains is None, the nominal gains in the amplifier object are used.
1182 correctionMode : `Bool`
1183 If True (default) the function applies correction for BFE. If False,
1184 the code can instead be used to generate a simulation of BFE (sign
1185 change in the direction of the effect)
1186
1187 Returns
1188 -------
1189 diff : `float`
1190 Final difference between iterations achieved in correction.
1191 iteration : `int`
1192 Number of iterations used to calculate correction.
1193
1194 Notes
1195 -----
1196 Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.
1197
1198 This correction takes a kernel that has been derived from flat
1199 field images to redistribute the charge. The gradient of the
1200 kernel is the deflection field due to the accumulated charge.
1201
1202 Given the original image I(x) and the kernel K(x) we can compute
1203 the corrected image Ic(x) using the following equation:
1204
1205 Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
1206
1207 Improved algorithm at this step applies the divergence theorem to
1208 obtain a pixelised correction.
1209
1210 Because we use the measured counts instead of the incident counts
1211 we apply the correction iteratively to reconstruct the original
1212 counts and the correction. We stop iterating when the summed
1213 difference between the current corrected image and the one from
1214 the previous iteration is below the threshold. We do not require
1215 convergence because the number of iterations is too large a
1216 computational cost. How we define the threshold still needs to be
1217 evaluated, the current default was shown to work reasonably well
1218 on a small set of images.
1219
1220 Edges are handled in the convolution by padding. This is still not
1221 a physical model for the edge, but avoids discontinuity in the correction.
1222
1223 Author of modified version: Lance.Miller@physics.ox.ac.uk
1224 (see DM-38555).
1225 """
1226 image = exposure.getMaskedImage().getImage()
1227
1228 # The image needs to be units of electrons/holes
1229 with gainContext(exposure, image, applyGain, gains):
1230
1231 # get kernel and its shape
1232 kLy, kLx = kernel.shape
1233 kernelImage = afwImage.ImageD(kLx, kLy)
1234 kernelImage.getArray()[:, :] = kernel
1235 tempImage = afwImage.ImageD(image, deep=True)
1236
1237 nanIndex = numpy.isnan(tempImage.getArray())
1238 tempImage.getArray()[nanIndex] = 0.
1239
1240 outImage = afwImage.ImageD(image.getDimensions())
1241 corr = numpy.zeros(image.array.shape, dtype=numpy.float64)
1242 prevImage = numpy.zeros(image.array.shape, dtype=numpy.float64)
1243 convCntrl = afwMath.ConvolutionControl(False, False, 1)
1244 fixedKernel = afwMath.FixedKernel(kernelImage)
1245
1246 # set the padding amount
1247 # ensure we pad by an even amount larger than the kernel
1248 kLy = 2 * ((1+kLy)//2)
1249 kLx = 2 * ((1+kLx)//2)
1250
1251 # The deflection potential only depends on the gradient of
1252 # the convolution, so we can subtract the mean, which then
1253 # allows us to pad the image with zeros and avoid wrap-around effects
1254 # (although still not handling the image edges with a physical model)
1255 # This wouldn't be great if there were a strong image gradient.
1256 imYdimension, imXdimension = tempImage.array.shape
1257 imean = numpy.mean(tempImage.getArray()[~nanIndex], dtype=numpy.float64)
1258 # subtract mean from image
1259 tempImage -= imean
1260 tempImage.array[nanIndex] = 0.0
1261 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
1262 outImage = afwImage.ImageD(numpy.pad(outImage.getArray(), ((0, kLy), (0, kLx))))
1263 # Convert array to afw image so afwMath.convolve works
1264 padImage = afwImage.ImageD(padArray.shape[1], padArray.shape[0])
1265 padImage.array[:] = padArray
1266
1267 for iteration in range(maxIter):
1268
1269 # create deflection potential, convolution of flux with kernel
1270 # using padded counts array
1271 afwMath.convolve(outImage, padImage, fixedKernel, convCntrl)
1272 tmpArray = tempImage.getArray()
1273 outArray = outImage.getArray()
1274
1275 # trim convolution output back to original shape
1276 outArray = outArray[:imYdimension, :imXdimension]
1277
1278 # generate the correction array, with correctionMode set as input
1279 corr[...] = transferFlux(outArray, tmpArray, correctionMode=correctionMode)
1280
1281 # update the arrays for the next iteration
1282 tmpArray[:, :] = image.getArray()[:, :]
1283 tmpArray += corr
1284 tmpArray[nanIndex] = 0.
1285 # update padded array
1286 # subtract mean
1287 tmpArray -= imean
1288 tempImage.array[nanIndex] = 0.
1289 padArray = numpy.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
1290
1291 if iteration > 0:
1292 diff = numpy.sum(numpy.abs(prevImage - tmpArray), dtype=numpy.float64)
1293
1294 if diff < threshold:
1295 break
1296 prevImage[:, :] = tmpArray[:, :]
1297
1298 image.getArray()[:] += corr[:]
1299
1300 return diff, iteration
1301
1302
1303@contextmanager
Parameters to control convolution.
A kernel created from an Image.
Definition Kernel.h:471
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.

◆ gainContext()

lsst.ip.isr.isrFunctions.gainContext ( exp,
image,
apply,
gains = None,
invert = False,
isTrimmed = True )
Context manager that applies and removes gain.

Parameters
----------
exp : `lsst.afw.image.Exposure`
    Exposure to apply/remove gain.
image : `lsst.afw.image.Image`
    Image to apply/remove gain.
apply : `bool`
    If True, apply and remove the amplifier gain.
gains : `dict` [`str`, `float`], optional
    A dictionary, keyed by amplifier name, of the gains to use.
    If gains is None, the nominal gains in the amplifier object are used.
invert : `bool`, optional
    Invert the gains (e.g. convert electrons to adu temporarily)?
isTrimmed : `bool`, optional
    Is this a trimmed exposure?

Yields
------
exp : `lsst.afw.image.Exposure`
    Exposure with the gain applied.

Definition at line 1304 of file isrFunctions.py.

1304def gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True):
1305 """Context manager that applies and removes gain.
1306
1307 Parameters
1308 ----------
1309 exp : `lsst.afw.image.Exposure`
1310 Exposure to apply/remove gain.
1311 image : `lsst.afw.image.Image`
1312 Image to apply/remove gain.
1313 apply : `bool`
1314 If True, apply and remove the amplifier gain.
1315 gains : `dict` [`str`, `float`], optional
1316 A dictionary, keyed by amplifier name, of the gains to use.
1317 If gains is None, the nominal gains in the amplifier object are used.
1318 invert : `bool`, optional
1319 Invert the gains (e.g. convert electrons to adu temporarily)?
1320 isTrimmed : `bool`, optional
1321 Is this a trimmed exposure?
1322
1323 Yields
1324 ------
1325 exp : `lsst.afw.image.Exposure`
1326 Exposure with the gain applied.
1327 """
1328 # check we have all of them if provided because mixing and matching would
1329 # be a real mess
1330 if gains and apply is True:
1331 ampNames = [amp.getName() for amp in exp.getDetector()]
1332 for ampName in ampNames:
1333 if ampName not in gains.keys():
1334 raise RuntimeError(f"Gains provided to gain context, but no entry found for amp {ampName}")
1335
1336 if apply:
1337 ccd = exp.getDetector()
1338 for amp in ccd:
1339 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1340 if gains:
1341 gain = gains[amp.getName()]
1342 else:
1343 gain = amp.getGain()
1344 if invert:
1345 sim /= gain
1346 else:
1347 sim *= gain
1348
1349 try:
1350 yield exp
1351 finally:
1352 if apply:
1353 ccd = exp.getDetector()
1354 for amp in ccd:
1355 sim = image.Factory(image, amp.getBBox() if isTrimmed else amp.getRawBBox())
1356 if gains:
1357 gain = gains[amp.getName()]
1358 else:
1359 gain = amp.getGain()
1360 if invert:
1361 sim *= gain
1362 else:
1363 sim /= gain
1364
1365

◆ getExposureGains()

lsst.ip.isr.isrFunctions.getExposureGains ( exposure)
Get the per-amplifier gains used for this exposure.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    The exposure to find gains for.

Returns
-------
gains : `dict` [`str` `float`]
    Dictionary of gain values, keyed by amplifier name.
    Returns empty dict when detector is None.

Definition at line 1648 of file isrFunctions.py.

1648def getExposureGains(exposure):
1649 """Get the per-amplifier gains used for this exposure.
1650
1651 Parameters
1652 ----------
1653 exposure : `lsst.afw.image.Exposure`
1654 The exposure to find gains for.
1655
1656 Returns
1657 -------
1658 gains : `dict` [`str` `float`]
1659 Dictionary of gain values, keyed by amplifier name.
1660 Returns empty dict when detector is None.
1661 """
1662 det = exposure.getDetector()
1663 if det is None:
1664 return dict()
1665
1666 metadata = exposure.getMetadata()
1667 gains = {}
1668 for amp in det:
1669 ampName = amp.getName()
1670 # The key may use the new LSST ISR or the old LSST prefix
1671 if (key1 := f"LSST ISR GAIN {ampName}") in metadata:
1672 gains[ampName] = metadata[key1]
1673 elif (key2 := f"LSST GAIN {ampName}") in metadata:
1674 gains[ampName] = metadata[key2]
1675 else:
1676 gains[ampName] = amp.getGain()
1677 return gains
1678
1679

◆ getExposureReadNoises()

lsst.ip.isr.isrFunctions.getExposureReadNoises ( exposure)
Get the per-amplifier read noise used for this exposure.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    The exposure to find read noise for.

Returns
-------
readnoises : `dict` [`str` `float`]
    Dictionary of read noise values, keyed by amplifier name.
    Returns empty dict when detector is None.

Definition at line 1680 of file isrFunctions.py.

1680def getExposureReadNoises(exposure):
1681 """Get the per-amplifier read noise used for this exposure.
1682
1683 Parameters
1684 ----------
1685 exposure : `lsst.afw.image.Exposure`
1686 The exposure to find read noise for.
1687
1688 Returns
1689 -------
1690 readnoises : `dict` [`str` `float`]
1691 Dictionary of read noise values, keyed by amplifier name.
1692 Returns empty dict when detector is None.
1693 """
1694 det = exposure.getDetector()
1695 if det is None:
1696 return dict()
1697
1698 metadata = exposure.getMetadata()
1699 readnoises = {}
1700 for amp in det:
1701 ampName = amp.getName()
1702 # The key may use the new LSST ISR or the old LSST prefix
1703 if (key1 := f"LSST ISR READNOISE {ampName}") in metadata:
1704 readnoises[ampName] = metadata[key1]
1705 elif (key2 := f"LSST READNOISE {ampName}") in metadata:
1706 readnoises[ampName] = metadata[key2]
1707 else:
1708 readnoises[ampName] = amp.getReadNoise()
1709 return readnoises
1710
1711

◆ getPhysicalFilter()

lsst.ip.isr.isrFunctions.getPhysicalFilter ( filterLabel,
log )
Get the physical filter label associated with the given filterLabel.

If ``filterLabel`` is `None` or there is no physicalLabel attribute
associated with the given ``filterLabel``, the returned label will be
"Unknown".

Parameters
----------
filterLabel : `lsst.afw.image.FilterLabel`
    The `lsst.afw.image.FilterLabel` object from which to derive the
    physical filter label.
log : `logging.Logger`
    Logger to handle messages.

Returns
-------
physicalFilter : `str`
    The value returned by the physicalLabel attribute of ``filterLabel`` if
    it exists, otherwise set to \"Unknown\".

Definition at line 1595 of file isrFunctions.py.

1595def getPhysicalFilter(filterLabel, log):
1596 """Get the physical filter label associated with the given filterLabel.
1597
1598 If ``filterLabel`` is `None` or there is no physicalLabel attribute
1599 associated with the given ``filterLabel``, the returned label will be
1600 "Unknown".
1601
1602 Parameters
1603 ----------
1604 filterLabel : `lsst.afw.image.FilterLabel`
1605 The `lsst.afw.image.FilterLabel` object from which to derive the
1606 physical filter label.
1607 log : `logging.Logger`
1608 Logger to handle messages.
1609
1610 Returns
1611 -------
1612 physicalFilter : `str`
1613 The value returned by the physicalLabel attribute of ``filterLabel`` if
1614 it exists, otherwise set to \"Unknown\".
1615 """
1616 if filterLabel is None:
1617 physicalFilter = "Unknown"
1618 log.warning("filterLabel is None. Setting physicalFilter to \"Unknown\".")
1619 else:
1620 try:
1621 physicalFilter = filterLabel.physicalLabel
1622 except RuntimeError:
1623 log.warning("filterLabel has no physicalLabel attribute. Setting physicalFilter to \"Unknown\".")
1624 physicalFilter = "Unknown"
1625 return physicalFilter
1626
1627

◆ growMasks()

lsst.ip.isr.isrFunctions.growMasks ( mask,
radius = 0,
maskNameList = ['BAD'],
maskValue = "BAD" )
Grow a mask by an amount and add to the requested plane.

Parameters
----------
mask : `lsst.afw.image.Mask`
    Mask image to process.
radius : scalar
    Amount to grow the mask.
maskNameList : `str` or `list` [`str`]
    Mask names that should be grown.
maskValue : `str`
    Mask plane to assign the newly masked pixels to.

Definition at line 206 of file isrFunctions.py.

206def growMasks(mask, radius=0, maskNameList=['BAD'], maskValue="BAD"):
207 """Grow a mask by an amount and add to the requested plane.
208
209 Parameters
210 ----------
211 mask : `lsst.afw.image.Mask`
212 Mask image to process.
213 radius : scalar
214 Amount to grow the mask.
215 maskNameList : `str` or `list` [`str`]
216 Mask names that should be grown.
217 maskValue : `str`
218 Mask plane to assign the newly masked pixels to.
219 """
220 if radius > 0:
221 spans = SpanSet.fromMask(mask, mask.getPlaneBitMask(maskNameList))
222 # Use MANHATTAN for equivalence with 'isotropic=False` footprint grows,
223 # but CIRCLE is probably better and might be just as fast.
224 spans = spans.dilated(radius, Stencil.MANHATTAN)
225 spans = spans.clippedTo(mask.getBBox())
226 spans.setMask(mask, mask.getPlaneBitMask(maskValue))
227
228

◆ illuminationCorrection()

lsst.ip.isr.isrFunctions.illuminationCorrection ( maskedImage,
illumMaskedImage,
illumScale,
trimToFit = True )
Apply illumination correction in place.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.  The image is modified.
illumMaskedImage : `lsst.afw.image.MaskedImage`
    Illumination correction image of the same size as ``maskedImage``.
illumScale : scalar
    Scale factor for the illumination correction.
trimToFit : `Bool`, optional
    If True, raw data is symmetrically trimmed to match
    calibration size.

Raises
------
RuntimeError
    Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
    the same size.

Definition at line 913 of file isrFunctions.py.

913def illuminationCorrection(maskedImage, illumMaskedImage, illumScale, trimToFit=True):
914 """Apply illumination correction in place.
915
916 Parameters
917 ----------
918 maskedImage : `lsst.afw.image.MaskedImage`
919 Image to process. The image is modified.
920 illumMaskedImage : `lsst.afw.image.MaskedImage`
921 Illumination correction image of the same size as ``maskedImage``.
922 illumScale : scalar
923 Scale factor for the illumination correction.
924 trimToFit : `Bool`, optional
925 If True, raw data is symmetrically trimmed to match
926 calibration size.
927
928 Raises
929 ------
930 RuntimeError
931 Raised if ``maskedImage`` and ``illumMaskedImage`` do not have
932 the same size.
933 """
934 if trimToFit:
935 maskedImage = trimToMatchCalibBBox(maskedImage, illumMaskedImage)
936
937 if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
938 raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" %
939 (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
940
941 maskedImage.scaledDivides(1.0/illumScale, illumMaskedImage)
942
943

◆ interpolateDefectList()

lsst.ip.isr.isrFunctions.interpolateDefectList ( maskedImage,
defectList,
fwhm,
fallbackValue = None,
maskNameList = None,
useLegacyInterp = True )
Interpolate over defects specified in a defect list.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.
defectList : `lsst.meas.algorithms.Defects`
    List of defects to interpolate over.
fwhm : `float`
    FWHM of double Gaussian smoothing kernel.
fallbackValue : scalar, optional
    Fallback value if an interpolated value cannot be determined.
    If None, then the clipped mean of the image is used.
maskNameList : `list [string]`
    List of the defects to interpolate over (used for GP interpolator).
useLegacyInterp : `bool`
    Use the legacy interpolation (polynomial interpolation) if True. Use
    Gaussian Process interpolation if False.

Notes
-----
The ``fwhm`` parameter is used to create a PSF, but the underlying
interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
not currently make use of this information in legacy Interpolation, but use
if for the Gaussian Process as an estimation of the correlation lenght.

Definition at line 114 of file isrFunctions.py.

115 maskNameList=None, useLegacyInterp=True):
116 """Interpolate over defects specified in a defect list.
117
118 Parameters
119 ----------
120 maskedImage : `lsst.afw.image.MaskedImage`
121 Image to process.
122 defectList : `lsst.meas.algorithms.Defects`
123 List of defects to interpolate over.
124 fwhm : `float`
125 FWHM of double Gaussian smoothing kernel.
126 fallbackValue : scalar, optional
127 Fallback value if an interpolated value cannot be determined.
128 If None, then the clipped mean of the image is used.
129 maskNameList : `list [string]`
130 List of the defects to interpolate over (used for GP interpolator).
131 useLegacyInterp : `bool`
132 Use the legacy interpolation (polynomial interpolation) if True. Use
133 Gaussian Process interpolation if False.
134
135 Notes
136 -----
137 The ``fwhm`` parameter is used to create a PSF, but the underlying
138 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
139 not currently make use of this information in legacy Interpolation, but use
140 if for the Gaussian Process as an estimation of the correlation lenght.
141 """
142 psf = createPsf(fwhm)
143 if fallbackValue is None:
144 fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
145 if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict():
146 maskedImage.getMask().addMaskPlane('INTRP')
147
148 # Hardcoded fwhm value. PSF estimated latter in step1,
149 # not in ISR.
150 if useLegacyInterp:
151 kwargs = {}
152 fwhm = fwhm
153 else:
154 # tested on a dozens of images and looks a good set of
155 # hyperparameters, but cannot guarrenty this is optimal,
156 # need further testing.
157 kwargs = {"bin_spacing": 20,
158 "threshold_dynamic_binning": 2000,
159 "threshold_subdivide": 20000}
160 fwhm = 15
161
162 measAlg.interpolateOverDefects(maskedImage, psf, defectList,
163 fallbackValue=fallbackValue,
164 useFallbackValueAtEdge=True,
165 fwhm=fwhm,
166 useLegacyInterp=useLegacyInterp,
167 maskNameList=maskNameList, **kwargs)
168 return maskedImage
169
170

◆ interpolateFromMask()

lsst.ip.isr.isrFunctions.interpolateFromMask ( maskedImage,
fwhm,
growSaturatedFootprints = 1,
maskNameList = ['SAT'],
fallbackValue = None,
useLegacyInterp = True )
Interpolate over defects identified by a particular set of mask planes.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.
fwhm : `float`
    FWHM of double Gaussian smoothing kernel.
growSaturatedFootprints : scalar, optional
    Number of pixels to grow footprints for saturated pixels.
maskNameList : `List` of `str`, optional
    Mask plane name.
fallbackValue : scalar, optional
    Value of last resort for interpolation.

Notes
-----
The ``fwhm`` parameter is used to create a PSF, but the underlying
interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
not currently make use of this information.

Definition at line 632 of file isrFunctions.py.

633 maskNameList=['SAT'], fallbackValue=None, useLegacyInterp=True):
634 """Interpolate over defects identified by a particular set of mask planes.
635
636 Parameters
637 ----------
638 maskedImage : `lsst.afw.image.MaskedImage`
639 Image to process.
640 fwhm : `float`
641 FWHM of double Gaussian smoothing kernel.
642 growSaturatedFootprints : scalar, optional
643 Number of pixels to grow footprints for saturated pixels.
644 maskNameList : `List` of `str`, optional
645 Mask plane name.
646 fallbackValue : scalar, optional
647 Value of last resort for interpolation.
648
649 Notes
650 -----
651 The ``fwhm`` parameter is used to create a PSF, but the underlying
652 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
653 not currently make use of this information.
654 """
655 mask = maskedImage.getMask()
656
657 if growSaturatedFootprints > 0 and "SAT" in maskNameList:
658 # If we are interpolating over an area larger than the original masked
659 # region, we need to expand the original mask bit to the full area to
660 # explain why we interpolated there.
661 growMasks(mask, radius=growSaturatedFootprints, maskNameList=['SAT'], maskValue="SAT")
662
663 thresh = afwDetection.Threshold(mask.getPlaneBitMask(maskNameList), afwDetection.Threshold.BITMASK)
664 fpSet = afwDetection.FootprintSet(mask, thresh)
665 defectList = Defects.fromFootprintList(fpSet.getFootprints())
666
667 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
668 maskNameList=maskNameList, useLegacyInterp=useLegacyInterp)
669
670 return maskedImage
671
672

◆ isTrimmedExposure()

lsst.ip.isr.isrFunctions.isTrimmedExposure ( exposure)
Check if the unused pixels (pre-/over-scan pixels) have
been trimmed from an exposure.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    The exposure to check.

Returns
-------
result : `bool`
    True if the image is trimmed, else False.

Definition at line 1712 of file isrFunctions.py.

1712def isTrimmedExposure(exposure):
1713 """Check if the unused pixels (pre-/over-scan pixels) have
1714 been trimmed from an exposure.
1715
1716 Parameters
1717 ----------
1718 exposure : `lsst.afw.image.Exposure`
1719 The exposure to check.
1720
1721 Returns
1722 -------
1723 result : `bool`
1724 True if the image is trimmed, else False.
1725 """
1726 return exposure.getDetector().getBBox() == exposure.getBBox()
1727
1728

◆ isTrimmedImage()

lsst.ip.isr.isrFunctions.isTrimmedImage ( image,
detector )
Check if the unused pixels (pre-/over-scan pixels) have
been trimmed from an image

Parameters
----------
image : `lsst.afw.image.Image`
    The image to check.
detector : `lsst.afw.cameraGeom.Detector`
    The detector associated with the image.

Returns
-------
result : `bool`
    True if the image is trimmed, else False.

Definition at line 1729 of file isrFunctions.py.

1729def isTrimmedImage(image, detector):
1730 """Check if the unused pixels (pre-/over-scan pixels) have
1731 been trimmed from an image
1732
1733 Parameters
1734 ----------
1735 image : `lsst.afw.image.Image`
1736 The image to check.
1737 detector : `lsst.afw.cameraGeom.Detector`
1738 The detector associated with the image.
1739
1740 Returns
1741 -------
1742 result : `bool`
1743 True if the image is trimmed, else False.
1744 """
1745 return detector.getBBox() == image.getBBox()
1746
1747

◆ makeThresholdMask()

lsst.ip.isr.isrFunctions.makeThresholdMask ( maskedImage,
threshold,
growFootprints = 1,
maskName = 'SAT' )
Mask pixels based on threshold detection.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.  Only the mask plane is updated.
threshold : scalar
    Detection threshold.
growFootprints : scalar, optional
    Number of pixels to grow footprints of detected regions.
maskName : str, optional
    Mask plane name, or list of names to convert

Returns
-------
defectList : `lsst.meas.algorithms.Defects`
    Defect list constructed from pixels above the threshold.

Definition at line 171 of file isrFunctions.py.

171def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName='SAT'):
172 """Mask pixels based on threshold detection.
173
174 Parameters
175 ----------
176 maskedImage : `lsst.afw.image.MaskedImage`
177 Image to process. Only the mask plane is updated.
178 threshold : scalar
179 Detection threshold.
180 growFootprints : scalar, optional
181 Number of pixels to grow footprints of detected regions.
182 maskName : str, optional
183 Mask plane name, or list of names to convert
184
185 Returns
186 -------
187 defectList : `lsst.meas.algorithms.Defects`
188 Defect list constructed from pixels above the threshold.
189 """
190 # find saturated regions
191 thresh = afwDetection.Threshold(threshold)
192 fs = afwDetection.FootprintSet(maskedImage, thresh)
193
194 if growFootprints > 0:
195 fs = afwDetection.FootprintSet(fs, rGrow=growFootprints, isotropic=False)
196 fpList = fs.getFootprints()
197
198 # set mask
199 mask = maskedImage.getMask()
200 bitmask = mask.getPlaneBitMask(maskName)
201 afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
202
203 return Defects.fromFootprintList(fpList)
204
205

◆ maskE2VEdgeBleed()

lsst.ip.isr.isrFunctions.maskE2VEdgeBleed ( exposure,
e2vEdgeBleedSatMinArea = 10000,
e2vEdgeBleedSatMaxArea = 100000,
e2vEdgeBleedYMax = 350,
saturatedMaskName = "SAT",
log = None )
Mask edge bleeds in E2V detectors.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to apply masking to.
e2vEdgeBleedSatMinArea : `int`, optional
    Minimum limit of saturated cores footprint area.
e2vEdgeBleedSatMaxArea : `int`, optional
    Maximum limit of saturated cores footprint area.
e2vEdgeBleedYMax: `float`, optional
    Height of edge bleed masking.
saturatedMaskName : `str`, optional
    Mask name for saturation.
log : `logging.Logger`, optional
    Logger to handle messages.

Definition at line 229 of file isrFunctions.py.

232 saturatedMaskName="SAT", log=None):
233 """Mask edge bleeds in E2V detectors.
234
235 Parameters
236 ----------
237 exposure : `lsst.afw.image.Exposure`
238 Exposure to apply masking to.
239 e2vEdgeBleedSatMinArea : `int`, optional
240 Minimum limit of saturated cores footprint area.
241 e2vEdgeBleedSatMaxArea : `int`, optional
242 Maximum limit of saturated cores footprint area.
243 e2vEdgeBleedYMax: `float`, optional
244 Height of edge bleed masking.
245 saturatedMaskName : `str`, optional
246 Mask name for saturation.
247 log : `logging.Logger`, optional
248 Logger to handle messages.
249 """
250
251 log = log if log else logging.getLogger(__name__)
252
253 maskedImage = exposure.maskedImage
254 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
255
256 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
257
258 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
259
260 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
261 largeAreas, = numpy.where((satAreas >= e2vEdgeBleedSatMinArea)
262 & (satAreas < e2vEdgeBleedSatMaxArea))
263 for largeAreasIndex in largeAreas:
264 fpCore = fpList[largeAreasIndex]
265 xCore, yCore = fpCore.getCentroid()
266 xCore = int(xCore)
267 yCore = int(yCore)
268
269 for amp in exposure.getDetector():
270 if amp.getBBox().contains(xCore, yCore):
271 ampName = amp.getName()
272 if ampName[:2] == 'C0':
273 # Check that the footprint reaches the bottom of the
274 # amplifier.
275 if fpCore.getBBox().getMinY() == 0:
276 # This is a large saturation footprint that hits the
277 # edge, and is thus classified as an edge bleed.
278
279 # TODO DM-50587: Optimize number of rows to mask by
280 # looking at the median signal level as a function of
281 # row number on the right side of the saturation trail.
282
283 log.info("Found E2V edge bleed in amp %s, column %d.", ampName, xCore)
284 maskedImage.mask[amp.getBBox()].array[:e2vEdgeBleedYMax, :] |= saturatedBit
285
286

◆ maskITLDip()

lsst.ip.isr.isrFunctions.maskITLDip ( exposure,
detectorConfig,
maskPlaneNames = ["SUSPECT", "ITL_DIP"],
log = None )
Add mask bits according to the ITL dip model.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to do ITL dip masking.
detectorConfig : `lsst.ip.isr.overscanAmpConfig.OverscanDetectorConfig`
    Configuration for this detector.
maskPlaneNames : `list [`str`], optional
    Name of the ITL Dip mask planes.
log : `logging.Logger`, optional
    If not set, a default logger will be used.

Definition at line 548 of file isrFunctions.py.

548def maskITLDip(exposure, detectorConfig, maskPlaneNames=["SUSPECT", "ITL_DIP"], log=None):
549 """Add mask bits according to the ITL dip model.
550
551 Parameters
552 ----------
553 exposure : `lsst.afw.image.Exposure`
554 Exposure to do ITL dip masking.
555 detectorConfig : `lsst.ip.isr.overscanAmpConfig.OverscanDetectorConfig`
556 Configuration for this detector.
557 maskPlaneNames : `list [`str`], optional
558 Name of the ITL Dip mask planes.
559 log : `logging.Logger`, optional
560 If not set, a default logger will be used.
561 """
562 if detectorConfig.itlDipBackgroundFraction == 0.0:
563 # Nothing to do.
564 return
565
566 if log is None:
567 log = logging.getLogger(__name__)
568
569 thresh = afwDetection.Threshold(
570 exposure.mask.getPlaneBitMask("SAT"),
571 afwDetection.Threshold.BITMASK,
572 )
573 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
574
575 heights = numpy.asarray([fp.getBBox().getHeight() for fp in fpList])
576
577 largeHeights, = numpy.where(heights >= detectorConfig.itlDipMinHeight)
578
579 if len(largeHeights) == 0:
580 return
581
582 # Get the approximate image background.
583 approxBackground = numpy.median(exposure.image.array)
584 maskValue = exposure.mask.getPlaneBitMask(maskPlaneNames)
585
586 maskBak = exposure.mask.array.copy()
587 nMaskedCols = 0
588
589 for index in largeHeights:
590 fp = fpList[index]
591 center = fp.getCentroid()
592
593 nSat = numpy.sum(fp.getSpans().asArray(), axis=0)
594 width = numpy.sum(nSat > detectorConfig.itlDipMinHeight)
595
596 if width < detectorConfig.itlDipMinWidth:
597 continue
598
599 width = numpy.clip(width, None, detectorConfig.itlDipMaxWidth)
600
601 dipMax = detectorConfig.itlDipBackgroundFraction * approxBackground * width
602
603 # Assume sky-noise dominated; we could add in read noise here.
604 if dipMax < detectorConfig.itlDipMinBackgroundNoiseFraction * numpy.sqrt(approxBackground):
605 continue
606
607 minCol = int(center.getX() - (detectorConfig.itlDipWidthScale * width) / 2.)
608 maxCol = int(center.getX() + (detectorConfig.itlDipWidthScale * width) / 2.)
609 minCol = numpy.clip(minCol, 0, None)
610 maxCol = numpy.clip(maxCol, None, exposure.mask.array.shape[1] - 1)
611
612 log.info(
613 "Found ITL dip (width %d; bkg %.2f); masking column %d to %d.",
614 width,
615 approxBackground,
616 minCol,
617 maxCol,
618 )
619
620 exposure.mask.array[:, minCol: maxCol + 1] |= maskValue
621
622 nMaskedCols += (maxCol - minCol + 1)
623
624 if nMaskedCols > detectorConfig.itlDipMaxColsPerImage:
625 log.warning(
626 "Too many (%d) columns would be masked on this image from dip masking; restoring original mask.",
627 nMaskedCols,
628 )
629 exposure.mask.array[:, :] = maskBak
630
631

◆ maskITLEdgeBleed()

lsst.ip.isr.isrFunctions.maskITLEdgeBleed ( ccdExposure,
badAmpDict,
fpCore,
itlEdgeBleedSatMinArea = 10000,
itlEdgeBleedSatMaxArea = 100000,
itlEdgeBleedThreshold = 5000.,
itlEdgeBleedModelConstant = 0.02,
saturatedMaskName = "SAT",
log = None )
Mask edge bleeds in ITL detectors.

Parameters
----------
ccdExposure : `lsst.afw.image.Exposure`
    Exposure to apply masking to.
badAmpDict : `dict` [`str`, `bool`]
    Dictionary of amplifiers, keyed by name, value is True if
    amplifier is fully masked.
fpCore : `lsst.afw.detection._detection.Footprint`
    Footprint of saturated core.
itlEdgeBleedThreshold : `float`, optional
    Threshold above median sky background for edge bleed detection
    (electron units).
itlEdgeBleedModelConstant : `float`, optional
    Constant in the decaying exponential in the edge bleed masking.
saturatedMaskName : `str`, optional
    Mask name for saturation.
log : `logging.Logger`, optional
    Logger to handle messages.

Definition at line 287 of file isrFunctions.py.

292 saturatedMaskName="SAT", log=None):
293 """Mask edge bleeds in ITL detectors.
294
295 Parameters
296 ----------
297 ccdExposure : `lsst.afw.image.Exposure`
298 Exposure to apply masking to.
299 badAmpDict : `dict` [`str`, `bool`]
300 Dictionary of amplifiers, keyed by name, value is True if
301 amplifier is fully masked.
302 fpCore : `lsst.afw.detection._detection.Footprint`
303 Footprint of saturated core.
304 itlEdgeBleedThreshold : `float`, optional
305 Threshold above median sky background for edge bleed detection
306 (electron units).
307 itlEdgeBleedModelConstant : `float`, optional
308 Constant in the decaying exponential in the edge bleed masking.
309 saturatedMaskName : `str`, optional
310 Mask name for saturation.
311 log : `logging.Logger`, optional
312 Logger to handle messages.
313 """
314
315 log = log if log else logging.getLogger(__name__)
316
317 # Get median of amplifier saturation level
318 satLevel = numpy.nanmedian([ccdExposure.metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
319 for amp in ccdExposure.getDetector() if not badAmpDict[amp.getName()]])
320
321 # 1. we check if there are several cores in the footprint:
322 # Get centroid of saturated core
323 xCore, yCore = fpCore.getCentroid()
324 # Turn the Y detector coordinate into Y footprint coordinate
325 yCoreFP = int(yCore) - fpCore.getBBox().getMinY()
326 # Now test if there is one or more cores by checking if the slice at the
327 # center is full of saturated pixels or has several segments of saturated
328 # columns (i.e. several cores with trails)
329 checkCoreNbRow = fpCore.getSpans().asArray()[yCoreFP, :]
330 nbCore = 0
331 indexSwitchTrue = []
332 indexSwitchFalse = []
333 if checkCoreNbRow[0]:
334 # If the slice starts with saturated pixels
335 inSatSegment = True
336 nbCore = 1
337 indexSwitchTrue.append(0)
338 else:
339 # If the slice starts with non saturated pixels
340 inSatSegment = False
341
342 for i, value in enumerate(checkCoreNbRow):
343 if value:
344 if not inSatSegment:
345 indexSwitchTrue.append(i)
346 # nbCore is the number of detected cores.
347 nbCore += 1
348 inSatSegment = True
349 elif inSatSegment:
350 indexSwitchFalse.append(i)
351 inSatSegment = False
352
353 # 1. we look for edge bleed in saturated cores in the footprint
354 if nbCore == 2:
355 # we now estimate the x coordinates of the edges of the subfootprint
356 # for each core
357 xEdgesCores = [0]
358 xEdgesCores.append(int((indexSwitchTrue[1] + indexSwitchFalse[0])/2))
359 xEdgesCores.append(fpCore.getSpans().asArray().shape[1])
360 # Get the X and Y footprint coordinates of the cores
361 for i in range(nbCore):
362 subfp = fpCore.getSpans().asArray()[:, xEdgesCores[i]:xEdgesCores[i+1]]
363 xCoreFP = int(xEdgesCores[i] + numpy.argmax(numpy.sum(subfp, axis=0)))
364 # turn into X coordinate in detector space
365 xCore = xCoreFP + fpCore.getBBox().getMinX()
366 # get Y footprint coordinate of the core
367 # by trimming the edges where edge bleeds are potentially dominant
368 if subfp.shape[0] <= 200:
369 yCoreFP = int(numpy.argmax(numpy.sum(subfp, axis=1)))
370 else:
371 yCoreFP = int(numpy.argmax(numpy.sum(subfp[100:-100, :],
372 axis=1)))
373 yCoreFP = 100+yCoreFP
374
375 # Estimate the width of the saturated core
376 widthSat = numpy.sum(subfp[int(yCoreFP), :])
377
378 subfpArea = numpy.sum(subfp)
379 if subfpArea > itlEdgeBleedSatMinArea and subfpArea < itlEdgeBleedSatMaxArea:
380 _applyMaskITLEdgeBleed(ccdExposure, xCore,
381 satLevel, widthSat,
382 itlEdgeBleedThreshold,
383 itlEdgeBleedModelConstant,
384 saturatedMaskName, log)
385 elif nbCore > 2:
386 # TODO DM-49736: support N cores in saturated footprint
387 log.warning(
388 "Too many (%d) cores in saturated footprint to mask edge bleeds.",
389 nbCore,
390 )
391 else:
392 # Get centroid of saturated core
393 xCore, yCore = fpCore.getCentroid()
394 # Turn the Y detector coordinate into Y footprint coordinate
395 yCoreFP = yCore - fpCore.getBBox().getMinY()
396 # Get the number of saturated columns around the centroid
397 widthSat = numpy.sum(fpCore.getSpans().asArray()[int(yCoreFP), :])
398 _applyMaskITLEdgeBleed(ccdExposure, xCore,
399 satLevel, widthSat, itlEdgeBleedThreshold,
400 itlEdgeBleedModelConstant, saturatedMaskName, log)
401
402

◆ maskITLSatSag()

lsst.ip.isr.isrFunctions.maskITLSatSag ( ccdExposure,
fpCore,
saturatedMaskName = "SAT" )
Mask columns presenting saturation sag in saturated footprints in
ITL detectors.

Parameters
----------
ccdExposure : `lsst.afw.image.Exposure`
    Exposure to apply masking to.
fpCore : `lsst.afw.detection._detection.Footprint`
    Footprint of saturated core.
saturatedMaskName : `str`, optional
    Mask name for saturation.

Definition at line 520 of file isrFunctions.py.

520def maskITLSatSag(ccdExposure, fpCore, saturatedMaskName="SAT"):
521 """Mask columns presenting saturation sag in saturated footprints in
522 ITL detectors.
523
524 Parameters
525 ----------
526 ccdExposure : `lsst.afw.image.Exposure`
527 Exposure to apply masking to.
528 fpCore : `lsst.afw.detection._detection.Footprint`
529 Footprint of saturated core.
530 saturatedMaskName : `str`, optional
531 Mask name for saturation.
532 """
533
534 # TODO DM-49736: add a flux level check to apply masking
535
536 maskedImage = ccdExposure.maskedImage
537 saturatedBit = maskedImage.mask.getPlaneBitMask(saturatedMaskName)
538
539 cc = numpy.sum(fpCore.getSpans().asArray(), axis=0)
540 # Mask full columns that have 20 percent of the height of the footprint
541 # saturated
542 columnsToMaskFP = numpy.where(cc > fpCore.getSpans().asArray().shape[0]/5.)
543
544 columnsToMask = [x + int(fpCore.getBBox().getMinX()) for x in columnsToMaskFP]
545 maskedImage.mask.array[:, columnsToMask] |= saturatedBit
546
547

◆ saturationCorrection()

lsst.ip.isr.isrFunctions.saturationCorrection ( maskedImage,
saturation,
fwhm,
growFootprints = 1,
interpolate = True,
maskName = 'SAT',
fallbackValue = None,
useLegacyInterp = True )
Mark saturated pixels and optionally interpolate over them

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.
saturation  : scalar
    Saturation level used as the detection threshold.
fwhm : `float`
    FWHM of double Gaussian smoothing kernel.
growFootprints : scalar, optional
    Number of pixels to grow footprints of detected regions.
interpolate : Bool, optional
    If True, saturated pixels are interpolated over.
maskName : str, optional
    Mask plane name.
fallbackValue : scalar, optional
    Value of last resort for interpolation.

Notes
-----
The ``fwhm`` parameter is used to create a PSF, but the underlying
interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
not currently make use of this information.

Definition at line 673 of file isrFunctions.py.

674 fallbackValue=None, useLegacyInterp=True):
675 """Mark saturated pixels and optionally interpolate over them
676
677 Parameters
678 ----------
679 maskedImage : `lsst.afw.image.MaskedImage`
680 Image to process.
681 saturation : scalar
682 Saturation level used as the detection threshold.
683 fwhm : `float`
684 FWHM of double Gaussian smoothing kernel.
685 growFootprints : scalar, optional
686 Number of pixels to grow footprints of detected regions.
687 interpolate : Bool, optional
688 If True, saturated pixels are interpolated over.
689 maskName : str, optional
690 Mask plane name.
691 fallbackValue : scalar, optional
692 Value of last resort for interpolation.
693
694 Notes
695 -----
696 The ``fwhm`` parameter is used to create a PSF, but the underlying
697 interpolation code (`lsst.meas.algorithms.interpolateOverDefects`) does
698 not currently make use of this information.
699 """
700 defectList = makeThresholdMask(
701 maskedImage=maskedImage,
702 threshold=saturation,
703 growFootprints=growFootprints,
704 maskName=maskName,
705 )
706 if interpolate:
707 interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue,
708 maskNameList=[maskName], useLegacyInterp=useLegacyInterp)
709
710 return maskedImage
711
712

◆ setBadRegions()

lsst.ip.isr.isrFunctions.setBadRegions ( exposure,
badStatistic = "MEDIAN" )
Set all BAD areas of the chip to the average of the rest of the exposure

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to mask.  The exposure mask is modified.
badStatistic : `str`, optional
    Statistic to use to generate the replacement value from the
    image data.  Allowed values are 'MEDIAN' or 'MEANCLIP'.

Returns
-------
badPixelCount : scalar
    Number of bad pixels masked.
badPixelValue : scalar
    Value substituted for bad pixels.

Raises
------
RuntimeError
    Raised if `badStatistic` is not an allowed value.

Definition at line 1509 of file isrFunctions.py.

1509def setBadRegions(exposure, badStatistic="MEDIAN"):
1510 """Set all BAD areas of the chip to the average of the rest of the exposure
1511
1512 Parameters
1513 ----------
1514 exposure : `lsst.afw.image.Exposure`
1515 Exposure to mask. The exposure mask is modified.
1516 badStatistic : `str`, optional
1517 Statistic to use to generate the replacement value from the
1518 image data. Allowed values are 'MEDIAN' or 'MEANCLIP'.
1519
1520 Returns
1521 -------
1522 badPixelCount : scalar
1523 Number of bad pixels masked.
1524 badPixelValue : scalar
1525 Value substituted for bad pixels.
1526
1527 Raises
1528 ------
1529 RuntimeError
1530 Raised if `badStatistic` is not an allowed value.
1531 """
1532 if badStatistic == "MEDIAN":
1533 statistic = afwMath.MEDIAN
1534 elif badStatistic == "MEANCLIP":
1535 statistic = afwMath.MEANCLIP
1536 else:
1537 raise RuntimeError("Impossible method %s of bad region correction" % badStatistic)
1538
1539 mi = exposure.getMaskedImage()
1540 mask = mi.getMask()
1541 BAD = mask.getPlaneBitMask("BAD")
1542 INTRP = mask.getPlaneBitMask("INTRP")
1543
1545 sctrl.setAndMask(BAD)
1546 value = afwMath.makeStatistics(mi, statistic, sctrl).getValue()
1547
1548 maskArray = mask.getArray()
1549 imageArray = mi.getImage().getArray()
1550 badPixels = numpy.logical_and((maskArray & BAD) > 0, (maskArray & INTRP) == 0)
1551 imageArray[:] = numpy.where(badPixels, value, imageArray)
1552
1553 return badPixels.sum(), value
1554
1555
Pass parameters to a Statistics object.
Definition Statistics.h:83

◆ transferFlux()

lsst.ip.isr.isrFunctions.transferFlux ( cFunc,
fStep,
correctionMode = True )
Take the input convolved deflection potential and the flux array
to compute and apply the flux transfer into the correction array.

Parameters
----------
cFunc: `numpy.array`
    Deflection potential, being the convolution of the flux F with the
    kernel K.
fStep: `numpy.array`
    The array of flux values which act as the source of the flux transfer.
correctionMode: `bool`
    Defines if applying correction (True) or generating sims (False).

Returns
-------
corr:
    BFE correction array

Definition at line 1070 of file isrFunctions.py.

1070def transferFlux(cFunc, fStep, correctionMode=True):
1071 """Take the input convolved deflection potential and the flux array
1072 to compute and apply the flux transfer into the correction array.
1073
1074 Parameters
1075 ----------
1076 cFunc: `numpy.array`
1077 Deflection potential, being the convolution of the flux F with the
1078 kernel K.
1079 fStep: `numpy.array`
1080 The array of flux values which act as the source of the flux transfer.
1081 correctionMode: `bool`
1082 Defines if applying correction (True) or generating sims (False).
1083
1084 Returns
1085 -------
1086 corr:
1087 BFE correction array
1088 """
1089
1090 if cFunc.shape != fStep.shape:
1091 raise RuntimeError(f'transferFlux: array shapes do not match: {cFunc.shape}, {fStep.shape}')
1092
1093 # set the sign of the correction and set its value for the
1094 # time averaged solution
1095 if correctionMode:
1096 # negative sign if applying BFE correction
1097 factor = -0.5
1098 else:
1099 # positive sign if generating BFE simulations
1100 factor = 0.5
1101
1102 # initialise the BFE correction image to zero
1103 corr = numpy.zeros(cFunc.shape, dtype=numpy.float64)
1104
1105 # Generate a 2D mesh of x,y coordinates
1106 yDim, xDim = cFunc.shape
1107 y = numpy.arange(yDim, dtype=int)
1108 x = numpy.arange(xDim, dtype=int)
1109 xc, yc = numpy.meshgrid(x, y)
1110
1111 # process each axis in turn
1112 for ax in [0, 1]:
1113
1114 # gradient of phi on right/upper edge of pixel
1115 diff = numpy.diff(cFunc, axis=ax)
1116
1117 # expand array back to full size with zero gradient at the end
1118 gx = numpy.zeros(cFunc.shape, dtype=numpy.float64)
1119 yDiff, xDiff = diff.shape
1120 gx[:yDiff, :xDiff] += diff
1121
1122 # select pixels with either positive gradients on the right edge,
1123 # flux flowing to the right/up
1124 # or negative gradients, flux flowing to the left/down
1125 for i, sel in enumerate([gx > 0, gx < 0]):
1126 xSelPixels = xc[sel]
1127 ySelPixels = yc[sel]
1128 # and add the flux into the pixel to the right or top
1129 # depending on which axis we are handling
1130 if ax == 0:
1131 xPix = xSelPixels
1132 yPix = ySelPixels+1
1133 else:
1134 xPix = xSelPixels+1
1135 yPix = ySelPixels
1136 # define flux as the either current pixel value or pixel
1137 # above/right
1138 # depending on whether positive or negative gradient
1139 if i == 0:
1140 # positive gradients, flux flowing to higher coordinate values
1141 flux = factor * fStep[sel]*gx[sel]
1142 else:
1143 # negative gradients, flux flowing to lower coordinate values
1144 flux = factor * fStep[yPix, xPix]*gx[sel]
1145 # change the fluxes of the donor and receiving pixels
1146 # such that flux is conserved
1147 corr[sel] -= flux
1148 corr[yPix, xPix] += flux
1149
1150 # return correction array
1151 return corr
1152
1153

◆ transposeMaskedImage()

lsst.ip.isr.isrFunctions.transposeMaskedImage ( maskedImage)
Make a transposed copy of a masked image.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.

Returns
-------
transposed : `lsst.afw.image.MaskedImage`
    The transposed copy of the input image.

Definition at line 94 of file isrFunctions.py.

94def transposeMaskedImage(maskedImage):
95 """Make a transposed copy of a masked image.
96
97 Parameters
98 ----------
99 maskedImage : `lsst.afw.image.MaskedImage`
100 Image to process.
101
102 Returns
103 -------
104 transposed : `lsst.afw.image.MaskedImage`
105 The transposed copy of the input image.
106 """
107 transposed = maskedImage.Factory(lsst.geom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
108 transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
109 transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
110 transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
111 return transposed
112
113

◆ trimToMatchCalibBBox()

lsst.ip.isr.isrFunctions.trimToMatchCalibBBox ( rawMaskedImage,
calibMaskedImage )
Compute number of edge trim pixels to match the calibration data.

Use the dimension difference between the raw exposure and the
calibration exposure to compute the edge trim pixels.  This trim
is applied symmetrically, with the same number of pixels masked on
each side.

Parameters
----------
rawMaskedImage : `lsst.afw.image.MaskedImage`
    Image to trim.
calibMaskedImage : `lsst.afw.image.MaskedImage`
    Calibration image to draw new bounding box from.

Returns
-------
replacementMaskedImage : `lsst.afw.image.MaskedImage`
    ``rawMaskedImage`` trimmed to the appropriate size.

Raises
------
RuntimeError
   Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
   match ``calibMaskedImage``.

Definition at line 713 of file isrFunctions.py.

713def trimToMatchCalibBBox(rawMaskedImage, calibMaskedImage):
714 """Compute number of edge trim pixels to match the calibration data.
715
716 Use the dimension difference between the raw exposure and the
717 calibration exposure to compute the edge trim pixels. This trim
718 is applied symmetrically, with the same number of pixels masked on
719 each side.
720
721 Parameters
722 ----------
723 rawMaskedImage : `lsst.afw.image.MaskedImage`
724 Image to trim.
725 calibMaskedImage : `lsst.afw.image.MaskedImage`
726 Calibration image to draw new bounding box from.
727
728 Returns
729 -------
730 replacementMaskedImage : `lsst.afw.image.MaskedImage`
731 ``rawMaskedImage`` trimmed to the appropriate size.
732
733 Raises
734 ------
735 RuntimeError
736 Raised if ``rawMaskedImage`` cannot be symmetrically trimmed to
737 match ``calibMaskedImage``.
738 """
739 nx, ny = rawMaskedImage.getBBox().getDimensions() - calibMaskedImage.getBBox().getDimensions()
740 if nx != ny:
741 raise RuntimeError("Raw and calib maskedImages are trimmed differently in X and Y.")
742 if nx % 2 != 0:
743 raise RuntimeError("Calibration maskedImage is trimmed unevenly in X.")
744 if nx < 0:
745 raise RuntimeError("Calibration maskedImage is larger than raw data.")
746
747 nEdge = nx//2
748 if nEdge > 0:
749 replacementMaskedImage = rawMaskedImage[nEdge:-nEdge, nEdge:-nEdge, afwImage.LOCAL]
750 SourceDetectionTask.setEdgeBits(
751 rawMaskedImage,
752 replacementMaskedImage.getBBox(),
753 rawMaskedImage.getMask().getPlaneBitMask("EDGE")
754 )
755 else:
756 replacementMaskedImage = rawMaskedImage
757
758 return replacementMaskedImage
759
760

◆ updateVariance()

lsst.ip.isr.isrFunctions.updateVariance ( maskedImage,
gain,
readNoise,
replace = True )
Set the variance plane based on the image plane.

The maskedImage must have units of `adu` (if gain != 1.0) or
electron (if gain == 1.0). This routine will always produce a
variance plane in the same units as the image.

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
    Image to process.  The variance plane is modified.
gain : scalar
    The amplifier gain in electron/adu.
readNoise : scalar
    The amplifier read noise in electron/pixel.
replace : `bool`, optional
    Replace the current variance?  If False, the image
    variance will be added to the current variance plane.

Definition at line 834 of file isrFunctions.py.

834def updateVariance(maskedImage, gain, readNoise, replace=True):
835 """Set the variance plane based on the image plane.
836
837 The maskedImage must have units of `adu` (if gain != 1.0) or
838 electron (if gain == 1.0). This routine will always produce a
839 variance plane in the same units as the image.
840
841 Parameters
842 ----------
843 maskedImage : `lsst.afw.image.MaskedImage`
844 Image to process. The variance plane is modified.
845 gain : scalar
846 The amplifier gain in electron/adu.
847 readNoise : scalar
848 The amplifier read noise in electron/pixel.
849 replace : `bool`, optional
850 Replace the current variance? If False, the image
851 variance will be added to the current variance plane.
852 """
853 var = maskedImage.variance
854 if replace:
855 var[:, :] = maskedImage.image
856 else:
857 var[:, :] += maskedImage.image
858 var /= gain
859 var += (readNoise/gain)**2
860
861

◆ widenSaturationTrails()

lsst.ip.isr.isrFunctions.widenSaturationTrails ( mask)
Grow the saturation trails by an amount dependent on the width of the
trail.

Parameters
----------
mask : `lsst.afw.image.Mask`
    Mask which will have the saturated areas grown.

Definition at line 1461 of file isrFunctions.py.

1461def widenSaturationTrails(mask):
1462 """Grow the saturation trails by an amount dependent on the width of the
1463 trail.
1464
1465 Parameters
1466 ----------
1467 mask : `lsst.afw.image.Mask`
1468 Mask which will have the saturated areas grown.
1469 """
1470
1471 extraGrowDict = {}
1472 for i in range(1, 6):
1473 extraGrowDict[i] = 0
1474 for i in range(6, 8):
1475 extraGrowDict[i] = 1
1476 for i in range(8, 10):
1477 extraGrowDict[i] = 3
1478 extraGrowMax = 4
1479
1480 if extraGrowMax <= 0:
1481 return
1482
1483 saturatedBit = mask.getPlaneBitMask("SAT")
1484
1485 xmin, ymin = mask.getBBox().getMin()
1486 width = mask.getWidth()
1487
1488 thresh = afwDetection.Threshold(saturatedBit, afwDetection.Threshold.BITMASK)
1489 fpList = afwDetection.FootprintSet(mask, thresh).getFootprints()
1490
1491 for fp in fpList:
1492 for s in fp.getSpans():
1493 x0, x1 = s.getX0(), s.getX1()
1494
1495 extraGrow = extraGrowDict.get(x1 - x0 + 1, extraGrowMax)
1496 if extraGrow > 0:
1497 y = s.getY() - ymin
1498 x0 -= xmin + extraGrow
1499 x1 -= xmin - extraGrow
1500
1501 if x0 < 0:
1502 x0 = 0
1503 if x1 >= width - 1:
1504 x1 = width - 1
1505
1506 mask.array[y, x0:x1+1] |= saturatedBit
1507
1508