Coverage for python/lsst/ip/isr/linearize.py: 12%

1# 

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 

48class Linearizer(IsrCalib): 

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 inputGain : `dict` [`str`, `float`] 

88 Dictionary keyed by amplifirer name containing the input 

89 gain from the PTC used to construct this linearizer. 

90 linearityCoeffs : `dict` [`str`, `np.ndarray`] 

91 Coefficients to use in correction. Indexed by amplifier 

92 names. The format of the array depends on the type of 

93 correction to apply. 

94 linearityType : `dict` [`str`, `str`] 

95 Type of correction to use, indexed by amplifier names. 

96 linearityBBox : `dict` [`str`, `lsst.geom.Box2I`] 

97 Bounding box the correction is valid over, indexed by 

98 amplifier names. 

99 fitParams : `dict` [`str`, `np.ndarray`], optional 

100 Linearity fit parameters used to construct the correction 

101 coefficients, indexed as above. 

102 fitParamsErr : `dict` [`str`, `np.ndarray`], optional 

103 Uncertainty values of the linearity fit parameters used to 

104 construct the correction coefficients, indexed as above. 

105 fitChiSq : `dict` [`str`, `float`], optional 

106 Chi-squared value of the linearity fit, indexed as above. 

107 fitResiduals : `dict` [`str`, `np.ndarray`], optional 

108 Residuals of the fit, indexed as above. Used for 

109 calculating photodiode corrections 

110 fitResidualsSigmaMad : `dict` [`str`, `float`], optional 

111 Robust median-absolute-deviation of fit residuals, scaled 

112 by the signal level. 

113 fitResidualsUnmasked : `dict` [`str`, `np.ndarray`], optional 

114 Same as fitResiduals, but all outliers are included and 

115 not masked as nans. 

116 fitResidualsModel : `dict` [`str`, `np.ndarray`], optional 

117 The model count level for each of the fitResiduals. 

118 linearFit : The linear fit to the low flux region of the curve. 

119 [intercept, slope]. 

120 tableData : `np.ndarray`, optional 

121 Lookup table data for the linearity correction. 

122 inputAbscissa : `dict` [`str`, `np.ndarray`], optional 

123 Input abscissa used to construct linearizer (usually photodiode 

124 or exposure time). 

125 inputOrdinate : `dict` [`str`, `np.ndarray`], optional 

126 Input ordinate used to construct linearizer (raw mean counts). 

127 inputMask : `dict` [`str`, `np.ndarray`], optional 

128 Input mask used for the fitting. 

129 inputGroupingIndex : `dict` [`str`, `np.ndarray`], optional 

130 Input grouping index used for fitting. 

131 inputNormalization : `dict` [`str`, `np.ndarray`], optional 

132 Input normalization that was applied to the abscissa for 

133 each pair from the PTC used for the linearization fit. 

134 absoluteReferenceAmplifier : `str`, optional 

135 Amplifier used for the reference for absolute linearization 

136 (DoubleSpline) mode. 

137 

138 Version 1.4 adds ``linearityTurnoff`` and ``linearityMaxSignal``. 

139 Version 1.5 adds ``fitResidualsUnmasked``, ``inputAbscissa``, 

140 ``inputOrdinate``, ``inputMask``, ``inputGroupingIndex``, 

141 ``fitResidualsModel``, and ``inputNormalization``. 

142 Version 1.6 adds ``absoluteReferenceAmplifier``. 

143 Version 1.7 adds ``inputGain``. 

