Coverage for python/lsst/ip/isr/electrostaticBrighterFatter.py: 7%

1# This file is part of ip_isr.

3# Developed for the LSST Data Management System.

4# This product includes software developed by the LSST Project

5# (https://www.lsst.org).

6# See the COPYRIGHT file at the top-level directory of this distribution

7# for details of code ownership.

9# This program is free software: you can redistribute it and/or modify

10# it under the terms of the GNU General Public License as published by

11# the Free Software Foundation, either version 3 of the License, or

12# (at your option) any later version.

13#

14# This program is distributed in the hope that it will be useful,

15# but WITHOUT ANY WARRANTY; without even the implied warranty of

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

17# GNU General Public License for more details.

18#

19# You should have received a copy of the GNU General Public License

20# along with this program. If not, see <https://www.gnu.org/licenses/>.

21#

22"""Brighter Fatter Kernel calibration definition."""

25__all__ = [

26 'ElectrostaticBrighterFatterDistortionMatrix',

27 'electrostaticBrighterFatterCorrection',

28]

31import pyfftw

32import numpy as np

33from astropy.table import Table

35from . import IsrCalib

36from .isrFunctions import gainContext

39class ElectrostaticBrighterFatterDistortionMatrix(IsrCalib):

40 """Calibration of brighter-fatter kernels for an instrument.

42 ampKernels are the kernels for each amplifier in a detector, as

43 generated by having ``level == 'AMP'``.

45 detectorKernel is the kernel generated for a detector as a

46 whole, as generated by having ``level == 'DETECTOR'``.

48 makeDetectorKernelFromAmpwiseKernels is a method to generate the

49 kernel for a detector, constructed by averaging together the

50 ampwise kernels in the detector. The existing application code is

51 only defined for kernels with ``level == 'DETECTOR'``, so this method

52 is used if the supplied kernel was built with ``level == 'AMP'``.

54 Parameters

55 ----------

56 camera : `lsst.afw.cameraGeom.Camera`

57 Camera describing detector geometry.

58 level : `str`

59 Level the kernels will be generated for.

60 log : `logging.Logger`, optional

61 Log to write messages to.

62 **kwargs :

63 Parameters to pass to parent constructor.

65 Notes

66 -----

67 Version 1.1 adds the `expIdMask` property, and substitutes

68 `means` and `variances` for `rawMeans` and `rawVariances`

69 from the PTC dataset.

71 inputRange : `int`

72 The size of the input aMatrix shape in each dimension.

73 fitRange : `int`

74 The size of the input aMatrix shape in each dimension that is

75 used for fitting the electrostatic model. Must be less than or

76 equal to inputRange.

77 badAmps : `list`, [`str`]

78 List of bad amplifiers names.

79 gain : `dict`, [`str`,`float`]

80 Dictionary keyed by amp names containing the gains inherited

81 from the inputPTC.

82 aMatrix : `numpy.ndarray`

83 The average aMatrix inherited from the inputPTC's good amplifiers.

84 aMatrixSigma : `numpy.ndarray`

85 The uncertainty matrix used to weight the fit.

86 aMatrixModel : `numpy.ndarray`

87 The modeled aMatrix based on the electrostatic fit parameters.

88 aMatrixSum : `float`

89 The sum of the symmetrized aMatrix.

90 aMatrixModelSum : `float`

91 The sum of the symmetrized aMatrixModel.

92 modelNormalization : `list`

93 A two element array of the multiplicative and additive

94 normalization to the aMatrixModel.

95 usedPixels : `numpy.ndarray`, [`bool`]

96 Array of shape like aMatrix containing the mask indicating which

97 elements of the input aMatrix were used to fit the electrostatic

98 model.

99 fitParamNames : `list`, [`str`]

100 List of all the parameter names in the electrostatic fit.

101 freeFitParamNames : `list`, [`str`]

102 List of the parameter names that were allowed to vary during

103 the electrostatic fit.

104 fitParams : `dict`, [`str`, `float`]

105 Dictionary containing each named parameter and its final

106 fitted value.

107 fitParamErrors : `dict`, [`str`, `float`]

108 Dictionary containing each named parameter and its

109 estimated fitting error.

110 fitChi2 : `float`

111 The computed chi squared between the data and the

112 final model.

113 fitReducedChi2 : `float`

114 The computed reduced chi squared between the data

115 and the final model.

116 fitParamCovMatrix : `numpy.ndarray`

117 The estimated covariance matrix between all fit parameters.

118 ath : `numpy.ndarray`

119 something...

120 athMinusBeta : `numpy.ndarray`

121 something...

122 aN : `numpy.ndarray`

123 Array of shape (fitRange, fitRange) containing the

124 computed `North` component of the pixel boundary shift.

125 aS : `numpy.ndarray`

126 Array of shape (fitRange, fitRange) containing the

127 computed `South` component of the pixel boundary shift.

128 aE : `numpy.ndarray`

129 Array of shape (fitRange, fitRange) containing the

130 computed `East` component of the pixel boundary shift.

131 aW : `numpy.ndarray`

132 Array of shape (fitRange, fitRange) containing the

133 computed `West` component of the pixel boundary shift.

134 aMatrixModelByFilter : `dict` [`str`, `numpy.ndarray`]

135 Dict keyed by physical filter name containing the modeled aMatrix

136 arrays (shape (fitRange, fitRange)) per filter. Only populated

137 when color correction is enabled.

138 aNByFilter : `dict` [`str`, `numpy.ndarray`]

139 Dict keyed by physical filter name containing the North pixel

140 boundary shift arrays (shape (fitRange, fitRange)) per filter.

141 Only populated when color correction is enabled.

142 aSByFilter : `dict` [`str`, `numpy.ndarray`]

143 Dict keyed by physical filter name containing the South pixel

144 boundary shift arrays per filter.

145 aEByFilter : `dict` [`str`, `numpy.ndarray`]

146 Dict keyed by physical filter name containing the East pixel

147 boundary shift arrays per filter.

148 aWByFilter : `dict` [`str`, `numpy.ndarray`]

149 Dict keyed by physical filter name containing the West pixel

150 boundary shift arrays per filter.

151 availableFilters : `list` of `str`

152 List of physical filter names for which per-filter data exists;

153 these are the keys used in the ``*ByFilter`` attributes.

154

155 Version 1.1 adds the `aMatrixModelByFilter`, `athByFilter`,

156 `athMinusBetaByFilter`, `aNByFilter`, `aSByFilter`,

157 `aEByFilter`, and `aWByFilter` attributes.

158 """

