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

162 self.fitResidualsSigmaMad = dict() 

163 self.fitResidualsUnmasked = dict() 

164 self.fitResidualsModel = dict() 

165 self.linearFit = dict() 

166 self.linearityTurnoff = dict() 

167 self.linearityMaxSignal = dict() 

168 self.absoluteReferenceAmplifier = "" 

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 

788class LinearizePolynomial(LinearizeBase): 

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 

847class LinearizeSquared(LinearizeBase): 

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]