electrostaticBrighterFatter.py

Go to the documentation of this file.

# This file is part of ip_isr.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
#
"""Brighter Fatter Kernel calibration definition."""
 
 
__all__ = [
    'ElectrostaticBrighterFatterDistortionMatrix',
    'electrostaticBrighterFatterCorrection',
]
 
 
import pyfftw
import numpy as np
from astropy.table import Table
 
from . import IsrCalib
from .isrFunctions import gainContext
 
 
class ElectrostaticBrighterFatterDistortionMatrix(IsrCalib):
    """Calibration of brighter-fatter kernels for an instrument.
 
    ampKernels are the kernels for each amplifier in a detector, as
    generated by having ``level == 'AMP'``.
 
    detectorKernel is the kernel generated for a detector as a
    whole, as generated by having ``level == 'DETECTOR'``.
 
    makeDetectorKernelFromAmpwiseKernels is a method to generate the
    kernel for a detector, constructed by averaging together the
    ampwise kernels in the detector.  The existing application code is
    only defined for kernels with ``level == 'DETECTOR'``, so this method
    is used if the supplied kernel was built with ``level == 'AMP'``.
 
    Parameters
    ----------
    camera : `lsst.afw.cameraGeom.Camera`
        Camera describing detector geometry.
    level : `str`
        Level the kernels will be generated for.
    log : `logging.Logger`, optional
        Log to write messages to.
    **kwargs :
        Parameters to pass to parent constructor.
 
    Notes
    -----
    Version 1.1 adds the `expIdMask` property, and substitutes
    `means` and `variances` for `rawMeans` and `rawVariances`
    from the PTC dataset.
 
    inputRange : `int`
        The size of the input aMatrix shape in each dimension.
    fitRange : `int`
        The size of the input aMatrix shape in each dimension that is
        used for fitting the electrostatic model. Must be less than or
        equal to inputRange.
    badAmps : `list`, [`str`]
        List of bad amplifiers names.
    gain : `dict`, [`str`,`float`]
        Dictionary keyed by amp names containing the gains inherited
        from the inputPTC.
    aMatrix : `numpy.ndarray`
        The average aMatrix inherited from the inputPTC's good amplifiers.
    aMatrixSigma : `numpy.ndarray`
        The uncertainty matrix used to weight the fit.
    aMatrixModel : `numpy.ndarray`
        The modeled aMatrix based on the electrostatic fit parameters.
    aMatrixSum : `float`
        The sum of the symmetrized aMatrix.
    aMatrixModelSum : `float`
        The sum of the symmetrized aMatrixModel.
    modelNormalization : `list`
        A two element array of the multiplicative and additive
        normalization to the aMatrixModel.
    usedPixels : `numpy.ndarray`, [`bool`]
        Array of shape like aMatrix containing the mask indicating which
        elements of the input aMatrix were used to fit the electrostatic
        model.
    fitParamNames : `list`, [`str`]
        List of all the parameter names in the electrostatic fit.
    freeFitParamNames : `list`, [`str`]
        List of the parameter names that were allowed to vary during
        the electrostatic fit.
    fitParams : `dict`, [`str`, `float`]
        Dictionary containing each named parameter and its final
        fitted value.
    fitParamErrors : `dict`, [`str`, `float`]
        Dictionary containing each named parameter and its
        estimated fitting error.
    fitChi2 : `float`
        The computed chi squared between the data and the
        final model.
    fitReducedChi2 : `float`
        The computed reduced chi squared between the data
        and the final model.
    fitParamCovMatrix : `numpy.ndarray`
        The estimated covariance matrix between all fit parameters.
    ath : `numpy.ndarray`
        something...
    athMinusBeta : `numpy.ndarray`
        something...
    aN : `numpy.ndarray`
        Array of shape (fitRange, fitRange) containing the
        computed `North` component of the pixel boundary shift.
    aS : `numpy.ndarray`
        Array of shape (fitRange, fitRange) containing the
        computed `South` component of the pixel boundary shift.
    aE : `numpy.ndarray`
        Array of shape (fitRange, fitRange) containing the
        computed `East` component of the pixel boundary shift.
    aW : `numpy.ndarray`
        Array of shape (fitRange, fitRange) containing the
        computed `West` component of the pixel boundary shift.
    """
    _OBSTYPE = 'BF_DISTORTION_MATRIX'
    _SCHEMA = 'ELectrostatic brighter-fatter pixel distortion model'
    _VERSION = 1.0
 
    def __init__(self, camera=None, inputRange=1, fitRange=None, **kwargs):
        """
        Filename  refers to an input tuple that contains the
        boundary shifts for one electron. This file is produced by an
        electrostatic fit to data extracted from flat-field statistics,
        implemented in https://gitlab.in2p3.fr/astier/bfptc/tools/fit_cov.py
        """
        self.ampNames = []
        self.inputRange = inputRange
        if fitRange is None:
            fitRange = inputRange
        self.fitRange = fitRange
        self.badAmps = list()
        self.gain = dict()
        self.aMatrix = np.full((inputRange, inputRange), np.nan)
        self.aMatrixSigma = np.full((inputRange, inputRange), np.nan)
        self.aMatrixModel = np.full((fitRange, fitRange), np.nan)
        self.aMatrixSum = np.nan
        self.aMatrixModelSum = np.nan
        self.modelNormalization = [np.nan, np.nan]
        self.usedPixels = np.full((inputRange, inputRange), np.nan)
        self.fitParamNames = list()
        self.freeFitParamNames = list()
        self.fitParams = dict()
        self.fitParamErrors = dict()
        self.fitChi2 = np.nan
        self.fitReducedChi2 = np.nan
        self.fitParamCovMatrix = np.array([])
        self.ath = np.full((fitRange, fitRange), np.nan)
        self.athMinusBeta = np.full((fitRange, fitRange), np.nan)
        self.aN = np.full((fitRange, fitRange), np.nan)
        self.aS = np.full((fitRange, fitRange), np.nan)
        self.aE = np.full((fitRange, fitRange), np.nan)
        self.aW = np.full((fitRange, fitRange), np.nan)
 
        super().__init__(**kwargs)
 
        if camera:
            self.initFromCamera(camera, detectorId=kwargs.get('detectorId', None))
 
        self.requiredAttributes.update([
            'inputRange', 'fitRange', 'badAmps', 'gain', 'aMatrix',
            'aMatrixSigma', 'aMatrixModel', 'aMatrixSum', 'aMatrixModelSum',
            'modelNormalization', 'usedPixels', 'fitParamNames',
            'freeFitParamNames', 'fitParams', 'fitParamErrors', 'fitChi2',
            'fitReducedChi2', 'fitParamCovMatrix', 'ath', 'athMinusBeta',
            'aN', 'aS', 'aE', 'aW', 'ampNames',
        ])
 
    def updateMetadata(self, setDate=False, **kwargs):
        """Update calibration metadata.
 
        This calls the base class's method after ensuring the required
        calibration keywords will be saved.
 
        Parameters
        ----------
        setDate : `bool`, optional
            Update the CALIBDATE fields in the metadata to the current
            time. Defaults to False.
        kwargs :
            Other keyword parameters to set in the metadata.
        """
        kwargs['INPUT_RANGE'] = self.inputRange
        kwargs['FIT_RANGE'] = self.fitRange
 
        super().updateMetadata(setDate=setDate, **kwargs)
 
    def initFromCamera(self, camera, detectorId=None):
        """Initialize kernel structure from camera.
 
        Parameters
        ----------
        camera : `lsst.afw.cameraGeom.Camera`
            Camera to use to define geometry.
        detectorId : `int`, optional
            Index of the detector to generate.
 
        Returns
        -------
        calib : `lsst.ip.isr.BrighterFatterKernel`
            The initialized calibration.
 
        Raises
        ------
        RuntimeError
            Raised if no detectorId is supplied for a calibration with
            ``level='AMP'``.
        """
        self._instrument = camera.getName()
 
        if detectorId is not None:
            detector = camera[detectorId]
            self._detectorId = detectorId
            self._detectorName = detector.getName()
            self._detectorSerial = detector.getSerial()
 
            for amp in detector:
                ampName = amp.getName()
                self.ampNames.append(ampName)
        else:
            raise RuntimeError("A detectorId must be supplied if "
                               "camera is supplied.")
 
        return self
 
    @classmethod
    def fromDict(cls, dictionary):
        """Construct a calibration from a dictionary of properties.
 
        Parameters
        ----------
        dictionary : `dict`
            Dictionary of properties.
 
        Returns
        -------
        calib : `lsst.ip.isr.BrighterFatterKernel`
            Constructed calibration.
 
        Raises
        ------
        RuntimeError
            Raised if the supplied dictionary is for a different
            calibration.
            Raised if the version of the supplied dictionary is 1.0.
        """
        calib = cls()
 
        if calib._OBSTYPE != (found := dictionary['metadata']['OBSTYPE']):
            raise RuntimeError(f"Incorrect brighter-fatter calibration supplied.  Expected {calib._OBSTYPE}, "
                               f"found {found}")
        calib.setMetadata(dictionary['metadata'])
        calib.calibInfoFromDict(dictionary)
 
        calib.ampNames = dictionary['ampNames']
        inputRange = dictionary['inputRange']
        fitRange = dictionary['fitRange']
        calib.inputRange = inputRange
        calib.fitRange = fitRange
        calib.gain = dictionary['gain']
        calib.badAmps = dictionary['badAmps']
        calib.fitParamNames = dictionary['fitParamNames']
        calib.freeFitParamNames = dictionary['freeFitParamNames']
        calib.aMatrix = np.array(
            dictionary['aMatrix'],
            dtype=np.float64,
        ).reshape(inputRange, inputRange)
        calib.aMatrixSigma = np.array(
            dictionary['aMatrixSigma'],
            dtype=np.float64,
        ).reshape(inputRange, inputRange)
        calib.aMatrixModel = np.array(
            dictionary['aMatrixSigma'],
            dtype=np.float64,
        ).reshape(fitRange, fitRange)
        calib.aMatrixSum = float(dictionary['aMatrixSum'])
        calib.aMatrixModelSum = float(dictionary['aMatrixModelSum'])
        calib.modelNormalization = dictionary['modelNormalization']
        calib.usedPixels = np.array(dictionary['usedPixels'])
        calib.fitParams = dictionary['fitParams']
        calib.fitParamErrors = dictionary['fitParamErrors']
        calib.fitChi2 = float(dictionary['fitChi2'])
        calib.fitReducedChi2 = float(dictionary['fitReducedChi2'])
        calib.fitParamCovMatrix = np.array(
            dictionary['fitParamCovMatrix'],
            dtype=np.float64,
        ).reshape(len(calib.freeFitParamNames), len(calib.freeFitParamNames))
        calib.ath = np.array(
            dictionary['ath'],
            dtype=np.float64,
        ).reshape(fitRange, fitRange)
        calib.athMinusBeta = np.array(
            dictionary['athMinusBeta'],
            dtype=np.float64,
        ).reshape(fitRange, fitRange)
        calib.aN = np.array(
            dictionary['aN'],
            dtype=np.float64,
        ).reshape(fitRange, fitRange)
        calib.aS = np.array(
            dictionary['aS'],
            dtype=np.float64,
        ).reshape(fitRange, fitRange)
        calib.aE = np.array(
            dictionary['aE'],
            dtype=np.float64,
        ).reshape(fitRange, fitRange)
        calib.aW = np.array(
            dictionary['aW'],
            dtype=np.float64,
        ).reshape(fitRange, fitRange)
 
        calib.updateMetadata()
        return calib
 
    def toDict(self):
        """Return a dictionary containing the calibration properties.
 
        The dictionary should be able to be round-tripped through
        `fromDict`.
 
        Returns
        -------
        dictionary : `dict`
            Dictionary of properties.
        """
        self.updateMetadata()
 
        outDict = dict()
        metadata = self.getMetadata()
        outDict['metadata'] = metadata
 
        outDict['ampNames'] = self.ampNames
        outDict['inputRange'] = self.inputRange
        outDict['fitRange'] = self.fitRange
        outDict['badAmps'] = self.badAmps
        outDict['fitParamNames'] = self.fitParamNames
        outDict['freeFitParamNames'] = self.freeFitParamNames
        outDict['gain'] = self.gain
        outDict['aMatrix'] = self.aMatrix.ravel().tolist()
        outDict['aMatrixSigma'] = self.aMatrixSigma.ravel().tolist()
        outDict['aMatrixModel'] = self.aMatrixModel.ravel().tolist()
        outDict['aMatrixSum'] = self.aMatrixSum
        outDict['aMatrixModelSum'] = self.aMatrixModelSum
        outDict['aMatrixModel'] = self.aMatrixModel.ravel().tolist()
        outDict['modelNormalization'] = self.modelNormalization
        outDict['fitParamCovMatrix'] = self.fitParamCovMatrix.ravel().tolist()
        outDict['usedPixels'] = self.usedPixels.ravel().tolist()
        outDict['fitParams'] = self.fitParams
        outDict['fitParamErrors'] = self.fitParamErrors
        outDict['fitChi2'] = self.fitChi2
        outDict['fitReducedChi2'] = self.fitReducedChi2
        outDict['ath'] = self.ath.ravel().tolist()
        outDict['athMinusBeta'] = self.athMinusBeta.ravel().tolist()
        outDict['aN'] = self.aN.ravel().tolist()
        outDict['aS'] = self.aS.ravel().tolist()
        outDict['aE'] = self.aE.ravel().tolist()
        outDict['aW'] = self.aW.ravel().tolist()
 
        return outDict
 
    @classmethod
    def fromTable(cls, tableList):
        """Construct calibration from a list of tables.
 
        This method uses the `fromDict` method to create the
        calibration, after constructing an appropriate dictionary from
        the input tables.
 
        Parameters
        ----------
        tableList : `list` [`astropy.table.Table`]
            List of tables to use to construct the brighter-fatter
            calibration.
 
        Returns
        -------
        calib : `lsst.ip.isr.BrighterFatterKernel`
            The calibration defined in the tables.
        """
        record = tableList[0][0]
 
        metadata = record.meta
        calibVersion = metadata['BF_DISTORTION_MATRIX_VERSION']
 
        # Initialize inDict with all expected keys
        # and empty values of the corresponding type
        inDict = dict()
        inDict['metadata'] = metadata
        inDict['ampNames'] = record['AMPNAMES']
        inDict['inputRange'] = record['INPUT_RANGE']
        inDict['fitRange'] = record['FIT_RANGE']
        inDict['badAmps'] = record['BAD_AMPS']
        inDict['fitParamNames'] = record['FIT_PARAM_NAMES']
        inDict['freeFitParamNames'] = record['FREE_FIT_PARAM_NAMES']
        inDict['gain'] = {str(n): v for n, v in zip(record['AMPNAMES'], record['GAIN'])}
        inDict['aMatrix'] = record['A_MATRIX']
        inDict['aMatrixSigma'] = record['A_MATRIX_SIGMA']
        inDict['aMatrixModel'] = record['A_MATRIX_MODEL']
        inDict['aMatrixSum'] = record['A_MATRIX_SUM']
        inDict['aMatrixModelSum'] = record['A_MATRIX_MODEL_SUM']
        inDict['modelNormalization'] = record['MODEL_NORMALIZATION']
        inDict['usedPixels'] = record['USED_PIXELS']
        inDict['fitParams'] = {str(n): v for n, v in zip(record['FIT_PARAM_NAMES'], record['FIT_PARAMS'])}
        inDict['fitParamErrors'] = {str(n): v for n, v in zip(record['FIT_PARAM_NAMES'],
                                                              record['FIT_PARAM_ERRORS'])}
        inDict['fitChi2'] = record['FIT_CHI2']
        inDict['fitReducedChi2'] = record['FIT_REDUCED_CHI2']
        inDict['fitParamCovMatrix'] = record['FIT_PARAM_COV_MATRIX']
        inDict['ath'] = record['ATH']
        inDict['athMinusBeta'] = record['ATH_MINUS_BETA']
        inDict['aN'] = record['A_N']
        inDict['aS'] = record['A_S']
        inDict['aE'] = record['A_E']
        inDict['aW'] = record['A_W']
 
        if not calibVersion == 1.0:
            cls().log.warning("Unkown version found for "
                              f"ElectrostaticBrighterFatterDistortionMatrix: {calibVersion}. ")
 
        # Check for newer versions, but there are none...
 
        return cls.fromDict(inDict)
 
    def toTable(self):
        """Construct a list of tables containing the information in this
        calibration.
 
        The list of tables should create an identical calibration
        after being passed to this class's fromTable method.
 
        Returns
        -------
        tableList : `list` [`lsst.afw.table.Table`]
            List of tables containing the crosstalk calibration
            information.
 
        """
        tableList = []
        self.updateMetadata()
 
        recordDict = {
            'AMPNAMES': self.ampNames,
            'INPUT_RANGE': self.inputRange,
            'FIT_RANGE': self.fitRange,
            'BAD_AMPS': self.badAmps,
            'GAIN': list(self.gain.values()),
            'A_MATRIX': self.aMatrix.ravel(),
            'A_MATRIX_SIGMA': self.aMatrixSigma.ravel(),
            'A_MATRIX_MODEL': self.aMatrixModel.ravel(),
            'A_MATRIX_SUM': self.aMatrixSum,
            'A_MATRIX_MODEL_SUM': self.aMatrixModelSum,
            'MODEL_NORMALIZATION': self.modelNormalization,
            'USED_PIXELS': self.usedPixels.ravel(),
            'FIT_PARAM_NAMES': self.fitParamNames,
            'FREE_FIT_PARAM_NAMES': self.freeFitParamNames,
            'FIT_PARAMS': list(self.fitParams.values()),
            'FIT_PARAM_ERRORS': list(self.fitParamErrors.values()),
            'FIT_CHI2': self.fitChi2,
            'FIT_REDUCED_CHI2': self.fitReducedChi2,
            'FIT_PARAM_COV_MATRIX': self.fitParamCovMatrix.ravel(),
            'ATH': self.ath.ravel(),
            'ATH_MINUS_BETA': self.athMinusBeta.ravel(),
            'A_N': self.aN.ravel(),
            'A_S': self.aS.ravel(),
            'A_E': self.aE.ravel(),
            'A_W': self.aW.ravel(),
        }
 
        catalogList = [recordDict]
        catalog = Table(catalogList)
 
        inMeta = self.getMetadata().toDict()
        outMeta = {k: v for k, v in inMeta.items() if v is not None}
        outMeta.update({k: "" for k, v in inMeta.items() if v is None})
        catalog.meta = outMeta
        tableList.append(catalog)
 
        return tableList
 
 