159 _OBSTYPE = 'BF_DISTORTION_MATRIX'

160 _SCHEMA = 'ELectrostatic brighter-fatter pixel distortion model'

161 _VERSION = 1.1

162

163 def __init__(self, camera=None, inputRange=1, fitRange=None, availableFilters=list(), **kwargs):

164 """

165 Filename refers to an input tuple that contains the

166 boundary shifts for one electron. This file is produced by an

167 electrostatic fit to data extracted from flat-field statistics,

168 implemented in https://gitlab.in2p3.fr/astier/bfptc/tools/fit_cov.py

169 """

170 self.ampNames = []

171 self.inputRange = inputRange

172 if fitRange is None:

173 fitRange = inputRange

174 self.fitRange = fitRange

175 self.badAmps = list()

176 self.gain = dict()

177 self.aMatrix = np.full((inputRange, inputRange), np.nan)

178 self.aMatrixSigma = np.full((inputRange, inputRange), np.nan)

179 self.aMatrixModel = np.full((fitRange, fitRange), np.nan)

180 self.aMatrixSum = np.nan

181 self.aMatrixModelSum = np.nan

182 self.modelNormalization = [np.nan, np.nan]

183 self.usedPixels = np.full((inputRange, inputRange), np.nan)

184 self.fitParamNames = list()

185 self.freeFitParamNames = list()

186 self.fitParams = dict()

187 self.fitParamErrors = dict()

188 self.fitChi2 = np.nan

189 self.fitReducedChi2 = np.nan

190 self.fitParamCovMatrix = np.array([])

191 self.ath = np.full((fitRange, fitRange), np.nan)

192 self.athMinusBeta = np.full((fitRange, fitRange), np.nan)

193 self.aN = np.full((fitRange, fitRange), np.nan)

194 self.aS = np.full((fitRange, fitRange), np.nan)

195 self.aE = np.full((fitRange, fitRange), np.nan)

196 self.aW = np.full((fitRange, fitRange), np.nan)

197 self.availableFilters = availableFilters

198 self.aMatrixModelByFilter = {k: np.full((fitRange, fitRange), np.nan) for k in availableFilters}

199 self.athByFilter = {k: None for k in availableFilters}

200 self.athMinusBetaByFilter = {k: None for k in availableFilters}

201 self.aNByFilter = {k: np.full((fitRange, fitRange), np.nan) for k in availableFilters}

