Coverage for python / lsst / ip / isr / electrostaticBrighterFatter.py: 9%
231 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-14 23:58 +0000
1# This file is part of ip_isr.
2#
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.
8#
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 """
135 _OBSTYPE = 'BF_DISTORTION_MATRIX'
136 _SCHEMA = 'ELectrostatic brighter-fatter pixel distortion model'
137 _VERSION = 1.0
139 def __init__(self, camera=None, inputRange=1, fitRange=None, **kwargs):
140 """
141 Filename refers to an input tuple that contains the
142 boundary shifts for one electron. This file is produced by an
143 electrostatic fit to data extracted from flat-field statistics,
144 implemented in https://gitlab.in2p3.fr/astier/bfptc/tools/fit_cov.py
145 """
146 self.ampNames = []
147 self.inputRange = inputRange
148 if fitRange is None:
149 fitRange = inputRange
150 self.fitRange = fitRange
151 self.badAmps = list()
152 self.gain = dict()
153 self.aMatrix = np.full((inputRange, inputRange), np.nan)
154 self.aMatrixSigma = np.full((inputRange, inputRange), np.nan)
155 self.aMatrixModel = np.full((fitRange, fitRange), np.nan)
156 self.aMatrixSum = np.nan
157 self.aMatrixModelSum = np.nan
158 self.modelNormalization = [np.nan, np.nan]
159 self.usedPixels = np.full((inputRange, inputRange), np.nan)
160 self.fitParamNames = list()
161 self.freeFitParamNames = list()
162 self.fitParams = dict()
163 self.fitParamErrors = dict()
164 self.fitChi2 = np.nan
165 self.fitReducedChi2 = np.nan
166 self.fitParamCovMatrix = np.array([])
167 self.ath = np.full((fitRange, fitRange), np.nan)
168 self.athMinusBeta = np.full((fitRange, fitRange), np.nan)
169 self.aN = np.full((fitRange, fitRange), np.nan)
170 self.aS = np.full((fitRange, fitRange), np.nan)
171 self.aE = np.full((fitRange, fitRange), np.nan)
172 self.aW = np.full((fitRange, fitRange), np.nan)
174 super().__init__(**kwargs)
176 if camera:
177 self.initFromCamera(camera, detectorId=kwargs.get('detectorId', None))
179 self.requiredAttributes.update([
180 'inputRange', 'fitRange', 'badAmps', 'gain', 'aMatrix',
181 'aMatrixSigma', 'aMatrixModel', 'aMatrixSum', 'aMatrixModelSum',
182 'modelNormalization', 'usedPixels', 'fitParamNames',
183 'freeFitParamNames', 'fitParams', 'fitParamErrors', 'fitChi2',
184 'fitReducedChi2', 'fitParamCovMatrix', 'ath', 'athMinusBeta',
185 'aN', 'aS', 'aE', 'aW', 'ampNames',
186 ])
188 def updateMetadata(self, setDate=False, **kwargs):
189 """Update calibration metadata.
191 This calls the base class's method after ensuring the required
192 calibration keywords will be saved.
194 Parameters
195 ----------
196 setDate : `bool`, optional
197 Update the CALIBDATE fields in the metadata to the current
198 time. Defaults to False.
199 kwargs :
200 Other keyword parameters to set in the metadata.
201 """
202 kwargs['INPUT_RANGE'] = self.inputRange
203 kwargs['FIT_RANGE'] = self.fitRange
205 super().updateMetadata(setDate=setDate, **kwargs)
207 def initFromCamera(self, camera, detectorId=None):
208 """Initialize kernel structure from camera.
210 Parameters
211 ----------
212 camera : `lsst.afw.cameraGeom.Camera`
213 Camera to use to define geometry.
214 detectorId : `int`, optional
215 Index of the detector to generate.
217 Returns
218 -------
219 calib : `lsst.ip.isr.BrighterFatterKernel`
220 The initialized calibration.
222 Raises
223 ------
224 RuntimeError
225 Raised if no detectorId is supplied for a calibration with
226 ``level='AMP'``.
227 """
228 self._instrument = camera.getName()
230 if detectorId is not None:
231 detector = camera[detectorId]
232 self._detectorId = detectorId
233 self._detectorName = detector.getName()
234 self._detectorSerial = detector.getSerial()
236 for amp in detector:
237 ampName = amp.getName()
238 self.ampNames.append(ampName)
239 else:
240 raise RuntimeError("A detectorId must be supplied if "
241 "camera is supplied.")
243 return self
245 @classmethod
246 def fromDict(cls, dictionary):
247 """Construct a calibration from a dictionary of properties.
249 Parameters
250 ----------
251 dictionary : `dict`
252 Dictionary of properties.
254 Returns
255 -------
256 calib : `lsst.ip.isr.BrighterFatterKernel`
257 Constructed calibration.
259 Raises
260 ------
261 RuntimeError
262 Raised if the supplied dictionary is for a different
263 calibration.
264 Raised if the version of the supplied dictionary is 1.0.
265 """
266 calib = cls()
268 if calib._OBSTYPE != (found := dictionary['metadata']['OBSTYPE']):
269 raise RuntimeError(f"Incorrect brighter-fatter calibration supplied. Expected {calib._OBSTYPE}, "
270 f"found {found}")
271 calib.setMetadata(dictionary['metadata'])
272 calib.calibInfoFromDict(dictionary)
274 calib.ampNames = dictionary['ampNames']
275 inputRange = dictionary['inputRange']
276 fitRange = dictionary['fitRange']
277 calib.inputRange = inputRange
278 calib.fitRange = fitRange
279 calib.gain = dictionary['gain']
280 calib.badAmps = dictionary['badAmps']
281 calib.fitParamNames = dictionary['fitParamNames']
282 calib.freeFitParamNames = dictionary['freeFitParamNames']
283 calib.aMatrix = np.array(
284 dictionary['aMatrix'],
285 dtype=np.float64,
286 ).reshape(inputRange, inputRange)
287 calib.aMatrixSigma = np.array(
288 dictionary['aMatrixSigma'],
289 dtype=np.float64,
290 ).reshape(inputRange, inputRange)
291 calib.aMatrixModel = np.array(
292 dictionary['aMatrixSigma'],
293 dtype=np.float64,
294 ).reshape(fitRange, fitRange)
295 calib.aMatrixSum = float(dictionary['aMatrixSum'])
296 calib.aMatrixModelSum = float(dictionary['aMatrixModelSum'])
297 calib.modelNormalization = dictionary['modelNormalization']
298 calib.usedPixels = np.array(dictionary['usedPixels'])
299 calib.fitParams = dictionary['fitParams']
300 calib.fitParamErrors = dictionary['fitParamErrors']
301 calib.fitChi2 = float(dictionary['fitChi2'])
302 calib.fitReducedChi2 = float(dictionary['fitReducedChi2'])
303 calib.fitParamCovMatrix = np.array(
304 dictionary['fitParamCovMatrix'],
305 dtype=np.float64,
306 ).reshape(len(calib.freeFitParamNames), len(calib.freeFitParamNames))
307 calib.ath = np.array(
308 dictionary['ath'],
309 dtype=np.float64,
310 ).reshape(fitRange, fitRange)
311 calib.athMinusBeta = np.array(
312 dictionary['athMinusBeta'],
313 dtype=np.float64,
314 ).reshape(fitRange, fitRange)
315 calib.aN = np.array(
316 dictionary['aN'],
317 dtype=np.float64,
318 ).reshape(fitRange, fitRange)
319 calib.aS = np.array(
320 dictionary['aS'],
321 dtype=np.float64,
322 ).reshape(fitRange, fitRange)
323 calib.aE = np.array(
324 dictionary['aE'],
325 dtype=np.float64,
326 ).reshape(fitRange, fitRange)
327 calib.aW = np.array(
328 dictionary['aW'],
329 dtype=np.float64,
330 ).reshape(fitRange, fitRange)
332 calib.updateMetadata()
333 return calib
335 def toDict(self):
336 """Return a dictionary containing the calibration properties.
338 The dictionary should be able to be round-tripped through
339 `fromDict`.
341 Returns
342 -------
343 dictionary : `dict`
344 Dictionary of properties.
345 """
346 self.updateMetadata()
348 outDict = dict()
349 metadata = self.getMetadata()
350 outDict['metadata'] = metadata
352 outDict['ampNames'] = self.ampNames
353 outDict['inputRange'] = self.inputRange
354 outDict['fitRange'] = self.fitRange
355 outDict['badAmps'] = self.badAmps
356 outDict['fitParamNames'] = self.fitParamNames
357 outDict['freeFitParamNames'] = self.freeFitParamNames
358 outDict['gain'] = self.gain
359 outDict['aMatrix'] = self.aMatrix.ravel().tolist()
360 outDict['aMatrixSigma'] = self.aMatrixSigma.ravel().tolist()
361 outDict['aMatrixModel'] = self.aMatrixModel.ravel().tolist()
362 outDict['aMatrixSum'] = self.aMatrixSum
363 outDict['aMatrixModelSum'] = self.aMatrixModelSum
364 outDict['aMatrixModel'] = self.aMatrixModel.ravel().tolist()
365 outDict['modelNormalization'] = self.modelNormalization
366 outDict['fitParamCovMatrix'] = self.fitParamCovMatrix.ravel().tolist()
367 outDict['usedPixels'] = self.usedPixels.ravel().tolist()
368 outDict['fitParams'] = self.fitParams
369 outDict['fitParamErrors'] = self.fitParamErrors
370 outDict['fitChi2'] = self.fitChi2
371 outDict['fitReducedChi2'] = self.fitReducedChi2
372 outDict['ath'] = self.ath.ravel().tolist()
373 outDict['athMinusBeta'] = self.athMinusBeta.ravel().tolist()
374 outDict['aN'] = self.aN.ravel().tolist()
375 outDict['aS'] = self.aS.ravel().tolist()
376 outDict['aE'] = self.aE.ravel().tolist()
377 outDict['aW'] = self.aW.ravel().tolist()
379 return outDict
381 @classmethod
382 def fromTable(cls, tableList):
383 """Construct calibration from a list of tables.
385 This method uses the `fromDict` method to create the
386 calibration, after constructing an appropriate dictionary from
387 the input tables.
389 Parameters
390 ----------
391 tableList : `list` [`astropy.table.Table`]
392 List of tables to use to construct the brighter-fatter
393 calibration.
395 Returns
396 -------
397 calib : `lsst.ip.isr.BrighterFatterKernel`
398 The calibration defined in the tables.
399 """
400 record = tableList[0][0]
402 metadata = record.meta
403 calibVersion = metadata['BF_DISTORTION_MATRIX_VERSION']
405 # Initialize inDict with all expected keys
406 # and empty values of the corresponding type
407 inDict = dict()
408 inDict['metadata'] = metadata
409 inDict['ampNames'] = record['AMPNAMES']
410 inDict['inputRange'] = record['INPUT_RANGE']
411 inDict['fitRange'] = record['FIT_RANGE']
412 inDict['badAmps'] = record['BAD_AMPS']
413 inDict['fitParamNames'] = record['FIT_PARAM_NAMES']
414 inDict['freeFitParamNames'] = record['FREE_FIT_PARAM_NAMES']
415 inDict['gain'] = {str(n): v for n, v in zip(record['AMPNAMES'], record['GAIN'])}
416 inDict['aMatrix'] = record['A_MATRIX']
417 inDict['aMatrixSigma'] = record['A_MATRIX_SIGMA']
418 inDict['aMatrixModel'] = record['A_MATRIX_MODEL']
419 inDict['aMatrixSum'] = record['A_MATRIX_SUM']
420 inDict['aMatrixModelSum'] = record['A_MATRIX_MODEL_SUM']
421 inDict['modelNormalization'] = record['MODEL_NORMALIZATION']
422 inDict['usedPixels'] = record['USED_PIXELS']
423 inDict['fitParams'] = {str(n): v for n, v in zip(record['FIT_PARAM_NAMES'], record['FIT_PARAMS'])}
424 inDict['fitParamErrors'] = {str(n): v for n, v in zip(record['FIT_PARAM_NAMES'],
425 record['FIT_PARAM_ERRORS'])}
426 inDict['fitChi2'] = record['FIT_CHI2']
427 inDict['fitReducedChi2'] = record['FIT_REDUCED_CHI2']
428 inDict['fitParamCovMatrix'] = record['FIT_PARAM_COV_MATRIX']
429 inDict['ath'] = record['ATH']
430 inDict['athMinusBeta'] = record['ATH_MINUS_BETA']
431 inDict['aN'] = record['A_N']
432 inDict['aS'] = record['A_S']
433 inDict['aE'] = record['A_E']
434 inDict['aW'] = record['A_W']
436 if not calibVersion == 1.0:
437 cls().log.warning("Unkown version found for "
438 f"ElectrostaticBrighterFatterDistortionMatrix: {calibVersion}. ")
440 # Check for newer versions, but there are none...
442 return cls.fromDict(inDict)
444 def toTable(self):
445 """Construct a list of tables containing the information in this
446 calibration.
448 The list of tables should create an identical calibration
449 after being passed to this class's fromTable method.
451 Returns
452 -------
453 tableList : `list` [`lsst.afw.table.Table`]
454 List of tables containing the crosstalk calibration
455 information.
457 """
458 tableList = []
459 self.updateMetadata()
461 recordDict = {
462 'AMPNAMES': self.ampNames,
463 'INPUT_RANGE': self.inputRange,
464 'FIT_RANGE': self.fitRange,
465 'BAD_AMPS': self.badAmps,
466 'GAIN': list(self.gain.values()),
467 'A_MATRIX': self.aMatrix.ravel(),
468 'A_MATRIX_SIGMA': self.aMatrixSigma.ravel(),
469 'A_MATRIX_MODEL': self.aMatrixModel.ravel(),
470 'A_MATRIX_SUM': self.aMatrixSum,
471 'A_MATRIX_MODEL_SUM': self.aMatrixModelSum,
472 'MODEL_NORMALIZATION': self.modelNormalization,
473 'USED_PIXELS': self.usedPixels.ravel(),
474 'FIT_PARAM_NAMES': self.fitParamNames,
475 'FREE_FIT_PARAM_NAMES': self.freeFitParamNames,
476 'FIT_PARAMS': list(self.fitParams.values()),
477 'FIT_PARAM_ERRORS': list(self.fitParamErrors.values()),
478 'FIT_CHI2': self.fitChi2,
479 'FIT_REDUCED_CHI2': self.fitReducedChi2,
480 'FIT_PARAM_COV_MATRIX': self.fitParamCovMatrix.ravel(),
481 'ATH': self.ath.ravel(),
482 'ATH_MINUS_BETA': self.athMinusBeta.ravel(),
483 'A_N': self.aN.ravel(),
484 'A_S': self.aS.ravel(),
485 'A_E': self.aE.ravel(),
486 'A_W': self.aW.ravel(),
487 }
489 catalogList = [recordDict]
490 catalog = Table(catalogList)
492 inMeta = self.getMetadata().toDict()
493 outMeta = {k: v for k, v in inMeta.items() if v is not None}
494 outMeta.update({k: "" for k, v in inMeta.items() if v is None})
495 catalog.meta = outMeta
496 tableList.append(catalog)
498 return tableList
501class CustomFFTConvolution(object):
502 """
503 A class that performs image convolutions in Fourier space, using pyfftw.
504 The constructor takes images as arguments and creates the FFTW plans.
505 The convolutions are performed by the __call__ routine.
506 This is faster than scipy.signal.fftconvolve, and it saves some transforms
507 by allowing the same image to be convolved with several kernels.
508 pyfftw does not accommodate float32 images, so everything
509 should be double precision.
511 Code adaped from :
512 https://stackoverflow.com/questions/14786920/convolution-of-two-three-dimensional-arrays-with-padding-on-one-side-too-slow
513 Code posted by Henry Gomersal
514 """
515 def __init__(self, im, kernel, threads=1):
516 # Compute the minimum size of the convolution result.
517 shape = (np.array(im.shape) + np.array(kernel.shape)) - 1
519 # Find the next larger "fast size" for FFT computation.
520 # This can provide a significant speedup.
521 shape = np.array([pyfftw.next_fast_len(s) for s in shape])
523 # Create FFTW building objects for the image and kernel.
524 self.fftImageObj = pyfftw.builders.rfftn(im, s=shape, threads=threads)
525 self.fftKernelObj = pyfftw.builders.rfftn(kernel, s=shape, threads=threads)
526 self.ifftObj = pyfftw.builders.irfftn(
527 self.fftImageObj.get_output_array(),
528 s=shape,
529 threads=threads,
530 )
532 def __call__(self, im, kernels):
533 """
534 Perform the convolution and trim the result to the
535 size of the input image. If kernels is a list, then
536 the routine returns a list of corresponding
537 convolutions.
538 """
539 # Handle both a list of kernels and a single kernel.
540 ks = [kernels] if type(kernels) is not list else kernels
541 convolutions = []
542 for k in ks:
543 # Transform the image and the kernel using FFT.
544 tim = self.fftImageObj(im)
545 tk = self.fftKernelObj(k)
547 # Multiply in Fourier space and perform the inverse FFT.
548 convolution = self.ifftObj(tim * tk)
549 # Trim the result to match the input image size, following
550 # the 'same' policy of scipy.signal.fftconvolve.
551 oy = k.shape[0] // 2
552 ox = k.shape[1] // 2
553 convolutions.append(convolution[oy:oy + im.shape[0], ox:ox + im.shape[1]].copy())
554 # Return a single convolution if kernels was
555 # not a list, otherwise return the list.
556 return convolutions[0] if type(kernels) is not list else convolutions
559def electrostaticBrighterFatterCorrection(exposure, electroBfDistortionMatrix, applyGain, gains=None):
560 """
561 Evaluates the correction of CCD images affected by the
562 brighter-fatter effect, as described in
563 https://arxiv.org/abs/2301.03274. Requires as input the result of
564 an electrostatic fit to flat covariance data (or any other
565 determination of pixel boundary shifts under the influence of a
566 single electron).
568 The filename refers to an input tuple that contains the
569 boundary shifts for one electron. This file is produced by an
570 electrostatic fit to data extracted from flat-field statistics,
571 implemented in https://gitlab.in2p3.fr/astier/bfptc/tools/fit_cov.py.
572 """
574 # Use the symmetrize function to fill the four quadrants for each kernel
575 r = electroBfDistortionMatrix.fitRange - 1
576 aN = electroBfDistortionMatrix.aN
577 aS = electroBfDistortionMatrix.aS
578 aE = electroBfDistortionMatrix.aE
579 aW = electroBfDistortionMatrix.aW
581 # Initialize kN and kE arrays
582 kN = np.zeros((2 * r + 1, 2 * r + 1))
583 kE = np.zeros_like(kN)
585 # Fill in the 4 quadrants for kN
586 kN[r:, r:] = aN # Quadrant 1 (bottom-right)
587 kN[:r+1, r:] = np.flipud(aN) # Quadrant 2 (top-right)
588 kN[r:, :r+1] = np.fliplr(aS) # Quadrant 3 (bottom-left)
589 kN[:r+1, :r+1] = np.flipud(np.fliplr(aS)) # Quadrant 4 (top-left)
591 # Fill in the 4 quadrants for kE
592 kE[r:, r:] = aE # Quadrant 1 (bottom-right)
593 kE[:r+1, r:] = np.flipud(aW) # Quadrant 2 (top-right)
594 kE[r:, :r+1] = np.fliplr(aE) # Quadrant 3 (bottom-left)
595 kE[:r+1, :r+1] = np.flipud(np.fliplr(aW)) # Quadrant 4 (top-left)
597 # Tweak the edges so that the sum rule applies.
598 kN[:, 0] = -kN[:, -1]
599 kE[0, :] = -kE[-1, :]
601 # We use the normalization of Guyonnet et al. (2015),
602 # which is compatible with the way the input file is produced.
603 # The factor 1/2 is due to the fact that the charge distribution at the end
604 # is twice the average, and the second 1/2 is due to
605 # charge interpolation.
606 kN *= 0.25
607 kE *= 0.25
609 # Indeed, i and j in the tuple refer to serial and parallel directions.
610 # In most Python code, the image reads im[j, i], so we transpose:
611 kN = kN.T
612 kE = kE.T
614 # Get the image to perform the correction
615 image = exposure.getMaskedImage().getImage()
617 # The image needs to be units of electrons/holes
618 with gainContext(exposure, image, applyGain, gains):
619 # Computes the correction and returns the "delta_image",
620 # which should be subtracted from "im" in order to undo the BF effect.
621 # The input image should be expressed in electrons
622 im = image.getArray().copy()
623 convolver = CustomFFTConvolution(im, kN)
624 convolutions = convolver(im, [kN, kE])
626 # The convolutions contain the boundary shifts (in pixel size units)
627 # for [horizontal, vertical] boundaries.
628 # We now compute the charge to move around.
629 delta = np.zeros_like(im)
630 boundaryCharge = np.zeros_like(im)
632 # Horizontal boundaries (parallel direction).
633 # We could use a more elaborate interpolator for estimating the
634 # charge on the boundary.
635 boundaryCharge[:-1, :] = im[1:, :] + im[:-1, :]
636 # boundaryCharge[1:-2,:] = (9./8.)*(I[2:-1,:]+I[1:-2,:] -
637 # (1./8.)*(I[0:-3,:]+I[3,:])
639 # The charge to move around is the
640 # product of the boundary shift (in pixel size units) times the
641 # charge on the boundary (in charge per pixel unit).
642 dq = boundaryCharge * convolutions[0]
643 delta += dq
645 # What is gained by a pixel is lost by its neighbor (the right one).
646 delta[1:, :] -= dq[:-1, :]
648 # Vertical boundaries.
649 boundaryCharge = np.zeros_like(im) # Reset to zero.
650 boundaryCharge[:, :-1] = im[:, 1:] + im[:, :-1]
651 dq = boundaryCharge * convolutions[1]
652 delta += dq
654 # What is gained by a pixel is lost by its neighbor.
655 delta[:, 1:] -= dq[:, :-1]
657 # TODO: One might check that delta.sum() ~ 0 (charge conservation).
659 # Apply the correction to the original image
660 exposure.image.array -= delta
662 return exposure