class CustomFFTConvolution(object):
    """
    A class that performs image convolutions in Fourier space, using pyfftw.
    The constructor takes images as arguments and creates the FFTW plans.
    The convolutions are performed by the __call__ routine.
    This is faster than scipy.signal.fftconvolve, and it saves some transforms
    by allowing the same image to be convolved with several kernels.
    pyfftw does not accommodate float32 images, so everything
    should be double precision.
 
    Code adaped from :
    https://stackoverflow.com/questions/14786920/convolution-of-two-three-dimensional-arrays-with-padding-on-one-side-too-slow
    Code posted by Henry Gomersal
    """
    def __init__(self, im, kernel, threads=1):
        # Compute the minimum size of the convolution result.
        shape = (np.array(im.shape) + np.array(kernel.shape)) - 1
 
        # Find the next larger "fast size" for FFT computation.
        # This can provide a significant speedup.
        shape = np.array([pyfftw.next_fast_len(s) for s in shape])
 
        # Create FFTW building objects for the image and kernel.
        self.fftImageObj = pyfftw.builders.rfftn(im, s=shape, threads=threads)
        self.fftKernelObj = pyfftw.builders.rfftn(kernel, s=shape, threads=threads)
        self.ifftObj = pyfftw.builders.irfftn(
            self.fftImageObj.get_output_array(),
            s=shape,
            threads=threads,
        )
 
    def __call__(self, im, kernels):
        """
        Perform the convolution and trim the result to the
        size of the input image. If kernels is a list, then
        the routine returns a list of corresponding
        convolutions.
        """
        # Handle both a list of kernels and a single kernel.
        ks = [kernels] if type(kernels) is not list else kernels
        convolutions = []
        for k in ks:
            # Transform the image and the kernel using FFT.
            tim = self.fftImageObj(im)
            tk = self.fftKernelObj(k)
 
            # Multiply in Fourier space and perform the inverse FFT.
            convolution = self.ifftObj(tim * tk)
            # Trim the result to match the input image size, following
            # the 'same' policy of scipy.signal.fftconvolve.
            oy = k.shape[0] // 2
            ox = k.shape[1] // 2
            convolutions.append(convolution[oy:oy + im.shape[0], ox:ox + im.shape[1]].copy())
        # Return a single convolution if kernels was
        # not a list, otherwise return the list.
        return convolutions[0] if type(kernels) is not list else convolutions
 
 