202 self.aSByFilter = {k: np.full((fitRange, fitRange), np.nan) for k in availableFilters}

203 self.aEByFilter = {k: np.full((fitRange, fitRange), np.nan) for k in availableFilters}

204 self.aWByFilter = {k: np.full((fitRange, fitRange), np.nan) for k in availableFilters}

205

206 super().__init__(**kwargs)

207

208 if camera:

209 self.initFromCamera(camera, detectorId=kwargs.get('detectorId', None))

210

211 self.requiredAttributes.update([

212 'inputRange', 'fitRange', 'badAmps', 'gain', 'aMatrix',

213 'aMatrixSigma', 'aMatrixModel', 'aMatrixSum', 'aMatrixModelSum',

214 'modelNormalization', 'usedPixels', 'fitParamNames',

215 'freeFitParamNames', 'fitParams', 'fitParamErrors', 'fitChi2',

216 'fitReducedChi2', 'fitParamCovMatrix', 'ath', 'athMinusBeta',

217 'aN', 'aS', 'aE', 'aW', 'ampNames', 'availableFilters',

218 'aMatrixModelByFilter', 'athByFilter', 'athMinusBetaByFilter',

219 'aNByFilter', 'aSByFilter', 'aEByFilter', 'aWByFilter',

220 ])

221

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

223 """Update calibration metadata.

224

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

226 calibration keywords will be saved.

227

228 Parameters

229 ----------

230 setDate : `bool`, optional

231 Update the CALIBDATE fields in the metadata to the current

232 time. Defaults to False.

233 kwargs :

234 Other keyword parameters to set in the metadata.

235 """

236 kwargs['INPUT_RANGE'] = self.inputRange

237 kwargs['FIT_RANGE'] = self.fitRange

238

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

240

241 def initFromCamera(self, camera, detectorId=None):

242 """Initialize kernel structure from camera.

243

244 Parameters

245 ----------

246 camera : `lsst.afw.cameraGeom.Camera`

247 Camera to use to define geometry.

248 detectorId : `int`, optional

249 Index of the detector to generate.

250

251 Returns

252 -------

253 calib : `lsst.ip.isr.BrighterFatterKernel`

254 The initialized calibration.

255

256 Raises

257 ------

258 RuntimeError

259 Raised if no detectorId is supplied for a calibration with

260 ``level='AMP'``.

261 """

262 self._instrument = camera.getName()

263

264 if detectorId is not None:

265 detector = camera[detectorId]

266 self._detectorId = detectorId

267 self._detectorName = detector.getName()

268 self._detectorSerial = detector.getSerial()

269

270 for amp in detector:

271 ampName = amp.getName()

272 self.ampNames.append(ampName)

273 else:

274 raise RuntimeError("A detectorId must be supplied if "

275 "camera is supplied.")

276

277 return self

278

279 @classmethod

280 def fromDict(cls, dictionary):

281 """Construct a calibration from a dictionary of properties.

282

283 Parameters

284 ----------

285 dictionary : `dict`

286 Dictionary of properties.

287

288 Returns

289 -------

290 calib : `lsst.ip.isr.BrighterFatterKernel`

291 Constructed calibration.

292

293 Raises

294 ------

295 RuntimeError

296 Raised if the supplied dictionary is for a different

297 calibration.

298 Raised if the version of the supplied dictionary is 1.0.

299 """

300 calib = cls()

301

302 if calib._OBSTYPE != (found := dictionary['metadata']['OBSTYPE']):

303 raise RuntimeError(f"Incorrect brighter-fatter calibration supplied. Expected {calib._OBSTYPE}, "

304 f"found {found}")

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

306 calib.calibInfoFromDict(dictionary)

307

308 calib.ampNames = dictionary['ampNames']

309 inputRange = dictionary['inputRange']

310 fitRange = dictionary['fitRange']

311 calib.inputRange = inputRange

312 calib.fitRange = fitRange

313 calib.gain = dictionary['gain']

314 calib.badAmps = dictionary['badAmps']

315 calib.fitParamNames = dictionary['fitParamNames']

316 calib.freeFitParamNames = dictionary['freeFitParamNames']

317 calib.aMatrix = np.array(

318 dictionary['aMatrix'],

319 dtype=np.float64,

320 ).reshape(inputRange, inputRange)

321 calib.aMatrixSigma = np.array(

322 dictionary['aMatrixSigma'],

323 dtype=np.float64,

324 ).reshape(inputRange, inputRange)