144 """ 

145 _OBSTYPE = "LINEARIZER" 

146 _SCHEMA = 'Gen3 Linearizer' 

147 _VERSION = 1.7 

148 

149 def __init__(self, table=None, **kwargs): 

150 self.hasLinearity = False 

151 self.override = False 

152 

153 self.ampNames = list() 

154 self.inputGain = dict() 

155 self.linearityCoeffs = dict() 

156 self.linearityType = dict() 

157 self.linearityBBox = dict() 

158 self.inputAbscissa = dict() 

159 self.inputOrdinate = dict() 

160 self.inputMask = dict() 

161 self.inputGroupingIndex = dict() 

162 self.inputNormalization = dict() 

163 self.fitParams = dict() 

164 self.fitParamsErr = dict() 

165 self.fitChiSq = dict() 

166 self.fitResiduals = dict() 

167 self.fitResidualsSigmaMad = dict() 

168 self.fitResidualsUnmasked = dict() 

169 self.fitResidualsModel = dict() 

170 self.linearFit = dict() 

171 self.linearityTurnoff = dict() 

172 self.linearityMaxSignal = dict() 

173 self.absoluteReferenceAmplifier = "" 

174 self.tableData = None 

175 if table is not None: 

176 if len(table.shape) != 2: 

177 raise RuntimeError("table shape = %s; must have two dimensions" % (table.shape,)) 

178 if table.shape[1] < table.shape[0]: 

179 raise RuntimeError("table shape = %s; indices are switched" % (table.shape,)) 

180 self.tableData = np.array(table, order="C") 

181 

182 # The linearizer is always natively in adu because it 

183 # is computed prior to computing gains. 

184 self.linearityUnits = 'adu' 

185 

186 super().__init__(**kwargs) 

187 self.requiredAttributes.update(['hasLinearity', 'override', 

188 'ampNames', 'inputGain', 

189 'linearityCoeffs', 'linearityType', 'linearityBBox', 

190 'fitParams', 'fitParamsErr', 'fitChiSq', 

191 'fitResiduals', 'fitResidualsSigmaMad', 'linearFit', 'tableData', 

192 'units', 'linearityTurnoff', 'linearityMaxSignal', 

193 'fitResidualsUnmasked', 'inputAbscissa', 'inputOrdinate', 

194 'inputMask', 'inputGroupingIndex', 'fitResidualsModel', 

195 'inputNormalization', 'absoluteReferenceAmplifier']) 

196 

197 def updateMetadata(self, setDate=False, **kwargs): 

198 """Update metadata keywords with new values. 

199 

200 This calls the base class's method after ensuring the required 

201 calibration keywords will be saved. 

202 

203 Parameters 

204 ---------- 

205 setDate : `bool`, optional 

206 Update the CALIBDATE fields in the metadata to the current 

207 time. Defaults to False. 

208 kwargs : 

209 Other keyword parameters to set in the metadata. 

210 """ 

211 kwargs['HAS_LINEARITY'] = self.hasLinearity 

212 kwargs['OVERRIDE'] = self.override 

213 kwargs['HAS_TABLE'] = self.tableData is not None 

214 kwargs['LINEARITY_UNITS'] = self.linearityUnits 

215 

216 super().updateMetadata(setDate=setDate, **kwargs) 

217 

218 def fromDetector(self, detector): 

219 """Read linearity parameters from a detector. 

220 

221 Parameters 

222 ---------- 

223 detector : `lsst.afw.cameraGeom.detector` 

224 Input detector with parameters to use. 

225 

226 Returns 

227 ------- 

228 calib : `lsst.ip.isr.Linearizer` 

229 The calibration constructed from the detector. 

230 """ 

231 self._detectorName = detector.getName() 

232 self._detectorSerial = detector.getSerial() 

233 self._detectorId = detector.getId() 

234 self.hasLinearity = True 

235 

236 # Do not translate Threshold, Maximum, Units. 

237 for amp in detector.getAmplifiers(): 

238 ampName = amp.getName() 

239 self.ampNames.append(ampName) 

240 self.linearityType[ampName] = amp.getLinearityType() 

241 self.linearityCoeffs[ampName] = amp.getLinearityCoeffs() 

242 self.linearityBBox[ampName] = amp.getBBox() 

243 

244 # Detector linearizers (legacy) are assumed to be adu units. 

245 self.linearityUnits = 'adu' 

246 

247 return self 

248 

249 @classmethod 

250 def fromDict(cls, dictionary): 

251 """Construct a calibration from a dictionary of properties 

252 

253 Parameters 

254 ---------- 

255 dictionary : `dict` 

256 Dictionary of properties 

257 

258 Returns 

259 ------- 

260 calib : `lsst.ip.isr.Linearity` 

261 Constructed calibration. 

262 

263 Raises 

264 ------ 

265 RuntimeError 

266 Raised if the supplied dictionary is for a different 

267 calibration. 