def electrostaticBrighterFatterCorrection(exposure, electroBfDistortionMatrix, applyGain, gains=None):
    """
    Evaluates the correction of CCD images affected by the
    brighter-fatter effect, as described in
    https://arxiv.org/abs/2301.03274. Requires as input the result of
    an electrostatic fit to flat covariance data (or any other
    determination of pixel boundary shifts under the influence of a
    single electron).
 
    The filename refers to an input tuple that contains the
    boundary shifts for one electron. This file is produced by an
    electrostatic fit to data extracted from flat-field statistics,
    implemented in https://gitlab.in2p3.fr/astier/bfptc/tools/fit_cov.py.
    """
 
    # Use the symmetrize function to fill the four quadrants for each kernel
    r = electroBfDistortionMatrix.fitRange - 1
    aN = electroBfDistortionMatrix.aN
    aS = electroBfDistortionMatrix.aS
    aE = electroBfDistortionMatrix.aE
    aW = electroBfDistortionMatrix.aW
 
    # Initialize kN and kE arrays
    kN = np.zeros((2 * r + 1, 2 * r + 1))
    kE = np.zeros_like(kN)
 
    # Fill in the 4 quadrants for kN
    kN[r:, r:] = aN  # Quadrant 1 (bottom-right)
    kN[:r+1, r:] = np.flipud(aN)  # Quadrant 2 (top-right)
    kN[r:, :r+1] = np.fliplr(aS)  # Quadrant 3 (bottom-left)
    kN[:r+1, :r+1] = np.flipud(np.fliplr(aS))  # Quadrant 4 (top-left)
 
    # Fill in the 4 quadrants for kE
    kE[r:, r:] = aE  # Quadrant 1 (bottom-right)
    kE[:r+1, r:] = np.flipud(aW)  # Quadrant 2 (top-right)
    kE[r:, :r+1] = np.fliplr(aE)  # Quadrant 3 (bottom-left)
    kE[:r+1, :r+1] = np.flipud(np.fliplr(aW))  # Quadrant 4 (top-left)
 
    # Tweak the edges so that the sum rule applies.
    kN[:, 0] = -kN[:, -1]
    kE[0, :] = -kE[-1, :]
 
    # We use the normalization of Guyonnet et al. (2015),
    # which is compatible with the way the input file is produced.
    # The factor 1/2 is due to the fact that the charge distribution at the end
    # is twice the average, and the second 1/2 is due to
    # charge interpolation.
    kN *= 0.25
    kE *= 0.25
 
    # Indeed, i and j in the tuple refer to serial and parallel directions.
    # In most Python code, the image reads im[j, i], so we transpose:
    kN = kN.T
    kE = kE.T
 
    # Get the image to perform the correction
    image = exposure.getMaskedImage().getImage()
 
    # The image needs to be units of electrons/holes
    with gainContext(exposure, image, applyGain, gains):
        # Computes the correction and returns the "delta_image",
        # which should be subtracted from "im" in order to undo the BF effect.
        # The input image should be expressed in electrons
        im = image.getArray().copy()
        convolver = CustomFFTConvolution(im, kN)
        convolutions = convolver(im, [kN, kE])
 
        # The convolutions contain the boundary shifts (in pixel size units)
        # for [horizontal, vertical] boundaries.
        # We now compute the charge to move around.
        delta = np.zeros_like(im)
        boundaryCharge = np.zeros_like(im)
 
        # Horizontal boundaries (parallel direction).
        # We could use a more elaborate interpolator for estimating the
        # charge on the boundary.
        boundaryCharge[:-1, :] = im[1:, :] + im[:-1, :]
        # boundaryCharge[1:-2,:] = (9./8.)*(I[2:-1,:]+I[1:-2,:] -
        # (1./8.)*(I[0:-3,:]+I[3,:])
 
        # The charge to move around is the
        # product of the boundary shift (in pixel size units) times the
        # charge on the boundary (in charge per pixel unit).
        dq = boundaryCharge * convolutions[0]
        delta += dq
 
        # What is gained by a pixel is lost by its neighbor (the right one).
        delta[1:, :] -= dq[:-1, :]
 
        # Vertical boundaries.
        boundaryCharge = np.zeros_like(im)  # Reset to zero.
        boundaryCharge[:, :-1] = im[:, 1:] + im[:, :-1]
        dq = boundaryCharge * convolutions[1]
        delta += dq
 
        # What is gained by a pixel is lost by its neighbor.
        delta[:, 1:] -= dq[:, :-1]
 
        # TODO: One might check that delta.sum() ~ 0 (charge conservation).
 
        # Apply the correction to the original image
        exposure.image.array -= delta
 
    return exposure