325 calib.aMatrixModel = np.array(

326 dictionary['aMatrixSigma'],

327 dtype=np.float64,

328 ).reshape(fitRange, fitRange)

329 calib.aMatrixSum = float(dictionary['aMatrixSum'])

330 calib.aMatrixModelSum = float(dictionary['aMatrixModelSum'])

331 calib.modelNormalization = dictionary['modelNormalization']

332 calib.usedPixels = np.array(dictionary['usedPixels'])

333 calib.fitParams = dictionary['fitParams']

334 calib.fitParamErrors = dictionary['fitParamErrors']

335 calib.fitChi2 = float(dictionary['fitChi2'])

336 calib.fitReducedChi2 = float(dictionary['fitReducedChi2'])

337 calib.fitParamCovMatrix = np.array(

338 dictionary['fitParamCovMatrix'],

339 dtype=np.float64,

340 ).reshape(len(calib.freeFitParamNames), len(calib.freeFitParamNames))

341 calib.ath = np.array(

342 dictionary['ath'],

343 dtype=np.float64,

344 ).reshape(fitRange, fitRange)

345 calib.athMinusBeta = np.array(

346 dictionary['athMinusBeta'],

347 dtype=np.float64,

348 ).reshape(fitRange, fitRange)

349 calib.aN = np.array(

350 dictionary['aN'],

351 dtype=np.float64,

352 ).reshape(fitRange, fitRange)

353 calib.aS = np.array(

354 dictionary['aS'],

355 dtype=np.float64,

356 ).reshape(fitRange, fitRange)

357 calib.aE = np.array(

358 dictionary['aE'],

359 dtype=np.float64,

360 ).reshape(fitRange, fitRange)

361 calib.aW = np.array(

362 dictionary['aW'],

363 dtype=np.float64,

364 ).reshape(fitRange, fitRange)

365

366 filters = dictionary['availableFilters']

367 calib.availableFilters = dictionary['availableFilters']

368 calib.aMatrixModelByFilter = {

369 k: np.array(v, dtype=np.float64).reshape(fitRange, fitRange)

370 for k, v in zip(filters, dictionary['aMatrixModelByFilter'])

371 }

372 calib.athByFilter = {

373 k: np.array(v, dtype=np.float64).reshape(fitRange, fitRange)

374 for k, v in zip(filters, dictionary['athByFilter'])

375 }

376 calib.athMinusBetaByFilter = {

377 k: np.array(v, dtype=np.float64).reshape(fitRange, fitRange)

378 for k, v in zip(filters, dictionary['athMinusBetaByFilter'])

379 }

380 calib.aNByFilter = {

381 k: np.array(v, dtype=np.float64).reshape(fitRange, fitRange)

382 for k, v in zip(filters, dictionary['aNByFilter'])

383 }

384 calib.aSByFilter = {

385 k: np.array(v, dtype=np.float64).reshape(fitRange, fitRange)

386 for k, v in zip(filters, dictionary['aSByFilter'])

387 }

388 calib.aEByFilter = {

389 k: np.array(v, dtype=np.float64).reshape(fitRange, fitRange)

390 for k, v in zip(filters, dictionary['aEByFilter'])

391 }

392 calib.aWByFilter = {

393 k: np.array(v, dtype=np.float64).reshape(fitRange, fitRange)

394 for k, v in zip(filters, dictionary['aWByFilter'])

395 }

396

397 calib.updateMetadata()

398 return calib

399

400 def toDict(self):

401 """Return a dictionary containing the calibration properties.

402

403 The dictionary should be able to be round-tripped through

404 `fromDict`.

405

406 Returns

407 -------

408 dictionary : `dict`

409 Dictionary of properties.

410 """

411 self.updateMetadata()

412

413 def _dictOfArraysToDictOfLists(dictOfArrays):

414 dictOfLists = {}

415 for key, value in dictOfArrays.items():

416 dictOfLists[key] = value.ravel().tolist()

417

418 return dictOfLists

419

420 outDict = dict()

421 metadata = self.getMetadata()

422 outDict['metadata'] = metadata

423

424 outDict['ampNames'] = self.ampNames

425 outDict['inputRange'] = self.inputRange

426 outDict['fitRange'] = self.fitRange

427 outDict['badAmps'] = self.badAmps

428 outDict['fitParamNames'] = self.fitParamNames

429 outDict['freeFitParamNames'] = self.freeFitParamNames

430 outDict['gain'] = self.gain

