Skip to content

Commit

Permalink
Add realistic brighter-fatter to IsrMockLSST and add it to IsrTaskLSST
Browse files Browse the repository at this point in the history
  • Loading branch information
Alex Broughton committed Aug 15, 2024
1 parent 23a948e commit 2dfb481
Show file tree
Hide file tree
Showing 5 changed files with 252 additions and 37 deletions.
1 change: 1 addition & 0 deletions python/lsst/ip/isr/isrMock.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,7 @@ def __init__(self, **kwargs):
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
# Generic gaussian BF kernel
self.bfKernel = np.array([[1., 4., 7., 4., 1.],
[4., 16., 26., 16., 4.],
[7., 26., 41., 26., 7.],
Expand Down
183 changes: 181 additions & 2 deletions python/lsst/ip/isr/isrMockLSST.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
"TransmissionMockLSST"]

import numpy as np
import galsim

import lsst.geom as geom
import lsst.pex.config as pexConfig
Expand All @@ -36,6 +37,7 @@
from .defects import Defects
from .assembleCcdTask import AssembleCcdTask
from .linearize import Linearizer
from .brighterFatterKernel import BrighterFatterKernel


class IsrMockLSSTConfig(IsrMockConfig):
Expand Down Expand Up @@ -68,6 +70,22 @@ class IsrMockLSSTConfig(IsrMockConfig):
default=30000.0,
doc="Bright defect level (electron).",
)
doAddBrighterFatter = pexConfig.Field(
dtype=bool,
default=False,
doc="Add brighter fatter and/or diffusion effects to image.",
)
bfStrength = pexConfig.Field(
dtype=float,
default=2.0,
doc="The brighter fatter effect scaling parameter (cannot be zero)."
"Nominally = 1, but = 2 is more realistic."
)
nRecalc = pexConfig.Field(
dtype=int,
default=10000,
doc="Number of electrons to accumulate before recalculating pixel shapes.",
)
doAddClockInjectedOffset = pexConfig.Field(
dtype=bool,
default=True,
Expand Down Expand Up @@ -168,9 +186,116 @@ class IsrMockLSST(IsrMock):
_DefaultName = "isrMockLSST"

def __init__(self, **kwargs):
# cross-talk coeffs, bf kernel are defined in the parent class.
super().__init__(**kwargs)

# Get kernel derived from imSim generated flats with BFE. The kernel

Check failure on line 191 in python/lsst/ip/isr/isrMockLSST.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

W291 

trailing whitespace
# was used for Ops Rehearsal 3 for LSSTCam-type sensors
# See https://rubinobs.atlassian.net/browse/DM-43059 for more details.
self.bfKernel = np.array([[4.83499829e-01, 8.10171823e-01, 5.31096720e-01,
3.54369868e-02, -8.44782871e-01, -1.64614462e+00,
-3.83933101e+00, -5.60243416e+00, -6.51691578e+00,
-5.60243416e+00, -3.83933101e+00, -1.64614462e+00,
-8.44782871e-01, 3.54369868e-02, 5.31096720e-01,
8.10171823e-01, 4.83499829e-01],
[1.12382749e+00, 2.22609074e+00, 1.27877807e+00,
4.55434098e-01, -1.76842385e+00, -1.90046460e+00,
-8.10874526e+00, -1.20534899e+01, -1.48627948e+01,
-1.20534899e+01, -8.10874526e+00, -1.90046460e+00,
-1.76842385e+00, 4.55434098e-01, 1.27877807e+00,
2.22609074e+00, 1.12382749e+00],
[1.78571940e+00, 4.38918110e+00, 3.95098587e+00,
3.70961649e-01, -3.48151981e+00, -9.61567736e+00,
-1.78621172e+01, -2.32278872e+01, -2.31833727e+01,
-2.32278872e+01, -1.78621172e+01, -9.61567736e+00,
-3.48151981e+00, 3.70961649e-01, 3.95098587e+00,
4.38918110e+00, 1.78571940e+00],
[1.62986900e+00, 3.67851228e+00, 5.68645252e+00,
2.15342566e-01, -8.89937202e+00, -1.44739813e+01,
-2.98952660e+01, -4.37420817e+01, -4.83160958e+01,
-4.37420817e+01, -2.98952660e+01, -1.44739813e+01,
-8.89937202e+00, 2.15342566e-01, 5.68645252e+00,
3.67851228e+00, 1.62986900e+00],
[1.05524430e+00, 1.71917897e+00, 1.73105590e+00,
-2.10088420e+00, -1.15118208e+01, -2.55007598e+01,
-4.73056159e+01, -6.97257685e+01, -8.09264433e+01,
-6.97257685e+01, -4.73056159e+01, -2.55007598e+01,
-1.15118208e+01, -2.10088420e+00, 1.73105590e+00,
1.71917897e+00, 1.05524430e+00],
[8.71929228e-01, 5.41025574e-01, 9.47560771e-01,
-5.75314708e-01, -7.46104027e+00, -4.42314481e+01,
-9.54126971e+01, -1.61603201e+02, -2.07520692e+02,
-1.61603201e+02, -9.54126971e+01, -4.42314481e+01,
-7.46104027e+00, -5.75314708e-01, 9.47560771e-01,
5.41025574e-01, 8.71929228e-01],
[1.89144704e+00, 3.57543979e+00, -6.91419168e-02,
-3.37950835e+00, -1.46695089e+01, -7.22850746e+01,
-1.65563055e+02, -3.10820425e+02, -4.70026655e+02,
-3.10820425e+02, -1.65563055e+02, -7.22850746e+01,
-1.46695089e+01, -3.37950835e+00, -6.91419168e-02,
3.57543979e+00, 1.89144704e+00],
[3.11841913e+00, 7.84024994e+00, 1.88495248e+00,
-7.69011009e+00, -2.71782400e+01, -1.04343326e+02,
-2.47561370e+02, -5.32959841e+02, -1.16529012e+03,
-5.32959841e+02, -2.47561370e+02, -1.04343326e+02,
-2.71782400e+01, -7.69011009e+00, 1.88495248e+00,
7.84024994e+00, 3.11841913e+00],
[2.74197956e+00, 4.73107997e+00, -9.48352966e-01,
-9.44822832e+00, -3.06477671e+01, -1.26788739e+02,
-3.22828411e+02, -8.47943472e+02, -3.87702420e+03,
-8.47943472e+02, -3.22828411e+02, -1.26788739e+02,
-3.06477671e+01, -9.44822832e+00, -9.48352966e-01,
4.73107997e+00, 2.74197956e+00],
[3.11841913e+00, 7.84024994e+00, 1.88495248e+00,
-7.69011009e+00, -2.71782400e+01, -1.04343326e+02,
-2.47561370e+02, -5.32959841e+02, -1.16529012e+03,
-5.32959841e+02, -2.47561370e+02, -1.04343326e+02,
-2.71782400e+01, -7.69011009e+00, 1.88495248e+00,
7.84024994e+00, 3.11841913e+00],
[1.89144704e+00, 3.57543979e+00, -6.91419168e-02,
-3.37950835e+00, -1.46695089e+01, -7.22850746e+01,
-1.65563055e+02, -3.10820425e+02, -4.70026655e+02,
-3.10820425e+02, -1.65563055e+02, -7.22850746e+01,
-1.46695089e+01, -3.37950835e+00, -6.91419168e-02,
3.57543979e+00, 1.89144704e+00],
[8.71929228e-01, 5.41025574e-01, 9.47560771e-01,
-5.75314708e-01, -7.46104027e+00, -4.42314481e+01,
-9.54126971e+01, -1.61603201e+02, -2.07520692e+02,
-1.61603201e+02, -9.54126971e+01, -4.42314481e+01,
-7.46104027e+00, -5.75314708e-01, 9.47560771e-01,
5.41025574e-01, 8.71929228e-01],
[1.05524430e+00, 1.71917897e+00, 1.73105590e+00,
-2.10088420e+00, -1.15118208e+01, -2.55007598e+01,
-4.73056159e+01, -6.97257685e+01, -8.09264433e+01,
-6.97257685e+01, -4.73056159e+01, -2.55007598e+01,
-1.15118208e+01, -2.10088420e+00, 1.73105590e+00,
1.71917897e+00, 1.05524430e+00],
[1.62986900e+00, 3.67851228e+00, 5.68645252e+00,
2.15342566e-01, -8.89937202e+00, -1.44739813e+01,
-2.98952660e+01, -4.37420817e+01, -4.83160958e+01,
-4.37420817e+01, -2.98952660e+01, -1.44739813e+01,
-8.89937202e+00, 2.15342566e-01, 5.68645252e+00,
3.67851228e+00, 1.62986900e+00],
[1.78571940e+00, 4.38918110e+00, 3.95098587e+00,
3.70961649e-01, -3.48151981e+00, -9.61567736e+00,
-1.78621172e+01, -2.32278872e+01, -2.31833727e+01,
-2.32278872e+01, -1.78621172e+01, -9.61567736e+00,
-3.48151981e+00, 3.70961649e-01, 3.95098587e+00,
4.38918110e+00, 1.78571940e+00],
[1.12382749e+00, 2.22609074e+00, 1.27877807e+00,
4.55434098e-01, -1.76842385e+00, -1.90046460e+00,
-8.10874526e+00, -1.20534899e+01, -1.48627948e+01,
-1.20534899e+01, -8.10874526e+00, -1.90046460e+00,
-1.76842385e+00, 4.55434098e-01, 1.27877807e+00,
2.22609074e+00, 1.12382749e+00],
[4.83499829e-01, 8.10171823e-01, 5.31096720e-01,
3.54369868e-02, -8.44782871e-01, -1.64614462+00,
-3.83933101e+00, -5.60243416e+00, -6.51691578e+00,
-5.60243416e+00, -3.83933101e+00, -1.64614462e+00,
-8.44782871e-01, 3.54369868e-02, 5.31096720e-01,
8.10171823e-01, 4.83499829e-01]]) * 1e-10

# cross-talk coeffs are defined in the parent class.

self.makeSubtask("assembleCcd")

def run(self):
Expand Down Expand Up @@ -220,6 +345,7 @@ def makeImage(self):
rng2DBias = np.random.RandomState(seed=self.config.rngSeed + 3)
rngOverscan = np.random.RandomState(seed=self.config.rngSeed + 4)
rngReadNoise = np.random.RandomState(seed=self.config.rngSeed + 5)
rngBrighterFatter = galsim.BaseDeviate(self.config.rngSeed + 6)

# Create the linearizer if we will need it.
if self.config.doAddHighSignalNonlinearity:
Expand Down Expand Up @@ -322,7 +448,11 @@ def makeImage(self):
self.amplifierAddNoise(ampImageData, darkLevel, darkNoise, rng=rngDark)

# 3. Add BF effect (electron) to imaging portion of the amp.
# TODO
if self.config.doAddBrighterFatter is True:
self.amplifierAddBrighterFatter(ampImageData,
rngBrighterFatter,
self.config.bfStrength,
self.config.nRecalc)

# 4. Add serial CTI (electron) to amplifier (imaging + overscan).
# TODO
Expand Down Expand Up @@ -536,6 +666,55 @@ def makeDefectList(self, isTrimmed=True):

return Defects.fromFootprintList(footprintSet.getFootprints())

def makeBfKernel(self):
"""Generate a simple simulated brighter-fatter kernel.
Returns
-------
kernel : `lsst.ip.isr.BrighterFatterKernel`
Simulated brighter-fatter kernel.
"""
bfkArray = super().makeBfKernel()
bfKernelObject = BrighterFatterKernel()
bfKernelObject.level = 'AMP'
bfKernelObject.gain = self.config.gainDict

for amp in self.getExposure().getDetector():
# Kernel must be in (y,x) orientation
bfKernelObject.ampKernels[amp.getName()] = bfkArray.T

return bfKernelObject

def amplifierAddBrighterFatter(self, ampImageData, rng, bfStrength, nRecalc):
"""Add brighter fatter effect and/or diffusion to the image.
Parameters
----------
ampImageData : `lsst.afw.image.ImageF`
Amplifier image to operate on.
rng : `galsim.BaseDeviate`
Random number generator.
bfStrength : `float`
Scaling parameter of the brighter fatter effect (nominally = 1)
nRecalc: 'int'
The number of electrons to accumulate before recalculating the
distortion of the pixel shapes.
"""

incidentImage = galsim.Image(ampImageData.array, scale=1)
measuredImage = galsim.ImageF(ampImageData.array.shape[1],
ampImageData.array.shape[0],
scale=1)
photons = galsim.PhotonArray.makeFromImage(incidentImage)

sensorModel = galsim.SiliconSensor(strength=bfStrength,
rng=rng,
diffusion_factor=0.0,
nrecalc=nRecalc)

totalFluxAdded = sensorModel.accumulate(photons, measuredImage)
ampImageData.array = measuredImage.array

return totalFluxAdded

def makeLinearizer(self):
# docstring inherited.

Expand Down
67 changes: 35 additions & 32 deletions python/lsst/ip/isr/isrTaskLSST.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -412,11 +412,6 @@ class IsrTaskLSSTConfig(pipeBase.PipelineTaskConfig,
"from the previous iteration summed over all the pixels.",
default=1000,
)
brighterFatterApplyGain = pexConfig.Field(
dtype=bool,
doc="Should the gain be applied when applying the brighter-fatter correction?",
default=True,
)
brighterFatterMaskListToInterpolate = pexConfig.ListField(
dtype=str,
doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
Expand Down Expand Up @@ -1061,39 +1056,29 @@ def getBrighterFatterKernel(self, detector, bfKernel):
else:
raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
elif bfKernel.level == 'AMP':
self.log.warning("Making DETECTOR level kernel from AMP based brighter "
"fatter kernels.")
self.log.info("Making DETECTOR level kernel from AMP based brighter "
"fatter kernels.")
bfKernel.makeDetectorKernelFromAmpwiseKernels(detName)
bfKernelOut = bfKernel.detKernels[detName]
return bfKernelOut, bfGains

def applyBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel, bfGains):
# We need to apply flats and darks before we can interpolate, and
# we need to interpolate before we do B-F, but we do B-F without
# the flats and darks applied so we can work in units of electrons
# or holes. This context manager applies and then removes the darks
# and flats.
#
# We also do not want to interpolate values here, so operate on
# temporary images so we can apply only the BF-correction and roll
# back the interpolation.
# This won't be necessary once the gain normalization
# is done appropriately.
def applyBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel, brighterFatterApplyGain,
bfGains):
# We need to interpolate before we do B-F. Note th
# brighterFatterFwhmForInterpolation is currently unused.
interpExp = ccdExposure.clone()
with self.flatContext(interpExp, flat, dark):
isrFunctions.interpolateFromMask(
maskedImage=interpExp.getMaskedImage(),
fwhm=self.config.brighterFatterFwhmForInterpolation,
growSaturatedFootprints=self.config.growSaturationFootprintSize,
maskNameList=list(self.config.brighterFatterMaskListToInterpolate)
)
isrFunctions.interpolateFromMask(
maskedImage=interpExp.getMaskedImage(),
fwhm=self.config.brighterFatterFwhmForInterpolation,
growSaturatedFootprints=self.config.growSaturationFootprintSize,
maskNameList=list(self.config.brighterFatterMaskListToInterpolate)
)
bfExp = interpExp.clone()
self.log.info("Applying brighter-fatter correction using kernel type %s / gains %s.",
type(bfKernel), type(bfGains))
self.log.info("Applying brighter-fatter correction.")
bfResults = isrFunctions.brighterFatterCorrection(bfExp, bfKernel,
self.config.brighterFatterMaxIter,
self.config.brighterFatterThreshold,
self.config.brighterFatterApplyGain,
brighterFatterApplyGain,
bfGains)
if bfResults[1] == self.config.brighterFatterMaxIter:
self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
Expand Down Expand Up @@ -1655,12 +1640,30 @@ def run(self, ccdExposure, *, dnlLUT=None, bias=None, deferredChargeCalib=None,
isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())