268 """ 

269 

270 calib = cls() 

271 

272 if calib._OBSTYPE != dictionary['metadata']['OBSTYPE']: 

273 raise RuntimeError(f"Incorrect linearity supplied. Expected {calib._OBSTYPE}, " 

274 f"found {dictionary['metadata']['OBSTYPE']}") 

275 

276 calib.setMetadata(dictionary['metadata']) 

277 

278 calib.hasLinearity = dictionary.get('hasLinearity', 

279 dictionary['metadata'].get('HAS_LINEARITY', False)) 

280 calib.override = dictionary.get('override', True) 

281 

282 # Old linearizers which do not have linearityUnits are 

283 # assumed to be adu because that's all that has been 

284 # supported. 

285 calib.linearityUnits = dictionary.get('linearityUnits', 'adu') 

286 

287 if calib.hasLinearity: 

288 for ampName in dictionary['amplifiers']: 

289 amp = dictionary['amplifiers'][ampName] 

290 calib.ampNames.append(ampName) 

291 # If gain = None on isrFunction.gainContext(), 

292 # the default amp geometry gains are used 

293 # in IsrTask* 

294 calib.inputGain[ampName] = amp.get('inputGain', np.nan) 

295 calib.linearityCoeffs[ampName] = np.array(amp.get('linearityCoeffs', [0.0])) 

296 calib.linearityType[ampName] = amp.get('linearityType', 'None') 

297 calib.linearityBBox[ampName] = amp.get('linearityBBox', None) 

298 

299 calib.inputAbscissa[ampName] = np.array(amp.get('inputAbscissa', [0.0])) 

300 calib.inputOrdinate[ampName] = np.array(amp.get('inputOrdinate', [0.0])) 

301 calib.inputMask[ampName] = np.array(amp.get('inputMask', [False])) 

302 calib.inputGroupingIndex[ampName] = np.array(amp.get('inputGroupingIndex', [0.0])) 

303 calib.inputNormalization[ampName] = np.array(amp.get('inputNormalization', [1.0])) 

304 

305 calib.fitParams[ampName] = np.array(amp.get('fitParams', [0.0])) 

306 calib.fitParamsErr[ampName] = np.array(amp.get('fitParamsErr', [0.0])) 

307 calib.fitChiSq[ampName] = amp.get('fitChiSq', np.nan) 

308 calib.fitResiduals[ampName] = np.array(amp.get('fitResiduals', [0.0])) 

309 calib.fitResidualsSigmaMad[ampName] = np.array(amp.get('fitResidualsSigmaMad', np.nan)) 

310 calib.fitResidualsUnmasked[ampName] = np.array(amp.get('fitResidualsUnmasked', [0.0])) 

311 calib.fitResidualsModel[ampName] = np.array(amp.get('fitResidualsModel', [0.0])) 

312 calib.linearFit[ampName] = np.array(amp.get('linearFit', [0.0])) 

313 

314 calib.linearityTurnoff[ampName] = np.array(amp.get('linearityTurnoff', np.nan)) 

315 calib.linearityMaxSignal[ampName] = np.array(amp.get('linearityMaxSignal', np.nan)) 

316 

317 calib.tableData = dictionary.get('tableData', None) 

318 if calib.tableData: 

319 calib.tableData = np.array(calib.tableData) 

320 

321 calib.absoluteReferenceAmplifier = dictionary.get( 

322 'absoluteReferenceAmplifier', calib.absoluteReferenceAmplifier, 

323 ) 

324 

325 return calib 

326 

327 def toDict(self): 

328 """Return linearity parameters as a dict. 

329 

330 Returns 

331 ------- 

332 outDict : `dict`: 

333 """ 

334 self.updateMetadata() 

335 

336 outDict = {'metadata': self.getMetadata(), 

337 'detectorName': self._detectorName, 

338 'detectorSerial': self._detectorSerial, 

339 'detectorId': self._detectorId, 

340 'hasTable': self.tableData is not None, 

341 'amplifiers': dict(), 

342 'linearityUnits': self.linearityUnits, 

343 } 

344 for ampName in self.linearityType: 

345 outDict['amplifiers'][ampName] = { 

346 'inputGain': self.inputGain[ampName], 

347 'linearityType': self.linearityType[ampName], 

348 'linearityCoeffs': self.linearityCoeffs[ampName].tolist(), 

349 'linearityBBox': self.linearityBBox[ampName], 

350 'inputAbscissa': self.inputAbscissa[ampName].tolist(), 

351 'inputOrdinate': self.inputOrdinate[ampName].tolist(), 

352 'inputMask': self.inputMask[ampName].tolist(), 

353 'inputGroupingIndex': self.inputGroupingIndex[ampName].tolist(), 

354 'inputNormalization': self.inputNormalization[ampName].tolist(), 

355 'fitParams': self.fitParams[ampName].tolist(), 

356 'fitParamsErr': self.fitParamsErr[ampName].tolist(), 

357 'fitChiSq': self.fitChiSq[ampName], 

358 'fitResiduals': self.fitResiduals[ampName].tolist(), 

359 'fitResidualsSigmaMad': self.fitResidualsSigmaMad[ampName], 

360 'fitResidualsUnmasked': self.fitResidualsUnmasked[ampName].tolist(), 

361 'fitResidualsModel': self.fitResidualsModel[ampName].tolist(), 

362 'linearFit': self.linearFit[ampName].tolist(), 

363 'linearityTurnoff': self.linearityTurnoff[ampName], 

364 'linearityMaxSignal': self.linearityMaxSignal[ampName], 

365 } 

366 if self.tableData is not None: 

367 outDict['tableData'] = self.tableData.tolist() 

368 

369 outDict['absoluteReferenceAmplifier'] = self.absoluteReferenceAmplifier 

370 

371 return outDict 

372 

373 @classmethod 

374 def fromTable(cls, tableList): 

375 """Read linearity from a FITS file. 