431 outDict['aMatrix'] = self.aMatrix.ravel().tolist()

432 outDict['aMatrixSigma'] = self.aMatrixSigma.ravel().tolist()

433 outDict['aMatrixModel'] = self.aMatrixModel.ravel().tolist()

434 outDict['aMatrixSum'] = self.aMatrixSum

435 outDict['aMatrixModelSum'] = self.aMatrixModelSum

436 outDict['aMatrixModel'] = self.aMatrixModel.ravel().tolist()

437 outDict['modelNormalization'] = self.modelNormalization

438 outDict['fitParamCovMatrix'] = self.fitParamCovMatrix.ravel().tolist()

439 outDict['usedPixels'] = self.usedPixels.ravel().tolist()

440 outDict['fitParams'] = self.fitParams

441 outDict['fitParamErrors'] = self.fitParamErrors

442 outDict['fitChi2'] = self.fitChi2

443 outDict['fitReducedChi2'] = self.fitReducedChi2

444 outDict['ath'] = self.ath.ravel().tolist()

445 outDict['athMinusBeta'] = self.athMinusBeta.ravel().tolist()

446 outDict['aN'] = self.aN.ravel().tolist()

447 outDict['aS'] = self.aS.ravel().tolist()

448 outDict['aE'] = self.aE.ravel().tolist()

449 outDict['aW'] = self.aW.ravel().tolist()

450 outDict['availableFilters'] = self.availableFilters

451 outDict['aMatrixModelByFilter'] = _dictOfArraysToDictOfLists(self.aMatrixModelByFilter)

452 outDict['aNByFilter'] = _dictOfArraysToDictOfLists(self.aNByFilter)

453 outDict['aSByFilter'] = _dictOfArraysToDictOfLists(self.aSByFilter)

454 outDict['aEByFilter'] = _dictOfArraysToDictOfLists(self.aEByFilter)

455 outDict['aWByFilter'] = _dictOfArraysToDictOfLists(self.aWByFilter)

456 outDict['athByFilter'] = _dictOfArraysToDictOfLists(self.athByFilter)

457 outDict['athMinusBetaByFilter'] = _dictOfArraysToDictOfLists(self.athMinusBetaByFilter)

458

459 return outDict

460

461 @classmethod

462 def fromTable(cls, tableList):

463 """Construct calibration from a list of tables.

464

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

466 calibration, after constructing an appropriate dictionary from

467 the input tables.

468

469 Parameters

470 ----------

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

472 List of tables to use to construct the brighter-fatter

473 calibration.

474

475 Returns

476 -------

477 calib : `lsst.ip.isr.BrighterFatterKernel`

478 The calibration defined in the tables.

479 """

480 record = tableList[0][0]

481

482 metadata = record.meta

483 calibVersion = metadata['BF_DISTORTION_MATRIX_VERSION']

484

485 # Initialize inDict with all expected keys

486 # and empty values of the corresponding type

487 inDict = dict()

488 inDict['metadata'] = metadata

489 inDict['ampNames'] = record['AMPNAMES']

490 inDict['inputRange'] = record['INPUT_RANGE']

491 inDict['fitRange'] = record['FIT_RANGE']

492 inDict['badAmps'] = record['BAD_AMPS']

493 inDict['fitParamNames'] = record['FIT_PARAM_NAMES']

494 inDict['freeFitParamNames'] = record['FREE_FIT_PARAM_NAMES']

495 inDict['gain'] = {str(n): v for n, v in zip(record['AMPNAMES'], record['GAIN'])}

496 inDict['aMatrix'] = record['A_MATRIX']

497 inDict['aMatrixSigma'] = record['A_MATRIX_SIGMA']

498 inDict['aMatrixModel'] = record['A_MATRIX_MODEL']

499 inDict['aMatrixSum'] = record['A_MATRIX_SUM']

500 inDict['aMatrixModelSum'] = record['A_MATRIX_MODEL_SUM']

501 inDict['modelNormalization'] = record['MODEL_NORMALIZATION']

502 inDict['usedPixels'] = record['USED_PIXELS']

503 inDict['fitParams'] = {str(n): v for n, v in zip(record['FIT_PARAM_NAMES'], record['FIT_PARAMS'])}

504 inDict['fitParamErrors'] = {str(n): v for n, v in zip(record['FIT_PARAM_NAMES'],

505 record['FIT_PARAM_ERRORS'])}

506 inDict['fitChi2'] = record['FIT_CHI2']

507 inDict['fitReducedChi2'] = record['FIT_REDUCED_CHI2']