# Brighter/Fatter
# FIXME: watch out (again) for gain units.
# Units: electron
if self.config.doBrighterFatter:
self.log.info("Applying Bright-Fatter kernels.")
self.log.info("Applying brighter-fatter correction.")

bfKernelOut, bfGains = self.getBrighterFatterKernel(detector, bfKernel)
ccdExposure = self.applyBrighterFatterCorrection(ccdExposure, flat, dark, bfKernelOut, bfGains)

# This needs to be done in electrons; applyBrighterFatterCorrection()

Check failure on line 1649 in python/lsst/ip/isr/isrTaskLSST.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

W505 

doc line too long (81 > 79 characters)
# will convert the image if necessary.
if exposureMetadata["LSST ISR UNITS"] == "electron":
brighterFatterApplyGain = False
else:
brighterFatterApplyGain = True
exposureMetadata["LSST ISR UNITS"] == "electron"

if brighterFatterApplyGain and (ptc is not None) and (not bfGains == gains):
# The supplied ptc should be the same as the ptc used to generate the bfKernel, in which case

Check failure on line 1658 in python/lsst/ip/isr/isrTaskLSST.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

W505 

doc line too long (109 > 79 characters)
# they will have the same stored gains. If not, there is a mismatch in the calibrations being

Check failure on line 1659 in python/lsst/ip/isr/isrTaskLSST.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

W505 

doc line too long (109 > 79 characters)
# used. This should not be always be a fatal error, but ideally, everything should to be

Check failure on line 1660 in python/lsst/ip/isr/isrTaskLSST.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

W505 

doc line too long (105 > 79 characters)

Check failure on line 1660 in python/lsst/ip/isr/isrTaskLSST.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

W291 

trailing whitespace
# consistent.
self.log.warning("Need to apply gain for brighter-fatter, but the stored gains in the kernel"
"are not the same as the gains stored in the PTC. Using the kernel gains.")

ccdExposure = self.applyBrighterFatterCorrection(ccdExposure, flat, dark, bfKernelOut,

Check failure on line 1665 in python/lsst/ip/isr/isrTaskLSST.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

W291 

trailing whitespace
brighterFatterApplyGain, bfGains)

# Variance plane creation
# Units: electron
Expand Down
Loading

0 comments on commit 2dfb481

Please sign in to comment.