LSST Applications g013ef56533+7c9321ec0f,g042eb84c57+c6cfa41bc3,g199a45376c+0ba108daf9,g1fd858c14a+fcad0d0313,g210f2d0738+c0f94c6586,g262e1987ae+a7e710680e,g29ae962dfc+fb55f2edb0,g2ac17093b6+61d6563b1e,g2b1d02342f+df6f932764,g2cef7863aa+aef1011c0b,g2f7ad74990+c0f94c6586,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53cf87ae69,g47891489e3+4316d04fff,g511e8cfd20+baa56acf6c,g53246c7159+8c5ae1fdc5,g54cd7ddccb+fd7ad03fde,g64539dfbff+c0f94c6586,g67b6fd64d1+4316d04fff,g67fd3c3899+c0f94c6586,g6985122a63+4316d04fff,g74acd417e5+ca833bee28,g786e29fd12+668abc6043,g81db2e9a8d+b2ec8e584f,g87389fa792+8856018cbb,g89139ef638+4316d04fff,g8d7436a09f+0a24083b20,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+11eb8fd5b8,gbf99507273+8c5ae1fdc5,gcdda8b9158+e4c84c9d5c,gce8aa8abaa+8c5ae1fdc5,gd7ef33dd92+4316d04fff,gdab6d2f7ff+ca833bee28,ge410e46f29+4316d04fff,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.40
LSST Data Management Base Package
Loading...
Searching...
No Matches
linearize.py
Go to the documentation of this file.
2# LSST Data Management System
3# Copyright 2016 AURA/LSST.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <http://www.lsstcorp.org/LegalNotices/>.
21#
22
23__all__ = [
24 "Linearizer",
25 "LinearizeBase",
26 "LinearizeLookupTable",
27 "LinearizeSquared",
28 "LinearizeProportional",
29 "LinearizePolynomial",
30 "LinearizeSpline",
31 "LinearizeDoubleSpline",
32 "LinearizeNone",
33]
34
35import abc
36import numpy as np
37from scipy.interpolate import Akima1DInterpolator
38
39from astropy.table import Table
40
41from lsst.pipe.base import Struct
42from lsst.geom import Box2I, Point2I, Extent2I
43from .applyLookupTable import applyLookupTable
44from .calibType import IsrCalib
45from .isrFunctions import isTrimmedImage
46
47
49 """Parameter set for linearization.
50
51 These parameters are included in `lsst.afw.cameraGeom.Amplifier`, but
52 should be accessible externally to allow for testing.
53
54 Parameters
55 ----------
56 table : `numpy.array`, optional
57 Lookup table; a 2-dimensional array of floats:
58
59 - one row for each row index (value of coef[0] in the amplifier)
60 - one column for each image value
61
62 To avoid copying the table the last index should vary fastest
63 (numpy default "C" order)
64 detector : `lsst.afw.cameraGeom.Detector`, optional
65 Detector object. Passed to self.fromDetector() on init.
66 log : `logging.Logger`, optional
67 Logger to handle messages.
68 kwargs : `dict`, optional
69 Other keyword arguments to pass to the parent init.
70
71 Raises
72 ------
73 RuntimeError
74 Raised if the supplied table is not 2D, or if the table has fewer
75 columns than rows (indicating that the indices are swapped).
76
77 Notes
78 -----
79 The linearizer attributes stored are:
80
81 hasLinearity : `bool`
82 Whether a linearity correction is defined for this detector.
83 override : `bool`
84 Whether the detector parameters should be overridden.
85 ampNames : `list` [`str`]
86 List of amplifier names to correct.
87 linearityCoeffs : `dict` [`str`, `np.ndarray`]
88 Coefficients to use in correction. Indexed by amplifier
89 names. The format of the array depends on the type of
90 correction to apply.
91 linearityType : `dict` [`str`, `str`]
92 Type of correction to use, indexed by amplifier names.
93 linearityBBox : `dict` [`str`, `lsst.geom.Box2I`]
94 Bounding box the correction is valid over, indexed by
95 amplifier names.
96 fitParams : `dict` [`str`, `np.ndarray`], optional
97 Linearity fit parameters used to construct the correction
98 coefficients, indexed as above.
99 fitParamsErr : `dict` [`str`, `np.ndarray`], optional
100 Uncertainty values of the linearity fit parameters used to
101 construct the correction coefficients, indexed as above.
102 fitChiSq : `dict` [`str`, `float`], optional
103 Chi-squared value of the linearity fit, indexed as above.
104 fitResiduals : `dict` [`str`, `np.ndarray`], optional
105 Residuals of the fit, indexed as above. Used for
106 calculating photodiode corrections
107 fitResidualsSigmaMad : `dict` [`str`, `float`], optional
108 Robust median-absolute-deviation of fit residuals, scaled
109 by the signal level.
110 fitResidualsUnmasked : `dict` [`str`, `np.ndarray`], optional
111 Same as fitResiduals, but all outliers are included and
112 not masked as nans.
113 fitResidualsModel : `dict` [`str`, `np.ndarray`], optional
114 The model count level for each of the fitResiduals.
115 linearFit : The linear fit to the low flux region of the curve.
116 [intercept, slope].
117 tableData : `np.ndarray`, optional
118 Lookup table data for the linearity correction.
119 inputAbscissa : `dict` [`str`, `np.ndarray`], optional
120 Input abscissa used to construct linearizer (usually photodiode
121 or exposure time).
122 inputOrdinate : `dict` [`str`, `np.ndarray`], optional
123 Input ordinate used to construct linearizer (raw mean counts).
124 inputMask : `dict` [`str`, `np.ndarray`], optional
125 Input mask used for the fitting.
126 inputGroupingIndex : `dict` [`str`, `np.ndarray`], optional
127 Input grouping index used for fitting.
128 inputNormalization : `dict` [`str`, `np.ndarray`], optional
129 Input normalization that was applied to the abscissa for
130 each pair from the PTC used for the linearization fit.
131 absoluteReferenceAmplifier : `str`, optional
132 Amplifier used for the reference for absolute linearization
133 (DoubleSpline) mode.
134
135 Version 1.4 adds ``linearityTurnoff`` and ``linearityMaxSignal``.
136 Version 1.5 adds ``fitResidualsUnmasked``, ``inputAbscissa``,
137 ``inputOrdinate``, ``inputMask``, ``inputGroupingIndex``,
138 ``fitResidualsModel``, and ``inputNormalization``.
139 Version 1.6 adds ``absoluteReferenceAmplifier``.
140 """
141 _OBSTYPE = "LINEARIZER"
142 _SCHEMA = 'Gen3 Linearizer'
143 _VERSION = 1.6
144
145 def __init__(self, table=None, **kwargs):
146 self.hasLinearity = False
147 self.override = False
148
149 self.ampNames = list()
150 self.linearityCoeffs = dict()
151 self.linearityType = dict()
152 self.linearityBBox = dict()
153 self.inputAbscissa = dict()
154 self.inputOrdinate = dict()
155 self.inputMask = dict()
156 self.inputGroupingIndex = dict()
157 self.inputNormalization = dict()
158 self.fitParams = dict()
159 self.fitParamsErr = dict()
160 self.fitChiSq = dict()
161 self.fitResiduals = dict()
164 self.fitResidualsModel = dict()
165 self.linearFit = dict()
166 self.linearityTurnoff = dict()
167 self.linearityMaxSignal = dict()
169 self.tableData = None
170 if table is not None:
171 if len(table.shape) != 2:
172 raise RuntimeError("table shape = %s; must have two dimensions" % (table.shape,))
173 if table.shape[1] < table.shape[0]:
174 raise RuntimeError("table shape = %s; indices are switched" % (table.shape,))
175 self.tableData = np.array(table, order="C")
176
177 # The linearizer is always natively in adu because it
178 # is computed prior to computing gains.
179 self.linearityUnits = 'adu'
180
181 super().__init__(**kwargs)
182 self.requiredAttributes.update(['hasLinearity', 'override',
183 'ampNames',
184 'linearityCoeffs', 'linearityType', 'linearityBBox',
185 'fitParams', 'fitParamsErr', 'fitChiSq',
186 'fitResiduals', 'fitResidualsSigmaMad', 'linearFit', 'tableData',
187 'units', 'linearityTurnoff', 'linearityMaxSignal',
188 'fitResidualsUnmasked', 'inputAbscissa', 'inputOrdinate',
189 'inputMask', 'inputGroupingIndex', 'fitResidualsModel',
190 'inputNormalization', 'absoluteReferenceAmplifier'])
191
192 def updateMetadata(self, setDate=False, **kwargs):
193 """Update metadata keywords with new values.
194
195 This calls the base class's method after ensuring the required
196 calibration keywords will be saved.
197
198 Parameters
199 ----------
200 setDate : `bool`, optional
201 Update the CALIBDATE fields in the metadata to the current
202 time. Defaults to False.
203 kwargs :
204 Other keyword parameters to set in the metadata.
205 """
206 kwargs['HAS_LINEARITY'] = self.hasLinearity
207 kwargs['OVERRIDE'] = self.override
208 kwargs['HAS_TABLE'] = self.tableData is not None
209 kwargs['LINEARITY_UNITS'] = self.linearityUnits
210
211 super().updateMetadata(setDate=setDate, **kwargs)
212
213 def fromDetector(self, detector):
214 """Read linearity parameters from a detector.
215
216 Parameters
217 ----------
218 detector : `lsst.afw.cameraGeom.detector`
219 Input detector with parameters to use.
220
221 Returns
222 -------
223 calib : `lsst.ip.isr.Linearizer`
224 The calibration constructed from the detector.
225 """
226 self._detectorName = detector.getName()
227 self._detectorSerial = detector.getSerial()
228 self._detectorId = detector.getId()
229 self.hasLinearity = True
230
231 # Do not translate Threshold, Maximum, Units.
232 for amp in detector.getAmplifiers():
233 ampName = amp.getName()
234 self.ampNames.append(ampName)
235 self.linearityType[ampName] = amp.getLinearityType()
236 self.linearityCoeffs[ampName] = amp.getLinearityCoeffs()
237 self.linearityBBox[ampName] = amp.getBBox()
238
239 # Detector linearizers (legacy) are assumed to be adu units.
240 self.linearityUnits = 'adu'
241
242 return self
243
244 @classmethod
245 def fromDict(cls, dictionary):
246 """Construct a calibration from a dictionary of properties
247
248 Parameters
249 ----------
250 dictionary : `dict`
251 Dictionary of properties
252
253 Returns
254 -------
255 calib : `lsst.ip.isr.Linearity`
256 Constructed calibration.
257
258 Raises
259 ------
260 RuntimeError
261 Raised if the supplied dictionary is for a different
262 calibration.
263 """
264
265 calib = cls()
266
267 if calib._OBSTYPE != dictionary['metadata']['OBSTYPE']:
268 raise RuntimeError(f"Incorrect linearity supplied. Expected {calib._OBSTYPE}, "
269 f"found {dictionary['metadata']['OBSTYPE']}")
270
271 calib.setMetadata(dictionary['metadata'])
272
273 calib.hasLinearity = dictionary.get('hasLinearity',
274 dictionary['metadata'].get('HAS_LINEARITY', False))
275 calib.override = dictionary.get('override', True)
276
277 # Old linearizers which do not have linearityUnits are
278 # assumed to be adu because that's all that has been
279 # supported.
280 calib.linearityUnits = dictionary.get('linearityUnits', 'adu')
281
282 if calib.hasLinearity:
283 for ampName in dictionary['amplifiers']:
284 amp = dictionary['amplifiers'][ampName]
285 calib.ampNames.append(ampName)
286 calib.linearityCoeffs[ampName] = np.array(amp.get('linearityCoeffs', [0.0]))
287 calib.linearityType[ampName] = amp.get('linearityType', 'None')
288 calib.linearityBBox[ampName] = amp.get('linearityBBox', None)
289
290 calib.inputAbscissa[ampName] = np.array(amp.get('inputAbscissa', [0.0]))
291 calib.inputOrdinate[ampName] = np.array(amp.get('inputOrdinate', [0.0]))
292 calib.inputMask[ampName] = np.array(amp.get('inputMask', [False]))
293 calib.inputGroupingIndex[ampName] = np.array(amp.get('inputGroupingIndex', [0.0]))
294 calib.inputNormalization[ampName] = np.array(amp.get('inputNormalization', [1.0]))
295
296 calib.fitParams[ampName] = np.array(amp.get('fitParams', [0.0]))
297 calib.fitParamsErr[ampName] = np.array(amp.get('fitParamsErr', [0.0]))
298 calib.fitChiSq[ampName] = amp.get('fitChiSq', np.nan)
299 calib.fitResiduals[ampName] = np.array(amp.get('fitResiduals', [0.0]))
300 calib.fitResidualsSigmaMad[ampName] = np.array(amp.get('fitResidualsSigmaMad', np.nan))
301 calib.fitResidualsUnmasked[ampName] = np.array(amp.get('fitResidualsUnmasked', [0.0]))
302 calib.fitResidualsModel[ampName] = np.array(amp.get('fitResidualsModel', [0.0]))
303 calib.linearFit[ampName] = np.array(amp.get('linearFit', [0.0]))
304
305 calib.linearityTurnoff[ampName] = np.array(amp.get('linearityTurnoff', np.nan))
306 calib.linearityMaxSignal[ampName] = np.array(amp.get('linearityMaxSignal', np.nan))
307
308 calib.tableData = dictionary.get('tableData', None)
309 if calib.tableData:
310 calib.tableData = np.array(calib.tableData)
311
312 calib.absoluteReferenceAmplifier = dictionary.get(
313 'absoluteReferenceAmplifier', calib.absoluteReferenceAmplifier,
314 )
315
316 return calib
317
318 def toDict(self):
319 """Return linearity parameters as a dict.
320
321 Returns
322 -------
323 outDict : `dict`:
324 """
325 self.updateMetadata()
326
327 outDict = {'metadata': self.getMetadata(),
328 'detectorName': self._detectorName,
329 'detectorSerial': self._detectorSerial,
330 'detectorId': self._detectorId,
331 'hasTable': self.tableData is not None,
332 'amplifiers': dict(),
333 'linearityUnits': self.linearityUnits,
334 }
335 for ampName in self.linearityType:
336 outDict['amplifiers'][ampName] = {
337 'linearityType': self.linearityType[ampName],
338 'linearityCoeffs': self.linearityCoeffs[ampName].tolist(),
339 'linearityBBox': self.linearityBBox[ampName],
340 'inputAbscissa': self.inputAbscissa[ampName].tolist(),
341 'inputOrdinate': self.inputOrdinate[ampName].tolist(),
342 'inputMask': self.inputMask[ampName].tolist(),
343 'inputGroupingIndex': self.inputGroupingIndex[ampName].tolist(),
344 'inputNormalization': self.inputNormalization[ampName].tolist(),
345 'fitParams': self.fitParams[ampName].tolist(),
346 'fitParamsErr': self.fitParamsErr[ampName].tolist(),
347 'fitChiSq': self.fitChiSq[ampName],
348 'fitResiduals': self.fitResiduals[ampName].tolist(),
349 'fitResidualsSigmaMad': self.fitResidualsSigmaMad[ampName],
350 'fitResidualsUnmasked': self.fitResidualsUnmasked[ampName].tolist(),
351 'fitResidualsModel': self.fitResidualsModel[ampName].tolist(),
352 'linearFit': self.linearFit[ampName].tolist(),
353 'linearityTurnoff': self.linearityTurnoff[ampName],
354 'linearityMaxSignal': self.linearityMaxSignal[ampName],
355 }
356 if self.tableData is not None:
357 outDict['tableData'] = self.tableData.tolist()
358
359 outDict['absoluteReferenceAmplifier'] = self.absoluteReferenceAmplifier
360
361 return outDict
362
363 @classmethod
364 def fromTable(cls, tableList):
365 """Read linearity from a FITS file.
366
367 This method uses the `fromDict` method to create the
368 calibration, after constructing an appropriate dictionary from
369 the input tables.
370
371 Parameters
372 ----------
373 tableList : `list` [`astropy.table.Table`]
374 afwTable read from input file name.
375
376 Returns
377 -------
378 linearity : `~lsst.ip.isr.linearize.Linearizer``
379 Linearity parameters.
380
381 Notes
382 -----
383 The method reads a FITS file with 1 or 2 extensions. The metadata is
384 read from the header of extension 1, which must exist. Then the table
385 is loaded, and the ['AMPLIFIER_NAME', 'TYPE', 'COEFFS', 'BBOX_X0',
386 'BBOX_Y0', 'BBOX_DX', 'BBOX_DY'] columns are read and used to set each
387 dictionary by looping over rows.
388 Extension 2 is then attempted to read in the try block (which only
389 exists for lookup tables). It has a column named 'LOOKUP_VALUES' that
390 contains a vector of the lookup entries in each row.
391 """
392 coeffTable = tableList[0]
393
394 metadata = coeffTable.meta
395 inDict = dict()
396 inDict['metadata'] = metadata
397 inDict['hasLinearity'] = metadata.get('HAS_LINEARITY', False)
398 inDict['amplifiers'] = dict()
399 inDict['linearityUnits'] = metadata.get('LINEARITY_UNITS', 'adu')
400
401 for record in coeffTable:
402 ampName = record['AMPLIFIER_NAME']
403
404 inputAbscissa = record['INP_ABSCISSA'] if 'INP_ABSCISSA' in record.columns else np.array([0.0])
405 inputOrdinate = record['INP_ORDINATE'] if 'INP_ORDINATE' in record.columns else np.array([0.0])
406 inputMask = record['INP_MASK'] if 'INP_MASK' in record.columns else np.array([False])
407 inputGroupingIndex = record['INP_GROUPING_INDEX'] if 'INP_GROUPING_INDEX' in record.columns \
408 else np.array([0])
409 inputNormalization = record['INP_NORMALIZATION'] if 'INP_NORMALIZATION' in record.columns \
410 else np.array([1.0])
411 fitParams = record['FIT_PARAMS'] if 'FIT_PARAMS' in record.columns else np.array([0.0])
412 fitParamsErr = record['FIT_PARAMS_ERR'] if 'FIT_PARAMS_ERR' in record.columns else np.array([0.0])
413 fitChiSq = record['RED_CHI_SQ'] if 'RED_CHI_SQ' in record.columns else np.nan
414 fitResiduals = record['FIT_RES'] if 'FIT_RES' in record.columns else np.array([0.0])
415 fitResidualsSigmaMad = record['FIT_RES_SIGMAD'] if 'FIT_RES_SIGMAD' in record.columns else np.nan
416 fitResidualsUnmasked = record['FIT_RES_UNMASKED'] \
417 if 'FIT_RES_UNMASKED' in record.columns else np.array([0.0])
418 fitResidualsModel = record['FIT_RES_MODEL'] \
419 if 'FIT_RES_MODEL' in record.columns else np.array([0.0])
420 linearFit = record['LIN_FIT'] if 'LIN_FIT' in record.columns else np.array([0.0])
421
422 linearityTurnoff = record['LINEARITY_TURNOFF'] if 'LINEARITY_TURNOFF' in record.columns \
423 else np.nan
424 linearityMaxSignal = record['LINEARITY_MAX_SIGNAL'] if 'LINEARITY_MAX_SIGNAL' in record.columns \
425 else np.nan
426
427 inDict['amplifiers'][ampName] = {
428 'linearityType': record['TYPE'],
429 'linearityCoeffs': record['COEFFS'],
430 'linearityBBox': Box2I(Point2I(record['BBOX_X0'], record['BBOX_Y0']),
431 Extent2I(record['BBOX_DX'], record['BBOX_DY'])),
432 'inputAbscissa': inputAbscissa,
433 'inputOrdinate': inputOrdinate,
434 'inputMask': inputMask,
435 'inputGroupingIndex': inputGroupingIndex,
436 'inputNormalization': inputNormalization,
437 'fitParams': fitParams,
438 'fitParamsErr': fitParamsErr,
439 'fitChiSq': fitChiSq,
440 'fitResiduals': fitResiduals,
441 'fitResidualsSigmaMad': fitResidualsSigmaMad,
442 'fitResidualsUnmasked': fitResidualsUnmasked,
443 'fitResidualsModel': fitResidualsModel,
444 'linearFit': linearFit,
445 'linearityTurnoff': linearityTurnoff,
446 'linearityMaxSignal': linearityMaxSignal,
447 }
448
449 if len(tableList) > 1:
450 tableData = tableList[1]
451 inDict['tableData'] = [record['LOOKUP_VALUES'] for record in tableData]
452
453 if 'ABS_REF_AMP' in coeffTable.columns:
454 inDict['absoluteReferenceAmplifier'] = str(coeffTable['ABS_REF_AMP'][0])
455 else:
456 inDict['absoluteReferenceAmplifier'] = ''
457
458 return cls().fromDict(inDict)
459
460 def toTable(self):
461 """Construct a list of tables containing the information in this
462 calibration.
463
464 The list of tables should create an identical calibration
465 after being passed to this class's fromTable method.
466
467 Returns
468 -------
469 tableList : `list` [`astropy.table.Table`]
470 List of tables containing the linearity calibration
471 information.
472 """
473
474 tableList = []
475 self.updateMetadata()
476 catalog = Table([{'AMPLIFIER_NAME': ampName,
477 'TYPE': self.linearityType[ampName],
478 'COEFFS': self.linearityCoeffs[ampName],
479 'BBOX_X0': self.linearityBBox[ampName].getMinX(),
480 'BBOX_Y0': self.linearityBBox[ampName].getMinY(),
481 'BBOX_DX': self.linearityBBox[ampName].getWidth(),
482 'BBOX_DY': self.linearityBBox[ampName].getHeight(),
483 'INP_ABSCISSA': self.inputAbscissa[ampName],
484 'INP_ORDINATE': self.inputOrdinate[ampName],
485 'INP_MASK': self.inputMask[ampName],
486 'INP_GROUPING_INDEX': self.inputGroupingIndex[ampName],
487 'INP_NORMALIZATION': self.inputNormalization[ampName],
488 'FIT_PARAMS': self.fitParams[ampName],
489 'FIT_PARAMS_ERR': self.fitParamsErr[ampName],
490 'RED_CHI_SQ': self.fitChiSq[ampName],
491 'FIT_RES': self.fitResiduals[ampName],
492 'FIT_RES_SIGMAD': self.fitResidualsSigmaMad[ampName],
493 'FIT_RES_UNMASKED': self.fitResidualsUnmasked[ampName],
494 'FIT_RES_MODEL': self.fitResidualsModel[ampName],
495 'LIN_FIT': self.linearFit[ampName],
496 'LINEARITY_TURNOFF': self.linearityTurnoff[ampName],
497 'LINEARITY_MAX_SIGNAL': self.linearityMaxSignal[ampName],
498 'ABS_REF_AMP': self.absoluteReferenceAmplifier,
499 } for ampName in self.ampNames])
500 catalog.meta = self.getMetadata().toDict()
501 tableList.append(catalog)
502
503 if self.tableData is not None:
504 catalog = Table([{'LOOKUP_VALUES': value} for value in self.tableData])
505 tableList.append(catalog)
506 return tableList
507
508 def getLinearityTypeByName(self, linearityTypeName):
509 """Determine the linearity class to use from the type name.
510
511 Parameters
512 ----------
513 linearityTypeName : str
514 String name of the linearity type that is needed.
515
516 Returns
517 -------
518 linearityType : `~lsst.ip.isr.linearize.LinearizeBase`
519 The appropriate linearity class to use. If no matching class
520 is found, `None` is returned.
521 """
522 for t in [LinearizeLookupTable,
523 LinearizeSquared,
524 LinearizePolynomial,
525 LinearizeProportional,
526 LinearizeSpline,
527 LinearizeDoubleSpline,
528 LinearizeNone]:
529 if t.LinearityType == linearityTypeName:
530 return t
531 return None
532
533 def validate(self, detector=None, amplifier=None):
534 """Validate linearity for a detector/amplifier.
535
536 Parameters
537 ----------
538 detector : `lsst.afw.cameraGeom.Detector`, optional
539 Detector to validate, along with its amplifiers.
540 amplifier : `lsst.afw.cameraGeom.Amplifier`, optional
541 Single amplifier to validate.
542
543 Raises
544 ------
545 RuntimeError
546 Raised if there is a mismatch in linearity parameters, and
547 the cameraGeom parameters are not being overridden.
548 """
549 amplifiersToCheck = []
550 if detector:
551 if self._detectorName != detector.getName():
552 raise RuntimeError("Detector names don't match: %s != %s" %
553 (self._detectorName, detector.getName()))
554 if int(self._detectorId) != int(detector.getId()):
555 raise RuntimeError("Detector IDs don't match: %s != %s" %
556 (int(self._detectorId), int(detector.getId())))
557 # TODO: DM-38778: This check fails on LATISS due to an
558 # error in the camera configuration.
559 # if self._detectorSerial != detector.getSerial():
560 # raise RuntimeError(
561 # "Detector serial numbers don't match: %s != %s" %
562 # (self._detectorSerial, detector.getSerial()))
563 if len(detector.getAmplifiers()) != len(self.linearityCoeffs.keys()):
564 raise RuntimeError("Detector number of amps = %s does not match saved value %s" %
565 (len(detector.getAmplifiers()),
566 len(self.linearityCoeffs.keys())))
567 amplifiersToCheck.extend(detector.getAmplifiers())
568
569 if amplifier:
570 amplifiersToCheck.extend(amplifier)
571
572 for amp in amplifiersToCheck:
573 ampName = amp.getName()
574 if ampName not in self.linearityCoeffs.keys():
575 raise RuntimeError("Amplifier %s is not in linearity data" %
576 (ampName, ))
577 if amp.getLinearityType() != self.linearityType[ampName]:
578 if self.override:
579 self.log.debug("Overriding amplifier defined linearityType (%s) for %s",
580 self.linearityType[ampName], ampName)
581 else:
582 raise RuntimeError("Amplifier %s type %s does not match saved value %s" %
583 (ampName, amp.getLinearityType(), self.linearityType[ampName]))
584 if (amp.getLinearityCoeffs().shape != self.linearityCoeffs[ampName].shape or not
585 np.allclose(amp.getLinearityCoeffs(), self.linearityCoeffs[ampName], equal_nan=True)):
586 if self.override:
587 self.log.debug("Overriding amplifier defined linearityCoeffs (%s) for %s",
588 self.linearityCoeffs[ampName], ampName)
589 else:
590 raise RuntimeError("Amplifier %s coeffs %s does not match saved value %s" %
591 (ampName, amp.getLinearityCoeffs(), self.linearityCoeffs[ampName]))
592
593 def applyLinearity(self, image, detector=None, log=None, gains=None):
594 """Apply the linearity to an image.
595
596 If the linearity parameters are populated, use those,
597 otherwise use the values from the detector.
598
599 Parameters
600 ----------
601 image : `~lsst.afw.image.image`
602 Image to correct.
603 detector : `~lsst.afw.cameraGeom.detector`, optional
604 Detector to use to determine exposure trimmed state. If
605 supplied, but no other linearity information is provided
606 by the calibration, then the static solution stored in the
607 detector will be used.
608 log : `~logging.Logger`, optional
609 Log object to use for logging.
610 gains : `dict` [`str`, `float`], optional
611 Dictionary of amp name to gain. If this is provided then
612 linearity terms will be converted from adu to electrons.
613 Only used for Spline linearity corrections.
614 """
615 if log is None:
616 log = self.log
617 if detector and not self.hasLinearity:
618 self.fromDetector(detector)
619
620 self.validate(detector)
621
622 # Check if the image is trimmed.
623 isTrimmed = None
624 if detector:
625 isTrimmed = isTrimmedImage(image, detector)
626
627 numAmps = 0
628 numLinearized = 0
629 numOutOfRange = 0
630 for ampName in self.linearityType.keys():
631 linearizer = self.getLinearityTypeByName(self.linearityType[ampName])
632 numAmps += 1
633
634 if gains and self.linearityUnits == 'adu':
635 gainValue = gains[ampName]
636 else:
637 gainValue = 1.0
638
639 if linearizer is not None:
640 match isTrimmed:
641 case True:
642 bbox = detector[ampName].getBBox()
643 case False:
644 bbox = detector[ampName].getRawBBox()
645 case None:
646 bbox = self.linearityBBox[ampName]
647
648 ampView = image.Factory(image, bbox)
649 success, outOfRange = linearizer()(ampView, **{'coeffs': self.linearityCoeffs[ampName],
650 'table': self.tableData,
651 'log': self.log,
652 'gain': gainValue})
653 numOutOfRange += outOfRange
654 if success:
655 numLinearized += 1
656 elif log is not None:
657 log.warning("Amplifier %s did not linearize.",
658 ampName)
659 return Struct(
660 numAmps=numAmps,
661 numLinearized=numLinearized,
662 numOutOfRange=numOutOfRange
663 )
664
665
666class LinearizeBase(metaclass=abc.ABCMeta):
667 """Abstract base class functor for correcting non-linearity.
668
669 Subclasses must define ``__call__`` and set class variable
670 LinearityType to a string that will be used for linearity type in
671 the cameraGeom.Amplifier.linearityType field.
672
673 All linearity corrections should be defined in terms of an
674 additive correction, such that:
675
676 corrected_value = uncorrected_value + f(uncorrected_value)
677 """
678 LinearityType = None # linearity type, a string used for AmpInfoCatalogs
679
680 @abc.abstractmethod
681 def __call__(self, image, **kwargs):
682 """Correct non-linearity.
683
684 Parameters
685 ----------
686 image : `lsst.afw.image.Image`
687 Image to be corrected
688 kwargs : `dict`
689 Dictionary of parameter keywords:
690
691 ``coeffs``
692 Coefficient vector (`list` or `np.ndarray`).
693 ``table``
694 Lookup table data (`np.ndarray`).
695 ``log``
696 Logger to handle messages (`logging.Logger`).
697
698 Returns
699 -------
700 output : `bool`
701 If `True`, a correction was applied successfully.
702
703 Raises
704 ------
705 RuntimeError:
706 Raised if the linearity type listed in the
707 detector does not match the class type.
708 """
709 pass
710
711
712class LinearizeLookupTable(LinearizeBase):
713 """Correct non-linearity with a persisted lookup table.
714
715 The lookup table consists of entries such that given
716 "coefficients" c0, c1:
717
718 for each i,j of image:
719 rowInd = int(c0)
720 colInd = int(c1 + uncorrImage[i,j])
721 corrImage[i,j] = uncorrImage[i,j] + table[rowInd, colInd]
722
723 - c0: row index; used to identify which row of the table to use
724 (typically one per amplifier, though one can have multiple
725 amplifiers use the same table)
726 - c1: column index offset; added to the uncorrected image value
727 before truncation; this supports tables that can handle
728 negative image values; also, if the c1 ends with .5 then
729 the nearest index is used instead of truncating to the
730 next smaller index
731 """
732 LinearityType = "LookupTable"
733
734 def __call__(self, image, **kwargs):
735 """Correct for non-linearity.
736
737 Parameters
738 ----------
739 image : `lsst.afw.image.Image`
740 Image to be corrected
741 kwargs : `dict`
742 Dictionary of parameter keywords:
743
744 ``coeffs``
745 Columnation vector (`list` or `np.ndarray`).
746 ``table``
747 Lookup table data (`np.ndarray`).
748 ``log``
749 Logger to handle messages (`logging.Logger`).
750
751 Returns
752 -------
753 output : `tuple` [`bool`, `int`]
754 If true, a correction was applied successfully. The
755 integer indicates the number of pixels that were
756 uncorrectable by being out of range.
757
758 Raises
759 ------
760 RuntimeError:
761 Raised if the requested row index is out of the table
762 bounds.
763 """
764 numOutOfRange = 0
765
766 rowInd, colIndOffset = kwargs['coeffs'][0:2]
767 table = kwargs['table']
768 log = kwargs['log']
769
770 numTableRows = table.shape[0]
771 rowInd = int(rowInd)
772 if rowInd < 0 or rowInd > numTableRows:
773 raise RuntimeError("LinearizeLookupTable rowInd=%s not in range[0, %s)" %
774 (rowInd, numTableRows))
775 tableRow = np.array(table[rowInd, :], dtype=image.getArray().dtype)
776
777 numOutOfRange += applyLookupTable(image, tableRow, colIndOffset)
778
779 if numOutOfRange > 0 and log is not None:
780 log.warning("%s pixels were out of range of the linearization table",
781 numOutOfRange)
782 if numOutOfRange < image.getArray().size:
783 return True, numOutOfRange
784 else:
785 return False, numOutOfRange
786
787
789 """Correct non-linearity with a polynomial mode.
790
791 .. code-block::
792
793 corrImage = uncorrImage + sum_i c_i uncorrImage^(2 + i)
794
795 where ``c_i`` are the linearity coefficients for each amplifier.
796 Lower order coefficients are not included as they duplicate other
797 calibration parameters:
798
799 ``k0``
800 A coefficient multiplied by ``uncorrImage**0`` is equivalent to
801 bias level. Irrelevant for correcting non-linearity.
802 ``k1``
803 A coefficient multiplied by ``uncorrImage**1`` is proportional
804 to the gain. Not necessary for correcting non-linearity.
805 """
806 LinearityType = "Polynomial"
807
808 def __call__(self, image, **kwargs):
809 """Correct non-linearity.
810
811 Parameters
812 ----------
813 image : `lsst.afw.image.Image`
814 Image to be corrected
815 kwargs : `dict`
816 Dictionary of parameter keywords:
817
818 ``coeffs``
819 Coefficient vector (`list` or `np.ndarray`).
820 If the order of the polynomial is n, this list
821 should have a length of n-1 ("k0" and "k1" are
822 not needed for the correction).
823 ``log``
824 Logger to handle messages (`logging.Logger`).
825
826 Returns
827 -------
828 output : `tuple` [`bool`, `int`]
829 If true, a correction was applied successfully. The
830 integer indicates the number of pixels that were
831 uncorrectable by being out of range.
832 """
833 if not np.any(np.isfinite(kwargs['coeffs'])):
834 return False, 0
835 if not np.any(kwargs['coeffs']):
836 return False, 0
837
838 ampArray = image.getArray()
839 correction = np.zeros_like(ampArray)
840 for order, coeff in enumerate(kwargs['coeffs'], start=2):
841 correction += coeff * np.power(ampArray, order)
842 ampArray += correction
843
844 return True, 0
845
846
848 """Correct non-linearity with a squared model.
849
850 corrImage = uncorrImage + c0*uncorrImage^2
851
852 where c0 is linearity coefficient 0 for each amplifier.
853 """
854 LinearityType = "Squared"
855
856 def __call__(self, image, **kwargs):
857 """Correct for non-linearity.
858
859 Parameters
860 ----------
861 image : `lsst.afw.image.Image`
862 Image to be corrected
863 kwargs : `dict`
864 Dictionary of parameter keywords:
865
866 ``coeffs``
867 Coefficient vector (`list` or `np.ndarray`).
868 ``log``
869 Logger to handle messages (`logging.Logger`).
870
871 Returns
872 -------
873 output : `tuple` [`bool`, `int`]
874 If true, a correction was applied successfully. The
875 integer indicates the number of pixels that were
876 uncorrectable by being out of range.
877 """
878
879 sqCoeff = kwargs['coeffs'][0]
880 if sqCoeff != 0:
881 ampArr = image.getArray()
882 ampArr *= (1 + sqCoeff*ampArr)
883 return True, 0
884 else:
885 return False, 0
886
887
889 """Correct non-linearity with a spline model.
890
891 corrImage = uncorrImage - Spline(coeffs, uncorrImage)
892
893 Notes
894 -----
895
896 The spline fit calculates a correction as a function of the
897 expected linear flux term. Because of this, the correction needs
898 to be subtracted from the observed flux.
899
900 """
901 LinearityType = "Spline"
902
903 def __call__(self, image, **kwargs):
904 """Correct for non-linearity.
905
906 Parameters
907 ----------
908 image : `lsst.afw.image.Image`
909 Image to be corrected
910 kwargs : `dict`
911 Dictionary of parameter keywords:
912
913 ``coeffs``
914 Coefficient vector (`list` or `np.ndarray`).
915 ``log``
916 Logger to handle messages (`logging.Logger`).
917 ``gain``
918 Gain value to apply.
919
920 Returns
921 -------
922 output : `tuple` [`bool`, `int`]
923 If true, a correction was applied successfully. The
924 integer indicates the number of pixels that were
925 uncorrectable by being out of range.
926 """
927 splineCoeff = kwargs['coeffs']
928 gain = kwargs.get('gain', 1.0)
929 centers, values = np.split(splineCoeff, 2)
930 values = values*gain
931 # If the spline is not anchored at zero, remove the offset
932 # found at the lowest flux available, and add an anchor at
933 # flux=0.0 if there's no entry at that point.
934 if values[0] != 0:
935 offset = values[0]
936 values -= offset
937 if centers[0] != 0.0:
938 centers = np.concatenate(([0.0], centers))
939 values = np.concatenate(([0.0], values))
940
941 ampArr = image.array
942
943 if np.any(~np.isfinite(values)):
944 # This cannot be used; turns everything into nans.
945 ampArr[:] = np.nan
946
947 return False, 0
948 else:
949 interp = Akima1DInterpolator(centers, values, method="akima")
950 # Clip to avoid extrapolation and hitting the top end.
951 ampArr -= interp(np.clip(ampArr, centers[0], centers[-1] - 0.1))
952
953 return True, 0
954
955
957 """Correct non-linearity with a spline model.
958
959 corrImage1 = uncorrImage - Spline1(coeffs, uncorrImage)
960 corrImage = corrImage1 - Spline2(coeffs, corrImage1)
961
962 Notes
963 -----
964
965 The spline fit calculates a correction as a function of the
966 expected linear flux term. Because of this, the correction needs
967 to be subtracted from the observed flux.
968
969 """
970 LinearityType = "DoubleSpline"
971
972 def __call__(self, image, **kwargs):
973 """Correct for non-linearity.
974
975 Parameters
976 ----------
977 image : `lsst.afw.image.Image`
978 Image to be corrected
979 kwargs : `dict`
980 Dictionary of parameter keywords:
981
982 ``coeffs``
983 Coefficient vector (`list` or `np.ndarray`).
984 ``log``
985 Logger to handle messages (`logging.Logger`).
986 ``gain``
987 Gain value to apply.
988
989 Returns
990 -------
991 output : `tuple` [`bool`, `int`]
992 If true, a correction was applied successfully. The
993 integer indicates the number of pixels that were
994 uncorrectable by being out of range.
995 """
996 coeffs = kwargs['coeffs']
997 gain = kwargs.get('gain', 1.0)
998
999 # The coeffs have [nNodes1, nNodes2, nodes1, values1, nodes2, values2]
1000
1001 nNodes1 = int(coeffs[0])
1002 nNodes2 = int(coeffs[1])
1003
1004 ampArr = image.array
1005
1006 if nNodes1 > 0:
1007 splineCoeff1 = coeffs[2: 2 + 2*nNodes1]
1008 centers1, values1 = np.split(splineCoeff1, 2)
1009
1010 if np.any(~np.isfinite(values1)):
1011 # Bad linearizer.
1012 ampArr[:] = np.nan
1013 return False, 0
1014
1015 if gain != 1.0:
1016 centers1 = centers1 * gain
1017 values1 = values1 * gain
1018 if nNodes2 > 0:
1019 splineCoeff2 = coeffs[2 + 2*nNodes1: 2 + 2*nNodes1 + 2*nNodes2]
1020 centers2, values2 = np.split(splineCoeff2, 2)
1021
1022 if np.any(~np.isfinite(values2)):
1023 # Bad linearizer.
1024 ampArr[:] = np.nan
1025 return False, 0
1026
1027 if gain != 1.0:
1028 centers2 = centers2 * gain
1029 values2 = values2 * gain
1030
1031 # Note the double-spline will always be anchored at zero.
1032
1033 if nNodes1 > 0:
1034 interp1 = Akima1DInterpolator(centers1, values1, method="akima")
1035 # Clip to avoid extrapolation and hitting the top end.
1036 ampArr -= interp1(np.clip(ampArr, centers1[0], centers1[-1] - 0.01))
1037 if nNodes2 > 0:
1038 interp2 = Akima1DInterpolator(centers2, values2, method="akima")
1039 # Clip to avoid extrapolation and hitting the top end.
1040 ampArr -= interp2(np.clip(ampArr, centers2[0], centers2[-1] - 0.01))
1041
1042 return True, 0
1043
1044
1046 """Do not correct non-linearity.
1047 """
1048 LinearityType = "Proportional"
1049
1050 def __call__(self, image, **kwargs):
1051 """Do not correct for non-linearity.
1052
1053 Parameters
1054 ----------
1055 image : `lsst.afw.image.Image`
1056 Image to be corrected
1057 kwargs : `dict`
1058 Dictionary of parameter keywords:
1059
1060 ``coeffs``
1061 Coefficient vector (`list` or `np.ndarray`).
1062 ``log``
1063 Logger to handle messages (`logging.Logger`).
1064
1065 Returns
1066 -------
1067 output : `tuple` [`bool`, `int`]
1068 If true, a correction was applied successfully. The
1069 integer indicates the number of pixels that were
1070 uncorrectable by being out of range.
1071 """
1072 return True, 0
1073
1074
1076 """Do not correct non-linearity.
1077 """
1078 LinearityType = "None"
1079
1080 def __call__(self, image, **kwargs):
1081 """Do not correct for non-linearity.
1082
1083 Parameters
1084 ----------
1085 image : `lsst.afw.image.Image`
1086 Image to be corrected
1087 kwargs : `dict`
1088 Dictionary of parameter keywords:
1089
1090 ``coeffs``
1091 Coefficient vector (`list` or `np.ndarray`).
1092 ``log``
1093 Logger to handle messages (`logging.Logger`).
1094
1095 Returns
1096 -------
1097 output : `tuple` [`bool`, `int`]
1098 If true, a correction was applied successfully. The
1099 integer indicates the number of pixels that were
1100 uncorrectable by being out of range.
1101 """
1102 return True, 0
An integer coordinate rectangle.
Definition Box.h:55
updateMetadata(self, camera=None, detector=None, filterName=None, setCalibId=False, setCalibInfo=False, setDate=False, **kwargs)
Definition calibType.py:210
__call__(self, image, **kwargs)
Definition linearize.py:681
__call__(self, image, **kwargs)
getLinearityTypeByName(self, linearityTypeName)
Definition linearize.py:508
__init__(self, table=None, **kwargs)
Definition linearize.py:145
updateMetadata(self, setDate=False, **kwargs)
Definition linearize.py:192
applyLinearity(self, image, detector=None, log=None, gains=None)
Definition linearize.py:593
validate(self, detector=None, amplifier=None)
Definition linearize.py:533
isTrimmedImage(image, detector)