508 inDict['fitParamCovMatrix'] = record['FIT_PARAM_COV_MATRIX']

509 inDict['ath'] = record['ATH']

510 inDict['athMinusBeta'] = record['ATH_MINUS_BETA']

511 inDict['aN'] = record['A_N']

512 inDict['aS'] = record['A_S']

513 inDict['aE'] = record['A_E']

514 inDict['aW'] = record['A_W']

515

516 if not calibVersion == 1.0:

517 cls().log.warning("Unkown version found for "

518 f"ElectrostaticBrighterFatterDistortionMatrix: {calibVersion}. ")

519

520 # Check for newer versions

521 if calibVersion < 1.1:

522 inDict['availableFilters'] = list()

523 inDict['aMatrixModelByFilter'] = dict()

524 inDict['athByFilter'] = dict()

525 inDict['athMinusBetaByFilter'] = dict()

526 inDict['aNByFilter'] = dict()

527 inDict['aSByFilter'] = dict()

528 inDict['aEByFilter'] = dict()

529 inDict['aWByFilter'] = dict()

530 else:

531 inDict['availableFilters'] = record['AVAILABLE_FILTERS']

532 inDict['aMatrixModelByFilter'] = record['A_MATRIX_MODEL_BY_FILTER']

533 inDict['athByFilter'] = record['ATH_BY_FILTER']

534 inDict['athMinusBetaByFilter'] = record['ATH_MINUS_BETA_BY_FILTER']

535 inDict['aNByFilter'] = record['A_N_BY_FILTER']

536 inDict['aSByFilter'] = record['A_S_BY_FILTER']

537 inDict['aEByFilter'] = record['A_E_BY_FILTER']

538 inDict['aWByFilter'] = record['A_W_BY_FILTER']

539

540 return cls.fromDict(inDict)

541

542 def toTable(self):

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

544 calibration.

545

546 The list of tables should create an identical calibration

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

548

549 Returns

550 -------

551 tableList : `list` [`lsst.afw.table.Table`]

552 List of tables containing the crosstalk calibration

553 information.

554

555 """

556 tableList = []

557 self.updateMetadata()

558

559 recordDict = {

560 'AMPNAMES': self.ampNames,

561 'INPUT_RANGE': self.inputRange,

562 'FIT_RANGE': self.fitRange,

563 'BAD_AMPS': self.badAmps,

564 'GAIN': list(self.gain.values()),

565 'A_MATRIX': self.aMatrix.ravel(),

566 'A_MATRIX_SIGMA': self.aMatrixSigma.ravel(),

567 'A_MATRIX_MODEL': self.aMatrixModel.ravel(),

568 'A_MATRIX_SUM': self.aMatrixSum,

569 'A_MATRIX_MODEL_SUM': self.aMatrixModelSum,

570 'MODEL_NORMALIZATION': self.modelNormalization,

571 'USED_PIXELS': self.usedPixels.ravel(),

572 'FIT_PARAM_NAMES': self.fitParamNames,

573 'FREE_FIT_PARAM_NAMES': self.freeFitParamNames,

574 'FIT_PARAMS': list(self.fitParams.values()),

575 'FIT_PARAM_ERRORS': list(self.fitParamErrors.values()),

576 'FIT_CHI2': self.fitChi2,

577 'FIT_REDUCED_CHI2': self.fitReducedChi2,

578 'FIT_PARAM_COV_MATRIX': self.fitParamCovMatrix.ravel(),

579 'ATH': self.ath.ravel(),

580 'ATH_MINUS_BETA': self.athMinusBeta.ravel(),

581 'A_N': self.aN.ravel(),

582 'A_S': self.aS.ravel(),

583 'A_E': self.aE.ravel(),

584 'A_W': self.aW.ravel(),

585 'AVAILABLE_FILTERS': self.availableFilters,

586 'A_MATRIX_MODEL_BY_FILTER': list(self.aMatrixModelByFilter.values()),

587 'ATH_BY_FILTER': list(self.athByFilter.values()),

588 'ATH_MINUS_BETA_BY_FILTER': list(self.athMinusBetaByFilter.values()),

589 'A_N_BY_FILTER': list(self.aNByFilter.values()),

590 'A_S_BY_FILTER': list(self.aSByFilter.values()),

591 'A_E_BY_FILTER': list(self.aEByFilter.values()),

592 'A_W_BY_FILTER': list(self.aWByFilter.values()),

593 }

594

595 catalogList = [recordDict]

596 catalog = Table(catalogList)

597

598 inMeta = self.getMetadata().toDict()

599 outMeta = {k: v for k, v in inMeta.items() if v is not None}

600 outMeta.update({k: "" for k, v in inMeta.items() if v is None})

601 catalog.meta = outMeta

602 tableList.append(catalog)

603

604 return tableList

605

606

607class CustomFFTConvolution(object):

608 """