376 

377 This method uses the `fromDict` method to create the 

378 calibration, after constructing an appropriate dictionary from 

379 the input tables. 

380 

381 Parameters 

382 ---------- 

383 tableList : `list` [`astropy.table.Table`] 

384 afwTable read from input file name. 

385 

386 Returns 

387 ------- 

388 linearity : `~lsst.ip.isr.linearize.Linearizer`` 

389 Linearity parameters. 

390 

391 Notes 

392 ----- 

393 The method reads a FITS file with 1 or 2 extensions. The metadata is 

394 read from the header of extension 1, which must exist. Then the table 

395 is loaded, and the ['AMPLIFIER_NAME', 'TYPE', 'COEFFS', 'BBOX_X0', 

396 'BBOX_Y0', 'BBOX_DX', 'BBOX_DY'] columns are read and used to set each 

397 dictionary by looping over rows. 

398 Extension 2 is then attempted to read in the try block (which only 

399 exists for lookup tables). It has a column named 'LOOKUP_VALUES' that 

400 contains a vector of the lookup entries in each row. 

401 """ 

402 coeffTable = tableList[0] 

403 

404 metadata = coeffTable.meta 

405 inDict = dict() 

406 inDict['metadata'] = metadata 

407 inDict['hasLinearity'] = metadata.get('HAS_LINEARITY', False) 

408 inDict['amplifiers'] = dict() 

409 inDict['linearityUnits'] = metadata.get('LINEARITY_UNITS', 'adu') 

410 

411 for record in coeffTable: 

412 ampName = record['AMPLIFIER_NAME'] 

413 inputGain = record['INPUT_GAIN'] if 'INPUT_GAIN' in record.columns else np.nan 

414 inputAbscissa = record['INP_ABSCISSA'] if 'INP_ABSCISSA' in record.columns else np.array([0.0]) 

415 inputOrdinate = record['INP_ORDINATE'] if 'INP_ORDINATE' in record.columns else np.array([0.0]) 

416 inputMask = record['INP_MASK'] if 'INP_MASK' in record.columns else np.array([False]) 

417 inputGroupingIndex = record['INP_GROUPING_INDEX'] if 'INP_GROUPING_INDEX' in record.columns \ 

418 else np.array([0]) 

419 inputNormalization = record['INP_NORMALIZATION'] if 'INP_NORMALIZATION' in record.columns \ 

420 else np.array([1.0]) 

421 fitParams = record['FIT_PARAMS'] if 'FIT_PARAMS' in record.columns else np.array([0.0]) 

422 fitParamsErr = record['FIT_PARAMS_ERR'] if 'FIT_PARAMS_ERR' in record.columns else np.array([0.0]) 

423 fitChiSq = record['RED_CHI_SQ'] if 'RED_CHI_SQ' in record.columns else np.nan 

424 fitResiduals = record['FIT_RES'] if 'FIT_RES' in record.columns else np.array([0.0]) 

425 fitResidualsSigmaMad = record['FIT_RES_SIGMAD'] if 'FIT_RES_SIGMAD' in record.columns else np.nan 

426 fitResidualsUnmasked = record['FIT_RES_UNMASKED'] \ 

427 if 'FIT_RES_UNMASKED' in record.columns else np.array([0.0]) 

428 fitResidualsModel = record['FIT_RES_MODEL'] \ 

429 if 'FIT_RES_MODEL' in record.columns else np.array([0.0]) 

430 linearFit = record['LIN_FIT'] if 'LIN_FIT' in record.columns else np.array([0.0]) 

431 

432 linearityTurnoff = record['LINEARITY_TURNOFF'] if 'LINEARITY_TURNOFF' in record.columns \ 

433 else np.nan 

434 linearityMaxSignal = record['LINEARITY_MAX_SIGNAL'] if 'LINEARITY_MAX_SIGNAL' in record.columns \ 

435 else np.nan 

436 

437 inDict['amplifiers'][ampName] = { 

438 'inputGain': inputGain, 

439 'linearityType': record['TYPE'], 

440 'linearityCoeffs': record['COEFFS'], 

441 'linearityBBox': Box2I(Point2I(record['BBOX_X0'], record['BBOX_Y0']), 

442 Extent2I(record['BBOX_DX'], record['BBOX_DY'])), 

443 'inputAbscissa': inputAbscissa, 

444 'inputOrdinate': inputOrdinate, 

445 'inputMask': inputMask, 

446 'inputGroupingIndex': inputGroupingIndex, 

447 'inputNormalization': inputNormalization, 

448 'fitParams': fitParams, 

449 'fitParamsErr': fitParamsErr, 

450 'fitChiSq': fitChiSq, 

451 'fitResiduals': fitResiduals, 

452 'fitResidualsSigmaMad': fitResidualsSigmaMad, 

453 'fitResidualsUnmasked': fitResidualsUnmasked, 

454 'fitResidualsModel': fitResidualsModel, 

455 'linearFit': linearFit, 

456 'linearityTurnoff': linearityTurnoff, 

457 'linearityMaxSignal': linearityMaxSignal, 

458 } 

459 

460 if len(tableList) > 1: 

461 tableData = tableList[1] 

462 inDict['tableData'] = [record['LOOKUP_VALUES'] for record in tableData] 

463 

464 if 'ABS_REF_AMP' in coeffTable.columns: 

465 inDict['absoluteReferenceAmplifier'] = str(coeffTable['ABS_REF_AMP'][0]) 

466 else: 

467 inDict['absoluteReferenceAmplifier'] = '' 

468 

469 return cls().fromDict(inDict) 

470 

471 def toTable(self): 

472 """Construct a list of tables containing the information in this 

473 calibration. 

474 

475 The list of tables should create an identical calibration 

476 after being passed to this class's fromTable method. 

477 

478 Returns 

479 ------- 

480 tableList : `list` [`astropy.table.Table`] 

481 List of tables containing the linearity calibration 

482 information. 

483 """ 

484 

485 tableList = [] 

486 self.updateMetadata() 

487 catalog = Table([{'AMPLIFIER_NAME': ampName, 

488 'INPUT_GAIN': self.inputGain[ampName], 

489 'TYPE': self.linearityType[ampName], 

490 'COEFFS': self.linearityCoeffs[ampName], 

491 'BBOX_X0': self.linearityBBox[ampName].getMinX(), 

492 'BBOX_Y0': self.linearityBBox[ampName].getMinY(), 

493 'BBOX_DX': self.linearityBBox[ampName].getWidth(), 

494 'BBOX_DY': self.linearityBBox[ampName].getHeight(), 

495 'INP_ABSCISSA': self.inputAbscissa[ampName], 

496 'INP_ORDINATE': self.inputOrdinate[ampName], 

497 'INP_MASK': self.inputMask[ampName], 

498 'INP_GROUPING_INDEX': self.inputGroupingIndex[ampName], 

499 'INP_NORMALIZATION': self.inputNormalization[ampName], 

500 'FIT_PARAMS': self.fitParams[ampName], 

501 'FIT_PARAMS_ERR': self.fitParamsErr[ampName], 

502 'RED_CHI_SQ': self.fitChiSq[ampName], 

503 'FIT_RES': self.fitResiduals[ampName], 

504 'FIT_RES_SIGMAD': self.fitResidualsSigmaMad[ampName], 

505 'FIT_RES_UNMASKED': self.fitResidualsUnmasked[ampName], 

506 'FIT_RES_MODEL': self.fitResidualsModel[ampName], 

507 'LIN_FIT': self.linearFit[ampName], 

508 'LINEARITY_TURNOFF': self.linearityTurnoff[ampName], 

509 'LINEARITY_MAX_SIGNAL': self.linearityMaxSignal[ampName], 

510 'ABS_REF_AMP': self.absoluteReferenceAmplifier, 

511 } for ampName in self.ampNames]) 

512 catalog.meta = self.getMetadata().toDict() 

513 tableList.append(catalog) 

514 

515 if self.tableData is not None: 

516 catalog = Table([{'LOOKUP_VALUES': value} for value in self.tableData]) 

517 tableList.append(catalog) 

518 return tableList 

519 

520 def getLinearityTypeByName(self, linearityTypeName): 

521 """Determine the linearity class to use from the type name. 

522 

523 Parameters 

524 ---------- 

525 linearityTypeName : str 

526 String name of the linearity type that is needed. 

527 

528 Returns 

529 ------- 

530 linearityType : `~lsst.ip.isr.linearize.LinearizeBase` 

531 The appropriate linearity class to use. If no matching class 

532 is found, `None` is returned. 

533 """ 

534 for t in [LinearizeLookupTable, 

535 LinearizeSquared, 

536 LinearizePolynomial, 

537 LinearizeProportional, 

538 LinearizeSpline, 

539 LinearizeDoubleSpline, 

540 LinearizeNone]: 

541 if t.LinearityType == linearityTypeName: 

542 return t 

543 return None 

544 

545 def validate(self, detector=None, amplifier=None): 

546 """Validate linearity for a detector/amplifier. 

547 

548 Parameters 

549 ---------- 

550 detector : `lsst.afw.cameraGeom.Detector`, optional 

551 Detector to validate, along with its amplifiers. 

552 amplifier : `lsst.afw.cameraGeom.Amplifier`, optional 

553 Single amplifier to validate. 

554 

555 Raises 

556 ------ 

557 RuntimeError 

558 Raised if there is a mismatch in linearity parameters, and 

559 the cameraGeom parameters are not being overridden. 

560 """ 

561 amplifiersToCheck = [] 

562 if detector: 

563 if self._detectorName != detector.getName(): 

564 raise RuntimeError("Detector names don't match: %s != %s" % 

565 (self._detectorName, detector.getName())) 

566 if int(self._detectorId) != int(detector.getId()): 

567 raise RuntimeError("Detector IDs don't match: %s != %s" % 

568 (int(self._detectorId), int(detector.getId()))) 

569 # TODO: DM-38778: This check fails on LATISS due to an 

570 # error in the camera configuration. 

571 # if self._detectorSerial != detector.getSerial(): 

572 # raise RuntimeError( 

573 # "Detector serial numbers don't match: %s != %s" % 

574 # (self._detectorSerial, detector.getSerial())) 

575 if len(detector.getAmplifiers()) != len(self.linearityCoeffs.keys()): 

576 raise RuntimeError("Detector number of amps = %s does not match saved value %s" % 

577 (len(detector.getAmplifiers()), 

578 len(self.linearityCoeffs.keys()))) 

579 amplifiersToCheck.extend(detector.getAmplifiers()) 

580 

581 if amplifier: 

582 amplifiersToCheck.extend(amplifier) 

583 

584 for amp in amplifiersToCheck: 

585 ampName = amp.getName() 

586 if ampName not in self.linearityCoeffs.keys(): 

587 raise RuntimeError("Amplifier %s is not in linearity data" % 

588 (ampName, )) 

589 if amp.getLinearityType() != self.linearityType[ampName]: 

590 if self.override: 

591 self.log.debug("Overriding amplifier defined linearityType (%s) for %s", 

592 self.linearityType[ampName], ampName) 

593 else: 

594 raise RuntimeError("Amplifier %s type %s does not match saved value %s" % 

595 (ampName, amp.getLinearityType(), self.linearityType[ampName])) 

596 if (amp.getLinearityCoeffs().shape != self.linearityCoeffs[ampName].shape or not 

597 np.allclose(amp.getLinearityCoeffs(), self.linearityCoeffs[ampName], equal_nan=True)): 

598 if self.override: 

599 self.log.debug("Overriding amplifier defined linearityCoeffs (%s) for %s", 

600 self.linearityCoeffs[ampName], ampName) 

601 else: 

602 raise RuntimeError("Amplifier %s coeffs %s does not match saved value %s" % 

603 (ampName, amp.getLinearityCoeffs(), self.linearityCoeffs[ampName])) 

604 

605 def applyLinearity(self, image, detector=None, log=None, gains=None): 

606 """Apply the linearity to an image. 

607 

608 If the linearity parameters are populated, use those, 

609 otherwise use the values from the detector. 

610 

611 Parameters 

612 ---------- 

613 image : `~lsst.afw.image.image` 

614 Image to correct. 

615 detector : `~lsst.afw.cameraGeom.detector`, optional 

616 Detector to use to determine exposure trimmed state. If 

617 supplied, but no other linearity information is provided 

618 by the calibration, then the static solution stored in the 

619 detector will be used. 

620 log : `~logging.Logger`, optional 

621 Log object to use for logging. 

622 gains : `dict` [`str`, `float`], optional 

623 Dictionary of amp name to gain. If this is provided then 

624 linearity terms will be converted from adu to electrons. 

625 Only used for Spline linearity corrections. 

626 """ 

627 if log is None: 

628 log = self.log 

629 if detector and not self.hasLinearity: 

630 self.fromDetector(detector) 

631 

632 self.validate(detector) 

633 

634 # Check if the image is trimmed. 

635 isTrimmed = None 

636 if detector: 

637 isTrimmed = isTrimmedImage(image, detector) 

638 

639 numAmps = 0 

640 numLinearized = 0 

641 numOutOfRange = 0 

642 for ampName in self.linearityType.keys(): 

643 linearizer = self.getLinearityTypeByName(self.linearityType[ampName]) 

644 numAmps += 1 

645 

646 if gains and self.linearityUnits == 'adu': 

647 gainValue = gains[ampName] 

648 else: 

649 gainValue = 1.0 

650 

651 if linearizer is not None: 

652 match isTrimmed: 

653 case True: 

654 bbox = detector[ampName].getBBox() 

655 case False: 

656 bbox = detector[ampName].getRawBBox() 

657 case None: 

658 bbox = self.linearityBBox[ampName] 

659 

660 ampView = image.Factory(image, bbox) 

661 success, outOfRange = linearizer()(ampView, **{'coeffs': self.linearityCoeffs[ampName], 

662 'table': self.tableData, 

663 'log': self.log, 

664 'gain': gainValue}) 

665 numOutOfRange += outOfRange 

666 if success: 

667 numLinearized += 1 

668 elif log is not None: 

669 log.warning("Amplifier %s did not linearize.", 

670 ampName) 

671 return Struct( 

672 numAmps=numAmps, 

673 numLinearized=numLinearized, 

674 numOutOfRange=numOutOfRange 

675 ) 

676 

677 

678class LinearizeBase(metaclass=abc.ABCMeta): 

679 """Abstract base class functor for correcting non-linearity. 

680 

681 Subclasses must define ``__call__`` and set class variable 

682 LinearityType to a string that will be used for linearity type in 

683 the cameraGeom.Amplifier.linearityType field. 

684 

685 All linearity corrections should be defined in terms of an 

686 additive correction, such that: 

687 

688 corrected_value = uncorrected_value + f(uncorrected_value) 

689 """ 

690 LinearityType = None # linearity type, a string used for AmpInfoCatalogs 

691 

692 @abc.abstractmethod 

693 def __call__(self, image, **kwargs): 

694 """Correct non-linearity. 

695 

696 Parameters 

697 ---------- 

698 image : `lsst.afw.image.Image` 

699 Image to be corrected 

700 kwargs : `dict` 

701 Dictionary of parameter keywords: 

702 

703 ``coeffs`` 

704 Coefficient vector (`list` or `np.ndarray`). 

705 ``table`` 

706 Lookup table data (`np.ndarray`). 

707 ``log`` 

708 Logger to handle messages (`logging.Logger`). 

709 

710 Returns 

711 ------- 

712 output : `bool` 

713 If `True`, a correction was applied successfully. 

714 

715 Raises 

716 ------ 

717 RuntimeError: 

718 Raised if the linearity type listed in the 

719 detector does not match the class type. 

720 """ 

721 pass 

722 

723 

724class LinearizeLookupTable(LinearizeBase): 

725 """Correct non-linearity with a persisted lookup table. 

726 

727 The lookup table consists of entries such that given 

728 "coefficients" c0, c1: 

729 

730 for each i,j of image: 

731 rowInd = int(c0) 

732 colInd = int(c1 + uncorrImage[i,j]) 

733 corrImage[i,j] = uncorrImage[i,j] + table[rowInd, colInd] 

734 

735 - c0: row index; used to identify which row of the table to use 

736 (typically one per amplifier, though one can have multiple 

737 amplifiers use the same table) 

738 - c1: column index offset; added to the uncorrected image value 

739 before truncation; this supports tables that can handle 

740 negative image values; also, if the c1 ends with .5 then 

741 the nearest index is used instead of truncating to the 

742 next smaller index 

743 """ 

744 LinearityType = "LookupTable" 

745 

746 def __call__(self, image, **kwargs): 

747 """Correct for non-linearity. 

748 

749 Parameters 

750 ---------- 

751 image : `lsst.afw.image.Image` 

752 Image to be corrected 

753 kwargs : `dict` 

754 Dictionary of parameter keywords: 

755 

756 ``coeffs`` 

757 Columnation vector (`list` or `np.ndarray`). 

758 ``table`` 

759 Lookup table data (`np.ndarray`). 

760 ``log`` 

761 Logger to handle messages (`logging.Logger`). 

762 

763 Returns 

764 ------- 

765 output : `tuple` [`bool`, `int`] 

766 If true, a correction was applied successfully. The 

767 integer indicates the number of pixels that were 

768 uncorrectable by being out of range. 

769 

770 Raises 

771 ------ 

772 RuntimeError: 

773 Raised if the requested row index is out of the table 

774 bounds. 

775 """ 

776 numOutOfRange = 0 

777 

778 rowInd, colIndOffset = kwargs['coeffs'][0:2] 

779 table = kwargs['table'] 

780 log = kwargs['log'] 

781 

782 numTableRows = table.shape[0] 

783 rowInd = int(rowInd) 

784 if rowInd < 0 or rowInd > numTableRows: 

785 raise RuntimeError("LinearizeLookupTable rowInd=%s not in range[0, %s)" % 

786 (rowInd, numTableRows)) 

787 tableRow = np.array(table[rowInd, :], dtype=image.getArray().dtype) 

788 

789 numOutOfRange += applyLookupTable(image, tableRow, colIndOffset) 

790 

791 if numOutOfRange > 0 and log is not None: 

792 log.warning("%s pixels were out of range of the linearization table", 

793 numOutOfRange) 

794 if numOutOfRange < image.getArray().size: 

795 return True, numOutOfRange 

796 else: 

797 return False, numOutOfRange 

798 

799 

800class LinearizePolynomial(LinearizeBase): 

801 """Correct non-linearity with a polynomial mode. 

802 

803 .. code-block:: 

804 

805 corrImage = uncorrImage + sum_i c_i uncorrImage^(2 + i) 

806 

807 where ``c_i`` are the linearity coefficients for each amplifier. 

808 Lower order coefficients are not included as they duplicate other 

809 calibration parameters: 

810 

811 ``k0`` 

812 A coefficient multiplied by ``uncorrImage**0`` is equivalent to 

813 bias level. Irrelevant for correcting non-linearity. 

814 ``k1`` 

815 A coefficient multiplied by ``uncorrImage**1`` is proportional 

816 to the gain. Not necessary for correcting non-linearity. 

817 """ 

818 LinearityType = "Polynomial" 

819 

820 def __call__(self, image, **kwargs): 

821 """Correct non-linearity. 

822 

823 Parameters 

824 ---------- 

825 image : `lsst.afw.image.Image` 

826 Image to be corrected 

827 kwargs : `dict` 

828 Dictionary of parameter keywords: 

829 

830 ``coeffs`` 

831 Coefficient vector (`list` or `np.ndarray`). 

832 If the order of the polynomial is n, this list 

833 should have a length of n-1 ("k0" and "k1" are 

834 not needed for the correction). 

835 ``log`` 

836 Logger to handle messages (`logging.Logger`). 

837 

838 Returns 

839 ------- 

840 output : `tuple` [`bool`, `int`] 

841 If true, a correction was applied successfully. The 

842 integer indicates the number of pixels that were 

843 uncorrectable by being out of range. 

844 """ 

845 if not np.any(np.isfinite(kwargs['coeffs'])): 

846 return False, 0 

847 if not np.any(kwargs['coeffs']): 

848 return False, 0 

849 

850 ampArray = image.getArray() 

851 correction = np.zeros_like(ampArray) 

852 for order, coeff in enumerate(kwargs['coeffs'], start=2): 

853 correction += coeff * np.power(ampArray, order) 

854 ampArray += correction 

855 

856 return True, 0 

857 

858 

859class LinearizeSquared(LinearizeBase): 

860 """Correct non-linearity with a squared model. 

861 

862 corrImage = uncorrImage + c0*uncorrImage^2 

863 

864 where c0 is linearity coefficient 0 for each amplifier. 

865 """ 

866 LinearityType = "Squared" 

867 

868 def __call__(self, image, **kwargs): 

869 """Correct for non-linearity.