22__all__ = (
'DeferredChargeConfig',
29 'FloatingOutputAmplifier',
30 'DeferredChargeCalib',
36from astropy.table
import Table
41from .isrFunctions
import gainContext
42from .calibType
import IsrCalib
44import scipy.interpolate
as interp
48 """Represents a serial register trap.
53 Size of the charge trap, in electrons.
54 emission_time : `float`
55 Trap emission time constant, in inverse transfers.
57 Serial pixel location of the trap, including the prescan.
59 Type of trap capture to use. Should be one of ``linear``,
60 ``logistic``, or ``spline``.
61 coeffs : `list` [`float`]
62 Coefficients for the capture process. Linear traps need one
63 coefficient, logistic traps need two, and spline based traps
64 need to have an even number of coefficients that can be split
65 into their spline locations and values.
70 Raised if the specified parameters are out of expected range.
73 def __init__(self, size, emission_time, pixel, trap_type, coeffs):
75 raise ValueError(
'Trap size must be greater than or equal to 0.')
78 if emission_time <= 0.0:
79 raise ValueError(
'Emission time must be greater than 0.')
80 if np.isnan(emission_time):
81 raise ValueError(
'Emission time must be real-valued, not NaN')
84 if int(pixel) != pixel:
85 raise ValueError(
'Fraction value for pixel not allowed.')
91 if self.
trap_type not in (
'linear',
'logistic',
'spline'):
92 raise ValueError(
'Unknown trap type: %s', self.
trap_type)
97 centers, values = np.split(np.array(self.
coeffs, dtype=np.float64), 2)
99 values = values[~np.isnan(centers)]
100 centers = centers[~np.isnan(centers)]
101 centers = centers[~np.isnan(values)]
102 values = values[~np.isnan(values)]
107 fill_value=(values[0], values[-1]),
73 def __init__(self, size, emission_time, pixel, trap_type, coeffs):
…
117 if self.
size != other.size:
121 if self.
pixel != other.pixel:
125 if self.
coeffs != other.coeffs:
138 """Initialize trapping arrays for simulated readout.
143 Number of rows to simulate.
145 Number of columns to simulate.
146 prescan_width : `int`
147 Additional transfers due to prescan.
152 Raised if the trap falls outside of the image.
154 if self.
pixel > nx+prescan_width:
155 raise ValueError(
'Trap location {0} must be less than {1}'.format(self.
pixel,
158 self.
_trap_array = np.zeros((ny, nx+prescan_width))
163 """Release charge through exponential decay.
167 released_charge : `float`
173 return released_charge
176 """Perform charge capture using a logistic function.
180 free_charge : `float`
181 Charge available to be trapped.
185 captured_charge : `float`
186 Amount of charge actually trapped.
192 return captured_charge
195 """Trap capture function.
199 pixel_signals : `list` [`float`]
204 captured_charge : `list` [`float`]
205 Amount of charge captured from each pixel.
210 Raised if the trap type is invalid.
214 return np.minimum(self.
size, pixel_signals*scaling)
217 return self.
size/(1.+np.exp(-k*(pixel_signals-f0)))
219 return self.
interp(pixel_signals)
221 raise RuntimeError(f
"Invalid trap capture type: {self.trap_type}.")
225 """Base class for handling model/data fit comparisons.
226 This handles all of the methods needed for the lmfit Minimizer to
232 """Generate a realization of the overscan model, using the specified
233 fit parameters and input signal.
237 params : `lmfit.Parameters`
238 Object containing the model parameters.
239 signal : `np.ndarray`, (nMeasurements)
240 Array of image means.
241 num_transfers : `int`
242 Number of serial transfers that the charge undergoes.
243 start : `int`, optional
244 First overscan column to fit. This number includes the
245 last imaging column, and needs to be adjusted by one when
246 using the overscan bounding box.
247 stop : `int`, optional
248 Last overscan column to fit. This number includes the
249 last imaging column, and needs to be adjusted by one when
250 using the overscan bounding box.
254 results : `np.ndarray`, (nMeasurements, nCols)
257 raise NotImplementedError(
"Subclasses must implement the model calculation.")
260 """Calculate log likelihood of the model.
264 params : `lmfit.Parameters`
265 Object containing the model parameters.
266 signal : `np.ndarray`, (nMeasurements)
267 Array of image means.
268 data : `np.ndarray`, (nMeasurements, nCols)
269 Array of overscan column means from each measurement.
273 Additional position arguments.
275 Additional keyword arguments.
280 The log-likelihood of the observed data given the model
283 model_results = self.
model_results(params, signal, *args, **kwargs)
285 inv_sigma2 = 1.0/(error**2.0)
286 diff = model_results - data
288 return -0.5*(np.sum(inv_sigma2*(diff)**2.))
291 """Calculate negative log likelihood of the model.
295 params : `lmfit.Parameters`
296 Object containing the model parameters.
297 signal : `np.ndarray`, (nMeasurements)
298 Array of image means.
299 data : `np.ndarray`, (nMeasurements, nCols)
300 Array of overscan column means from each measurement.
304 Additional position arguments.
306 Additional keyword arguments.
310 negativelogL : `float`
311 The negative log-likelihood of the observed data given the
314 ll = self.
loglikelihood(params, signal, data, error, *args, **kwargs)
318 def rms_error(self, params, signal, data, error, *args, **kwargs):
319 """Calculate RMS error between model and data.
323 params : `lmfit.Parameters`
324 Object containing the model parameters.
325 signal : `np.ndarray`, (nMeasurements)
326 Array of image means.
327 data : `np.ndarray`, (nMeasurements, nCols)
328 Array of overscan column means from each measurement.
332 Additional position arguments.
334 Additional keyword arguments.
339 The rms error between the model and input data.
341 model_results = self.
model_results(params, signal, *args, **kwargs)
343 diff = model_results - data
344 rms = np.sqrt(np.mean(np.square(diff)))
318 def rms_error(self, params, signal, data, error, *args, **kwargs):
…
348 def difference(self, params, signal, data, error, *args, **kwargs):
349 """Calculate the flattened difference array between model and data.
353 params : `lmfit.Parameters`
354 Object containing the model parameters.
355 signal : `np.ndarray`, (nMeasurements)
356 Array of image means.
357 data : `np.ndarray`, (nMeasurements, nCols)
358 Array of overscan column means from each measurement.
362 Additional position arguments.
364 Additional keyword arguments.
368 difference : `np.ndarray`, (nMeasurements*nCols)
369 The rms error between the model and input data.
371 model_results = self.
model_results(params, signal, *args, **kwargs)
372 diff = (model_results-data).flatten()
348 def difference(self, params, signal, data, error, *args, **kwargs):
…
378 """Simple analytic overscan model."""
382 """Generate a realization of the overscan model, using the specified
383 fit parameters and input signal.
387 params : `lmfit.Parameters`
388 Object containing the model parameters.
389 signal : `np.ndarray`, (nMeasurements)
390 Array of image means.
391 num_transfers : `int`
392 Number of serial transfers that the charge undergoes.
393 start : `int`, optional
394 First overscan column to fit. This number includes the
395 last imaging column, and needs to be adjusted by one when
396 using the overscan bounding box.
397 stop : `int`, optional
398 Last overscan column to fit. This number includes the
399 last imaging column, and needs to be adjusted by one when
400 using the overscan bounding box.
404 res : `np.ndarray`, (nMeasurements, nCols)
407 v = params.valuesdict()
408 v[
'cti'] = 10**v[
'ctiexp']
414 x = np.arange(start, stop+1)
415 res = np.zeros((signal.shape[0], x.shape[0]))
417 for i, s
in enumerate(signal):
425 res[i, :] = (np.minimum(v[
'trapsize'], s*v[
'scaling'])
426 * (np.exp(1/v[
'emissiontime']) - 1.0)
427 * np.exp(-x/v[
'emissiontime'])
428 + s*num_transfers*v[
'cti']**x
429 + v[
'driftscale']*s*np.exp(-x/float(v[
'decaytime'])))
435 """Simulated overscan model."""
438 def model_results(params, signal, num_transfers, amp, start=1, stop=10, trap_type=None):
439 """Generate a realization of the overscan model, using the specified
440 fit parameters and input signal.
444 params : `lmfit.Parameters`
445 Object containing the model parameters.
446 signal : `np.ndarray`, (nMeasurements)
447 Array of image means.
448 num_transfers : `int`
449 Number of serial transfers that the charge undergoes.
450 amp : `lsst.afw.cameraGeom.Amplifier`
451 Amplifier to use for geometry information.
452 start : `int`, optional
453 First overscan column to fit. This number includes the
454 last imaging column, and needs to be adjusted by one when
455 using the overscan bounding box.
456 stop : `int`, optional
457 Last overscan column to fit. This number includes the
458 last imaging column, and needs to be adjusted by one when
459 using the overscan bounding box.
460 trap_type : `str`, optional
461 Type of trap model to use.
465 results : `np.ndarray`, (nMeasurements, nCols)
468 v = params.valuesdict()
478 v[
'cti'] = 10**v[
'ctiexp']
481 if trap_type
is None:
483 elif trap_type ==
'linear':
484 trap =
SerialTrap(v[
'trapsize'], v[
'emissiontime'], 1,
'linear',
486 elif trap_type ==
'logistic':
487 trap =
SerialTrap(v[
'trapsize'], v[
'emissiontime'], 1,
'logistic',
490 raise ValueError(
'Trap type must be linear or logistic or None')
493 imarr = np.zeros((signal.shape[0], amp.getRawDataBBox().getWidth()))
494 ramp =
SegmentSimulator(imarr, amp.getRawSerialPrescanBBox().getWidth(), output_amplifier,
495 cti=v[
'cti'], traps=trap)
496 ramp.ramp_exp(signal)
497 model_results = ramp.readout(serial_overscan_width=amp.getRawSerialOverscanBBox().getWidth(),
498 parallel_overscan_width=0)
500 ncols = amp.getRawSerialPrescanBBox().getWidth() + amp.getRawDataBBox().getWidth()
502 return model_results[:, ncols+start-1:ncols+stop]
438 def model_results(params, signal, num_transfers, amp, start=1, stop=10, trap_type=None):
…
506 """Controls the creation of simulated segment images.
510 imarr : `np.ndarray` (nx, ny)
512 prescan_width : `int`
513 Number of serial prescan columns.
514 output_amplifier : `lsst.cp.pipe.FloatingOutputAmplifier`
515 An object holding some deferred charge parameters.
518 traps : `list` [`lsst.ip.isr.SerialTrap`]
519 Serial traps to simulate.
522 def __init__(self, imarr, prescan_width, output_amplifier, cti=0.0, traps=None):
528 self.
segarr[:, prescan_width:] = imarr
532 if isinstance(cti, np.ndarray):
533 raise ValueError(
"cti must be single value, not an array.")
538 if traps
is not None:
539 if not isinstance(traps, list):
522 def __init__(self, imarr, prescan_width, output_amplifier, cti=0.0, traps=None):
…
545 """Add a trap to the serial register.
549 serial_trap : `lsst.ip.isr.SerialTrap`
554 except AttributeError:
559 """Simulate an image with varying flux illumination per row.
561 This method simulates a segment image where the signal level
562 increases along the horizontal direction, according to the
563 provided list of signal levels.
567 signal_list : `list` [`float`]
568 List of signal levels.
573 Raised if the length of the signal list does not equal the
576 if len(signal_list) != self.
ny:
577 raise ValueError(
"Signal list does not match row count.")
579 ramp = np.tile(signal_list, (self.
nx, 1)).T
582 def readout(self, serial_overscan_width=10, parallel_overscan_width=0):
583 """Simulate serial readout of the segment image.
585 This method performs the serial readout of a segment image
586 given the appropriate SerialRegister object and the properties
587 of the ReadoutAmplifier. Additional arguments can be provided
588 to account for the number of desired overscan transfers. The
589 result is a simulated final segment image, in ADU.
593 serial_overscan_width : `int`, optional
594 Number of serial overscan columns.
595 parallel_overscan_width : `int`, optional
596 Number of parallel overscan rows.
600 result : `np.ndarray` (nx, ny)
601 Simulated image, including serial prescan, serial
602 overscan, and parallel overscan regions. Result in electrons.
605 iy = int(self.
ny + parallel_overscan_width)
608 image = np.random.default_rng().normal(
614 free_charge = copy.deepcopy(self.
segarr)
620 offset = np.zeros(self.
ny)
630 captured_charge = trap.trap_charge(free_charge)
631 free_charge -= captured_charge
634 transferred_charge = free_charge*cte
635 deferred_charge = free_charge*cti
639 transferred_charge[:, 0])
640 image[:iy-parallel_overscan_width, i] += transferred_charge[:, 0] + offset
642 free_charge = np.pad(transferred_charge, ((0, 0), (0, 1)),
643 mode=
'constant')[:, 1:] + deferred_charge
648 released_charge = trap.release_charge()
649 free_charge += released_charge
582 def readout(self, serial_overscan_width=10, parallel_overscan_width=0):
…
655 """Object representing the readout amplifier of a single channel.
660 Gain of the amplifier. Currently not used.
662 Drift scale for the amplifier.
664 Decay time for the bias drift.
665 noise : `float`, optional
666 Amplifier read noise.
667 offset : `float`, optional
671 def __init__(self, gain, scale, decay_time, noise=0.0, offset=0.0):
671 def __init__(self, gain, scale, decay_time, noise=0.0, offset=0.0):
…
680 """Calculate local offset hysteresis.
684 old : `np.ndarray`, (,)
686 signal : `np.ndarray`, (,)
687 Current column measurements.
690 offset : `np.ndarray`
693 new = self.
scale*signal
698 """Update parameter values, if within acceptable values.
703 Drift scale for the amplifier.
705 Decay time for the bias drift.
710 Raised if the input parameters are out of range.
713 raise ValueError(
"Scale must be greater than or equal to 0.")
715 raise ValueError(
"Scale must be real-valued number, not NaN.")
717 if decay_time <= 0.0:
718 raise ValueError(
"Decay time must be greater than 0.")
719 if np.isnan(decay_time):
720 raise ValueError(
"Decay time must be real-valued number, not NaN.")
725 r"""Calibration containing deferred charge/CTI parameters.
727 This includes, parameters from Snyder+2021 and exstimates of
728 the serial and parallel CTI using the extended pixel edge
729 response (EPER) method (also defined in Snyder+2021).
734 Additional parameters to pass to parent constructor.
738 The charge transfer inefficiency attributes stored are:
740 driftScale : `dict` [`str`, `float`]
741 A dictionary, keyed by amplifier name, of the local electronic
742 offset drift scale parameter, A_L in Snyder+2021.
743 decayTime : `dict` [`str`, `float`]
744 A dictionary, keyed by amplifier name, of the local electronic
745 offset decay time, \tau_L in Snyder+2021.
746 globalCti : `dict` [`str`, `float`]
747 A dictionary, keyed by amplifier name, of the mean global CTI
748 paramter, b in Snyder+2021.
749 serialTraps : `dict` [`str`, `lsst.ip.isr.SerialTrap`]
750 A dictionary, keyed by amplifier name, containing a single
751 serial trap for each amplifier.
752 signals : `dict` [`str`, `np.ndarray`]
753 A dictionary, keyed by amplifier name, of the mean signal
754 level for each input measurement.
755 serialEper : `dict` [`str`, `np.ndarray`, `float`]
756 A dictionary, keyed by amplifier name, of the serial EPER
757 estimator of serial CTI, given in a list for each input
759 parallelEper : `dict` [`str`, `np.ndarray`, `float`]
760 A dictionary, keyed by amplifier name, of the parallel
761 EPER estimator of parallel CTI, given in a list for each
763 serialCtiTurnoff : `dict` [`str`, `float`]
764 A dictionary, keyed by amplifier name, of the serial CTI
765 turnoff (unit: electrons).
766 parallelCtiTurnoff : `dict` [`str`, `float`]
767 A dictionary, keyed by amplifier name, of the parallel CTI
768 turnoff (unit: electrons).
769 serialCtiTurnoffSamplingErr : `dict` [`str`, `float`]
770 A dictionary, keyed by amplifier name, of the serial CTI
771 turnoff sampling error (unit: electrons).
772 parallelCtiTurnoffSamplingErr : `dict` [`str`, `float`]
773 A dictionary, keyed by amplifier name, of the parallel CTI
774 turnoff sampling error (unit: electrons).
776 Also, the values contained in this calibration are all derived
777 from and image and overscan in units of electron as these are
778 the most natural units in which to compute deferred charge.
779 However, this means the the user should supply a reliable set
780 of gains when computing the CTI statistics during ISR.
782 Version 1.1 deprecates the USEGAINS attribute and standardizes
783 everything to electron units.
784 Version 1.2 adds the signal, serialEper, parallelEper,
785 serialCtiTurnoff, parallelCtiTurnoff,
786 serialCtiTurnoffSamplingErr, parallelCtiTurnoffSamplingErr
790 _SCHEMA =
'Deferred Charge'
807 if kwargs.pop(
"useGains",
None)
is not None:
808 warnings.warn(
"useGains is deprecated, and will be removed "
809 "after v28.", FutureWarning)
816 self.
requiredAttributes.update([
'driftScale',
'decayTime',
'globalCti',
'serialTraps',
817 'signals',
'serialEper',
'parallelEper',
'serialCtiTurnoff',
818 'parallelCtiTurnoff',
'serialCtiTurnoffSamplingErr',
819 'parallelCtiTurnoffSamplingErr'])
822 """Read metadata parameters from a detector.
826 detector : `lsst.afw.cameraGeom.detector`
827 Input detector with parameters to use.
831 calib : `lsst.ip.isr.Linearizer`
832 The calibration constructed from the detector.
839 """Construct a calibration from a dictionary of properties.
844 Dictionary of properties.
848 calib : `lsst.ip.isr.CalibType`
849 Constructed calibration.
854 Raised if the supplied dictionary is for a different
859 if calib._OBSTYPE != dictionary[
'metadata'][
'OBSTYPE']:
860 raise RuntimeError(f
"Incorrect CTI supplied. Expected {calib._OBSTYPE}, "
861 f
"found {dictionary['metadata']['OBSTYPE']}")
863 calib.setMetadata(dictionary[
'metadata'])
865 calib.driftScale = dictionary[
'driftScale']
866 calib.decayTime = dictionary[
'decayTime']
867 calib.globalCti = dictionary[
'globalCti']
868 calib.serialCtiTurnoff = dictionary[
'serialCtiTurnoff']
869 calib.parallelCtiTurnoff = dictionary[
'parallelCtiTurnoff']
870 calib.serialCtiTurnoffSamplingErr = dictionary[
'serialCtiTurnoffSamplingErr']
871 calib.parallelCtiTurnoffSamplingErr = dictionary[
'parallelCtiTurnoffSamplingErr']
873 allAmpNames = dictionary[
'driftScale'].keys()
878 for ampName
in dictionary[
'serialTraps']:
879 ampTraps = dictionary[
'serialTraps'][ampName]
880 calib.serialTraps[ampName] =
SerialTrap(ampTraps[
'size'], ampTraps[
'emissionTime'],
881 ampTraps[
'pixel'], ampTraps[
'trap_type'],
884 for ampName
in allAmpNames:
885 calib.signals[ampName] = np.array(dictionary[
'signals'][ampName], dtype=np.float64)
886 calib.serialEper[ampName] = np.array(dictionary[
'serialEper'][ampName], dtype=np.float64)
887 calib.parallelEper[ampName] = np.array(dictionary[
'parallelEper'][ampName], dtype=np.float64)
889 calib.updateMetadata()
893 """Return a dictionary containing the calibration properties.
894 The dictionary should be able to be round-tripped through
900 Dictionary of properties.
909 outDict[
'signals'] = self.
signals
917 outDict[
'serialTraps'] = {}
920 'emissionTime': self.
serialTraps[ampName].emission_time,
924 outDict[
'serialTraps'][ampName] = ampTrap
930 """Construct calibration from a list of tables.
932 This method uses the ``fromDict`` method to create the
933 calibration, after constructing an appropriate dictionary from
938 tableList : `list` [`lsst.afw.table.Table`]
939 List of tables to use to construct the CTI
940 calibration. Two tables are expected in this list, the
941 first containing the per-amplifier CTI parameters, and the
942 second containing the parameters for serial traps.
946 calib : `lsst.ip.isr.DeferredChargeCalib`
947 The calibration defined in the tables.
952 Raised if the trap type or trap coefficients are not
955 ampTable = tableList[0]
958 inDict[
'metadata'] = ampTable.meta
959 calibVersion = inDict[
'metadata'][
'CTI_VERSION']
961 amps = ampTable[
'AMPLIFIER']
962 driftScale = ampTable[
'DRIFT_SCALE']
963 decayTime = ampTable[
'DECAY_TIME']
964 globalCti = ampTable[
'GLOBAL_CTI']
966 inDict[
'driftScale'] = {amp: value
for amp, value
in zip(amps, driftScale)}
967 inDict[
'decayTime'] = {amp: value
for amp, value
in zip(amps, decayTime)}
968 inDict[
'globalCti'] = {amp: value
for amp, value
in zip(amps, globalCti)}
971 if calibVersion < 1.1:
975 raise RuntimeError(f
"Using old version of CTI calibration (ver. {calibVersion} < 1.1), "
976 "which is no longer supported.")
977 elif calibVersion < 1.2:
978 inDict[
'signals'] = {amp: np.array([np.nan])
for amp
in amps}
979 inDict[
'serialEper'] = {amp: np.array([np.nan])
for amp
in amps}
980 inDict[
'parallelEper'] = {amp: np.array([np.nan])
for amp
in amps}
981 inDict[
'serialCtiTurnoff'] = {amp: np.nan
for amp
in amps}
982 inDict[
'parallelCtiTurnoff'] = {amp: np.nan
for amp
in amps}
983 inDict[
'serialCtiTurnoffSamplingErr'] = {amp: np.nan
for amp
in amps}
984 inDict[
'parallelCtiTurnoffSamplingErr'] = {amp: np.nan
for amp
in amps}
986 signals = ampTable[
'SIGNALS']
987 serialEper = ampTable[
'SERIAL_EPER']
988 parallelEper = ampTable[
'PARALLEL_EPER']
989 serialCtiTurnoff = ampTable[
'SERIAL_CTI_TURNOFF']
990 parallelCtiTurnoff = ampTable[
'PARALLEL_CTI_TURNOFF']
991 serialCtiTurnoffSamplingErr = ampTable[
'SERIAL_CTI_TURNOFF_SAMPLING_ERR']
992 parallelCtiTurnoffSamplingErr = ampTable[
'PARALLEL_CTI_TURNOFF_SAMPLING_ERR']
993 inDict[
'signals'] = {amp: value
for amp, value
in zip(amps, signals)}
994 inDict[
'serialEper'] = {amp: value
for amp, value
in zip(amps, serialEper)}
995 inDict[
'parallelEper'] = {amp: value
for amp, value
in zip(amps, parallelEper)}
996 inDict[
'serialCtiTurnoff'] = {amp: value
for amp, value
in zip(amps, serialCtiTurnoff)}
997 inDict[
'parallelCtiTurnoff'] = {amp: value
for amp, value
in zip(amps, parallelCtiTurnoff)}
998 inDict[
'serialCtiTurnoffSamplingErr'] = {
999 amp: value
for amp, value
in zip(amps, serialCtiTurnoffSamplingErr)
1001 inDict[
'parallelCtiTurnoffSamplingErr'] = {
1002 amp: value
for amp, value
in zip(amps, parallelCtiTurnoffSamplingErr)
1005 inDict[
'serialTraps'] = {}
1006 trapTable = tableList[1]
1008 amps = trapTable[
'AMPLIFIER']
1009 sizes = trapTable[
'SIZE']
1010 emissionTimes = trapTable[
'EMISSION_TIME']
1011 pixels = trapTable[
'PIXEL']
1012 trap_type = trapTable[
'TYPE']
1013 coeffs = trapTable[
'COEFFS']
1015 for index, amp
in enumerate(amps):
1017 ampTrap[
'size'] = sizes[index]
1018 ampTrap[
'emissionTime'] = emissionTimes[index]
1019 ampTrap[
'pixel'] = pixels[index]
1020 ampTrap[
'trap_type'] = trap_type[index]
1024 inCoeffs = coeffs[index]
1026 nanValues = np.where(np.isnan(inCoeffs))[0]
1027 if nanValues
is not None:
1028 coeffLength = len(inCoeffs)
1029 while breakIndex < coeffLength:
1030 if coeffLength - breakIndex
in nanValues:
1036 outCoeffs = inCoeffs[0: coeffLength - breakIndex]
1038 outCoeffs = inCoeffs
1039 ampTrap[
'coeffs'] = outCoeffs.tolist()
1041 if ampTrap[
'trap_type'] ==
'linear':
1042 if len(ampTrap[
'coeffs']) < 1:
1043 raise ValueError(
"CTI Amplifier %s coefficients for trap has illegal length %d.",
1044 amp, len(ampTrap[
'coeffs']))
1045 elif ampTrap[
'trap_type'] ==
'logistic':
1046 if len(ampTrap[
'coeffs']) < 2:
1047 raise ValueError(
"CTI Amplifier %s coefficients for trap has illegal length %d.",
1048 amp, len(ampTrap[
'coeffs']))
1049 elif ampTrap[
'trap_type'] ==
'spline':
1050 if len(ampTrap[
'coeffs']) % 2 != 0:
1051 raise ValueError(
"CTI Amplifier %s coefficients for trap has illegal length %d.",
1052 amp, len(ampTrap[
'coeffs']))
1054 raise ValueError(
'Unknown trap type: %s', ampTrap[
'trap_type'])
1056 inDict[
'serialTraps'][amp] = ampTrap
1061 """Construct a list of tables containing the information in this
1064 The list of tables should create an identical calibration
1065 after being passed to this class's fromTable method.
1069 tableList : `list` [`lsst.afw.table.Table`]
1070 List of tables containing the crosstalk calibration
1071 information. Two tables are generated for this list, the
1072 first containing the per-amplifier CTI parameters, and the
1073 second containing the parameters for serial traps.
1085 serialCtiTurnoff = []
1086 parallelCtiTurnoff = []
1087 serialCtiTurnoffSamplingErr = []
1088 parallelCtiTurnoffSamplingErr = []
1095 signals.append(self.
signals[amp])
1100 serialCtiTurnoffSamplingErr.append(
1103 parallelCtiTurnoffSamplingErr.append(
1108 'AMPLIFIER': ampList,
1109 'DRIFT_SCALE': driftScale,
1110 'DECAY_TIME': decayTime,
1111 'GLOBAL_CTI': globalCti,
1113 'SERIAL_EPER': serialEper,
1114 'PARALLEL_EPER': parallelEper,
1115 'SERIAL_CTI_TURNOFF': serialCtiTurnoff,
1116 'PARALLEL_CTI_TURNOFF': parallelCtiTurnoff,
1117 'SERIAL_CTI_TURNOFF_SAMPLING_ERR': serialCtiTurnoffSamplingErr,
1118 'PARALLEL_CTI_TURNOFF_SAMPLING_ERR': parallelCtiTurnoffSamplingErr,
1122 tableList.append(ampTable)
1134 maxCoeffLength = np.maximum(maxCoeffLength, len(trap.coeffs))
1139 sizeList.append(trap.size)
1140 timeList.append(trap.emission_time)
1141 pixelList.append(trap.pixel)
1142 typeList.append(trap.trap_type)
1144 coeffs = trap.coeffs
1145 if len(coeffs) != maxCoeffLength:
1146 coeffs = np.pad(coeffs, (0, maxCoeffLength - len(coeffs)),
1147 constant_values=np.nan).tolist()
1148 coeffList.append(coeffs)
1150 trapTable = Table({
'AMPLIFIER': ampList,
1152 'EMISSION_TIME': timeList,
1155 'COEFFS': coeffList})
1157 tableList.append(trapTable)
1163 """Settings for deferred charge correction.
1167 doc=
"Number of prior pixels to use for local offset correction.",
1172 doc=
"Number of prior pixels to use for trap correction.",
1177 doc=
"If true, scale by the gain.",
1180 deprecated=
"This field is no longer used. Will be removed after v28.",
1184 doc=
"If true, set serial prescan and parallel overscan to zero before correction.",
1190 """Task to correct an exposure for charge transfer inefficiency.
1192 This uses the methods described by Snyder et al. 2021, Journal of
1193 Astronimcal Telescopes, Instruments, and Systems, 7,
1194 048002. doi:10.1117/1.JATIS.7.4.048002 (Snyder+21).
1196 ConfigClass = DeferredChargeConfig
1197 _DefaultName =
'isrDeferredCharge'
1199 def run(self, exposure, ctiCalib, gains=None):
1200 """Correct deferred charge/CTI issues.
1204 exposure : `lsst.afw.image.Exposure`
1205 Exposure to correct the deferred charge on.
1206 ctiCalib : `lsst.ip.isr.DeferredChargeCalib`
1207 Calibration object containing the charge transfer
1209 gains : `dict` [`str`, `float`]
1210 A dictionary, keyed by amplifier name, of the gains to
1211 use. If gains is None, the nominal gains in the amplifier
1216 exposure : `lsst.afw.image.Exposure`
1217 The corrected exposure.
1221 This task will read the exposure metadata and determine if
1222 applying gains if necessary. The correction takes place in
1223 units of electrons. If bootstrapping, the gains used
1224 will just be 1.0. and the input/output units will stay in
1225 adu. If the input image is in adu, the output image will be
1226 in units of electrons. If the input image is in electron,
1227 the output image will be in electron.
1229 image = exposure.getMaskedImage().image
1230 detector = exposure.getDetector()
1233 imageUnits = exposure.getMetadata().get(
"LSST ISR UNITS")
1238 if imageUnits ==
"adu":
1244 if applyGains
and gains
is None:
1245 raise RuntimeError(
"No gains supplied for deferred charge correction.")
1247 with gainContext(exposure, image, apply=applyGains, gains=gains, isTrimmed=
False):
1249 for amp
in detector.getAmplifiers():
1250 ampName = amp.getName()
1252 ampImage = image[amp.getRawBBox()]
1253 if self.config.zeroUnusedPixels:
1256 ampImage[amp.getRawParallelOverscanBBox()].array[:, :] = 0.0
1257 ampImage[amp.getRawSerialPrescanBBox()].array[:, :] = 0.0
1261 ampData = self.
flipData(ampImage.array, amp)
1263 if ctiCalib.driftScale[ampName] > 0.0:
1265 ctiCalib.driftScale[ampName],
1266 ctiCalib.decayTime[ampName],
1267 self.config.nPixelOffsetCorrection)
1269 correctedAmpData = ampData.copy()
1272 ctiCalib.serialTraps[ampName],
1273 ctiCalib.globalCti[ampName],
1274 self.config.nPixelTrapCorrection)
1277 correctedAmpData = self.
flipData(correctedAmpData, amp)
1278 image[amp.getRawBBox()].array[:, :] = correctedAmpData[:, :]
1199 def run(self, exposure, ctiCalib, gains=None):
…
1284 """Flip data array such that readout corner is at lower-left.
1288 ampData : `numpy.ndarray`, (nx, ny)
1290 amp : `lsst.afw.cameraGeom.Amplifier`
1291 Amplifier to get readout corner information.
1295 ampData : `numpy.ndarray`, (nx, ny)
1298 X_FLIP = {ReadoutCorner.LL:
False,
1299 ReadoutCorner.LR:
True,
1300 ReadoutCorner.UL:
False,
1301 ReadoutCorner.UR:
True}
1302 Y_FLIP = {ReadoutCorner.LL:
False,
1303 ReadoutCorner.LR:
False,
1304 ReadoutCorner.UL:
True,
1305 ReadoutCorner.UR:
True}
1307 if X_FLIP[amp.getReadoutCorner()]:
1308 ampData = np.fliplr(ampData)
1309 if Y_FLIP[amp.getReadoutCorner()]:
1310 ampData = np.flipud(ampData)
1316 r"""Remove CTI effects from local offsets.
1318 This implements equation 10 of Snyder+21. For an image with
1319 CTI, s'(m, n), the correction factor is equal to the maximum
1320 value of the set of:
1324 {A_L s'(m, n - j) exp(-j t / \tau_L)}_j=0^jmax
1328 inputArr : `numpy.ndarray`, (nx, ny)
1329 Input image data to correct.
1330 drift_scale : `float`
1331 Drift scale (Snyder+21 A_L value) to use in correction.
1332 decay_time : `float`
1333 Decay time (Snyder+21 \tau_L) of the correction.
1334 num_previous_pixels : `int`, optional
1335 Number of previous pixels to use for correction. As the
1336 CTI has an exponential decay, this essentially truncates
1337 the correction where that decay scales the input charge to
1342 outputArr : `numpy.ndarray`, (nx, ny)
1343 Corrected image data.
1345 r = np.exp(-1/decay_time)
1346 Ny, Nx = inputArr.shape
1349 offset = np.zeros((num_previous_pixels, Ny, Nx))
1350 offset[0, :, :] = drift_scale*np.maximum(0, inputArr)
1353 for n
in range(1, num_previous_pixels):
1354 offset[n, :, n:] = drift_scale*np.maximum(0, inputArr[:, :-n])*(r**n)
1356 Linv = np.amax(offset, axis=0)
1357 outputArr = inputArr - Linv
1363 r"""Apply localized trapping inverse operator to pixel signals.
1365 This implements equation 13 of Snyder+21. For an image with
1366 CTI, s'(m, n), the correction factor is equal to the maximum
1367 value of the set of:
1371 {A_L s'(m, n - j) exp(-j t / \tau_L)}_j=0^jmax
1375 inputArr : `numpy.ndarray`, (nx, ny)
1376 Input image data to correct.
1377 trap : `lsst.ip.isr.SerialTrap`
1378 Serial trap describing the capture and release of charge.
1380 Mean charge transfer inefficiency, b from Snyder+21.
1381 num_previous_pixels : `int`, optional
1382 Number of previous pixels to use for correction.
1386 outputArr : `numpy.ndarray`, (nx, ny)
1387 Corrected image data.
1390 Ny, Nx = inputArr.shape
1392 r = np.exp(-1/trap.emission_time)
1395 trap_occupancy = np.zeros((num_previous_pixels, Ny, Nx))
1396 for n
in range(num_previous_pixels):
1397 trap_occupancy[n, :, n+1:] = trap.capture(np.maximum(0, inputArr))[:, :-(n+1)]*(r**n)
1398 trap_occupancy = np.amax(trap_occupancy, axis=0)
1401 C = trap.capture(np.maximum(0, inputArr)) - trap_occupancy*r
1405 R = np.zeros(inputArr.shape)
1406 R[:, 1:] = trap_occupancy[:, 1:]*(1-r)
1409 outputArr = inputArr - a*T
fromDict(cls, dictionary, **kwargs)
updateMetadata(self, camera=None, detector=None, filterName=None, setCalibId=False, setCalibInfo=False, setDate=False, **kwargs)
dict serialCtiTurnoffSamplingErr
dict parallelCtiTurnoffSamplingErr
fromTable(cls, tableList)
fromDict(cls, dictionary)
fromDetector(self, detector)
local_trap_inverse(inputArr, trap, global_cti=0.0, num_previous_pixels=6)
local_offset_inverse(inputArr, drift_scale, decay_time, num_previous_pixels=15)
run(self, exposure, ctiCalib, gains=None)
__init__(self, gain, scale, decay_time, noise=0.0, offset=0.0)
local_offset(self, old, signal)
update_parameters(self, scale, decay_time)
negative_loglikelihood(self, params, signal, data, error, *args, **kwargs)
loglikelihood(self, params, signal, data, error, *args, **kwargs)
rms_error(self, params, signal, data, error, *args, **kwargs)
difference(self, params, signal, data, error, *args, **kwargs)
model_results(params, signal, num_transfers, start=1, stop=10)
ramp_exp(self, signal_list)
add_trap(self, serial_trap)
readout(self, serial_overscan_width=10, parallel_overscan_width=0)
__init__(self, imarr, prescan_width, output_amplifier, cti=0.0, traps=None)
__init__(self, size, emission_time, pixel, trap_type, coeffs)
capture(self, pixel_signals)
initialize(self, ny, nx, prescan_width)
trap_charge(self, free_charge)
model_results(params, signal, num_transfers, start=1, stop=10)
model_results(params, signal, num_transfers, amp, start=1, stop=10, trap_type=None)
gainContext(exp, image, apply, gains=None, invert=False, isTrimmed=True)