609 A class that performs image convolutions in Fourier space, using pyfftw.

610 The constructor takes images as arguments and creates the FFTW plans.

611 The convolutions are performed by the __call__ routine.

612 This is faster than scipy.signal.fftconvolve, and it saves some transforms

613 by allowing the same image to be convolved with several kernels.

614 pyfftw does not accommodate float32 images, so everything

615 should be double precision.

616

617 Code adaped from :

618 https://stackoverflow.com/questions/14786920/convolution-of-two-three-dimensional-arrays-with-padding-on-one-side-too-slow

619 Code posted by Henry Gomersal

620 """

621 def __init__(self, im, kernel, threads=1):

622 # Compute the minimum size of the convolution result.

623 shape = (np.array(im.shape) + np.array(kernel.shape)) - 1

624

625 # Find the next larger "fast size" for FFT computation.

626 # This can provide a significant speedup.

627 shape = np.array([pyfftw.next_fast_len(s) for s in shape])

628

629 # Create FFTW building objects for the image and kernel.

630 self.fftImageObj = pyfftw.builders.rfftn(im, s=shape, threads=threads)

631 self.fftKernelObj = pyfftw.builders.rfftn(kernel, s=shape, threads=threads)

632 self.ifftObj = pyfftw.builders.irfftn(

633 self.fftImageObj.output_array,

634 s=shape,

635 threads=threads,

636 )

637

638 def __call__(self, im, kernels):

639 """

640 Perform the convolution and trim the result to the

641 size of the input image. If kernels is a list, then

642 the routine returns a list of corresponding

643 convolutions.

644 """

645 # Handle both a list of kernels and a single kernel.

646 ks = [kernels] if type(kernels) is not list else kernels

647 convolutions = []

648 for k in ks:

649 # Transform the image and the kernel using FFT.

650 tim = self.fftImageObj(im)

651 tk = self.fftKernelObj(k)

652

653 # Multiply in Fourier space and perform the inverse FFT.

654 convolution = self.ifftObj(tim * tk)

655 # Trim the result to match the input image size, following

656 # the 'same' policy of scipy.signal.fftconvolve.

657 oy = k.shape[0] // 2

658 ox = k.shape[1] // 2

659 convolutions.append(convolution[oy:oy + im.shape[0], ox:ox + im.shape[1]].copy())

660 # Return a single convolution if kernels was

661 # not a list, otherwise return the list.

662 return convolutions[0] if type(kernels) is not list else convolutions

663

664

665def electrostaticBrighterFatterCorrection(log, exposure, electroBfDistortionMatrix, applyGain,

666 gains=None, applyFilterCorrection=False):

667 """

668 Evaluates the correction of CCD images affected by the

669 brighter-fatter effect, as described in

670 https://arxiv.org/abs/2301.03274. Requires as input the result of

671 an electrostatic fit to flat covariance data (or any other

672 determination of pixel boundary shifts under the influence of a

673 single electron).

674

675 The filename refers to an input tuple that contains the

676 boundary shifts for one electron. This file is produced by an

677 electrostatic fit to data extracted from flat-field statistics,

678 implemented in https://gitlab.in2p3.fr/astier/bfptc/tools/fit_cov.py.

679 """

680

681 # Use the symmetrize function to fill the four quadrants for each kernel

682 r = electroBfDistortionMatrix.fitRange - 1

683 aN = electroBfDistortionMatrix.aN

684 aS = electroBfDistortionMatrix.aS

685 aE = electroBfDistortionMatrix.aE

686 aW = electroBfDistortionMatrix.aW

687

688 # Find the physicalFilter associated with the exposure

689 physicalFilter = exposure.metadata["FILTER"]

690

691 if applyFilterCorrection:

692 if (physicalFilter is not None) and \

693 (physicalFilter in electroBfDistortionMatrix.availableFilters):

694 # We have a pre-computed distortion matrix for

695 # the exposure's filter

696 aN = electroBfDistortionMatrix.aNByFilter[physicalFilter]

697 aS = electroBfDistortionMatrix.aSByFilter[physicalFilter]

698 aE = electroBfDistortionMatrix.aEByFilter[physicalFilter]

699 aW = electroBfDistortionMatrix.aWByFilter[physicalFilter]

700 else:

701 # There is no pre-computed distortion matrix

702 # for the exposure's filter

703 log.warning(f"Filter {exposure.metadata['FILTER']} not found in "

704 "electroBfDistortionMatrix.availableFilters. "

705 "Defaulting to the using the `generic` distortion "

706 "matrix, computed assuming all photons are converted "

707 "to electrons at the surface of the silicon detector.")

708

709 # Initialize kN and kE arrays

710 kN = np.zeros((2 * r + 1, 2 * r + 1))

711 kE = np.zeros_like(kN)

712

713 # Fill in the 4 quadrants for kN

714 kN[r:, r:] = aN # Quadrant 1 (bottom-right)

715 kN[:r+1, r:] = np.flipud(aN) # Quadrant 2 (top-right)

716 kN[r:, :r+1] = np.fliplr(aS) # Quadrant 3 (bottom-left)

717 kN[:r+1, :r+1] = np.flipud(np.fliplr(aS)) # Quadrant 4 (top-left)

718

719 # Fill in the 4 quadrants for kE

720 kE[r:, r:] = aE # Quadrant 1 (bottom-right)

721 kE[:r+1, r:] = np.flipud(aW) # Quadrant 2 (top-right)

722 kE[r:, :r+1] = np.fliplr(aE) # Quadrant 3 (bottom-left)

723 kE[:r+1, :r+1] = np.flipud(np.fliplr(aW)) # Quadrant 4 (top-left)

724

725 # Tweak the edges so that the sum rule applies.

726 kN[:, 0] = -kN[:, -1]

727 kE[0, :] = -kE[-1, :]

728

729 # We use the normalization of Guyonnet et al. (2015),

730 # which is compatible with the way the input file is produced.

731 # The factor 1/2 is due to the fact that the charge distribution at the end

732 # is twice the average, and the second 1/2 is due to

733 # charge interpolation.

734 kN *= 0.25

735 kE *= 0.25

736

737 # Indeed, i and j in the tuple refer to serial and parallel directions.

738 # In most Python code, the image reads im[j, i], so we transpose:

739 kN = kN.T

740 kE = kE.T

741

742 # Get the image to perform the correction

743 image = exposure.getMaskedImage().getImage()

744

745 # The image needs to be units of electrons/holes

746 with gainContext(exposure, image, applyGain, gains):

747 # Computes the correction and returns the "delta_image",

748 # which should be subtracted from "im" in order to undo the BF effect.

749 # The input image should be expressed in electrons

750 im = image.getArray().copy()

751 convolver = CustomFFTConvolution(im, kN)

752 convolutions = convolver(im, [kN, kE])

753

754 # The convolutions contain the boundary shifts (in pixel size units)

755 # for [horizontal, vertical] boundaries.

756 # We now compute the charge to move around.

757 delta = np.zeros_like(im)

758 boundaryCharge = np.zeros_like(im)

759

760 # Horizontal boundaries (parallel direction).

761 # We could use a more elaborate interpolator for estimating the

762 # charge on the boundary.

763 boundaryCharge[:-1, :] = im[1:, :] + im[:-1, :]

764 # boundaryCharge[1:-2,:] = (9./8.)*(I[2:-1,:]+I[1:-2,:] -

765 # (1./8.)*(I[0:-3,:]+I[3,:])

766

767 # The charge to move around is the

768 # product of the boundary shift (in pixel size units) times the

769 # charge on the boundary (in charge per pixel unit).

770 dq = boundaryCharge * convolutions[0]

771 delta += dq

772

773 # What is gained by a pixel is lost by its neighbor (the right one).

774 delta[1:, :] -= dq[:-1, :]

775

776 # Vertical boundaries.

777 boundaryCharge = np.zeros_like(im) # Reset to zero.

778 boundaryCharge[:, :-1] = im[:, 1:] + im[:, :-1]

779 dq = boundaryCharge * convolutions[1]

780 delta += dq

781

782 # What is gained by a pixel is lost by its neighbor.

783 delta[:, 1:] -= dq[:, :-1]

784

785 # TODO: One might check that delta.sum() ~ 0 (charge conservation).

786

787 # Apply the correction to the original image

788 exposure.image.array -= delta

789

790 return exposure

Coverage for python / lsst / ip / isr / electrostaticBrighterFatter.py: 7%

286 statements