lsst.ip.isr.brighterFatterKernel Namespace Reference

Classes
class	BrighterFatterKernel

Functions
	brighterFatterCorrection (exposure, kernel, maxIter, threshold, applyGain, gains=None)
	transferFlux (cFunc, fStep, correctionMode=True)
	fluxConservingBrighterFatterCorrection (exposure, kernel, maxIter, threshold, applyGain, gains=None, correctionMode=True)

Detailed Description

Brighter Fatter Kernel calibration definition.

Function Documentation

brighterFatterCorrection()

lsst.ip.isr.brighterFatterKernel.brighterFatterCorrection	(		exposure,
			kernel,
			maxIter,
			threshold,
			applyGain,
			gains = None )

Apply brighter fatter correction in place for the image.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to have brighter-fatter correction applied.  Modified
    by this method.
kernel : `numpy.ndarray`
    Brighter-fatter kernel to apply.
maxIter : scalar
    Number of correction iterations to run.
threshold : scalar
    Convergence threshold in terms of the sum of absolute
    deviations between an iteration and the previous one.
applyGain : `Bool`
    If True, then the exposure values are scaled by the gain prior
    to correction.
gains : `dict` [`str`, `float`]
    A dictionary, keyed by amplifier name, of the gains to use.
    If gains is None, the nominal gains in the amplifier object are used.

Returns
-------
diff : `float`
    Final difference between iterations achieved in correction.
iteration : `int`
    Number of iterations used to calculate correction.

Notes
-----
This correction takes a kernel that has been derived from flat
field images to redistribute the charge.  The gradient of the
kernel is the deflection field due to the accumulated charge.

Given the original image I(x) and the kernel K(x) we can compute
the corrected image Ic(x) using the following equation:

Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))

To evaluate the derivative term we expand it as follows:

0.5 * ( d/dx(I(x))*d/dx(int(dy*K(x-y)*I(y)))
    + I(x)*d^2/dx^2(int(dy* K(x-y)*I(y))) )

Because we use the measured counts instead of the incident counts
we apply the correction iteratively to reconstruct the original
counts and the correction.  We stop iterating when the summed
difference between the current corrected image and the one from
the previous iteration is below the threshold.  We do not require
convergence because the number of iterations is too large a
computational cost.  How we define the threshold still needs to be
evaluated, the current default was shown to work reasonably well
on a small set of images.  For more information on the method see
DocuShare Document-19407.

The edges as defined by the kernel are not corrected because they
have spurious values due to the convolution.

Definition at line 589 of file brighterFatterKernel.py.

def brighterFatterCorrection(exposure, kernel, maxIter, threshold, applyGain, gains=None):
    """Apply brighter fatter correction in place for the image.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure to have brighter-fatter correction applied.  Modified
        by this method.
    kernel : `numpy.ndarray`
        Brighter-fatter kernel to apply.
    maxIter : scalar
        Number of correction iterations to run.
    threshold : scalar
        Convergence threshold in terms of the sum of absolute
        deviations between an iteration and the previous one.
    applyGain : `Bool`
        If True, then the exposure values are scaled by the gain prior
        to correction.
    gains : `dict` [`str`, `float`]
        A dictionary, keyed by amplifier name, of the gains to use.
        If gains is None, the nominal gains in the amplifier object are used.
 
    Returns
    -------
    diff : `float`
        Final difference between iterations achieved in correction.
    iteration : `int`
        Number of iterations used to calculate correction.
 
    Notes
    -----
    This correction takes a kernel that has been derived from flat
    field images to redistribute the charge.  The gradient of the
    kernel is the deflection field due to the accumulated charge.
 
    Given the original image I(x) and the kernel K(x) we can compute
    the corrected image Ic(x) using the following equation:
 
    Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
 
    To evaluate the derivative term we expand it as follows:
 
    0.5 * ( d/dx(I(x))*d/dx(int(dy*K(x-y)*I(y)))
        + I(x)*d^2/dx^2(int(dy* K(x-y)*I(y))) )
 
    Because we use the measured counts instead of the incident counts
    we apply the correction iteratively to reconstruct the original
    counts and the correction.  We stop iterating when the summed
    difference between the current corrected image and the one from
    the previous iteration is below the threshold.  We do not require
    convergence because the number of iterations is too large a
    computational cost.  How we define the threshold still needs to be
    evaluated, the current default was shown to work reasonably well
    on a small set of images.  For more information on the method see
    DocuShare Document-19407.
 
    The edges as defined by the kernel are not corrected because they
    have spurious values due to the convolution.
    """
    image = exposure.getMaskedImage().getImage()
 
    # The image needs to be units of electrons/holes
    with gainContext(exposure, image, applyGain, gains):
 
        kLx = np.shape(kernel)[0]
        kLy = np.shape(kernel)[1]
        kernelImage = afwImage.ImageD(kLx, kLy)
        kernelImage.getArray()[:, :] = kernel
        tempImage = afwImage.ImageD(image, deep=True)
 
        nanIndex = np.isnan(tempImage.getArray())
        tempImage.getArray()[nanIndex] = 0.
 
        corr = np.zeros(image.array.shape, dtype=np.float64)
        prev_image = np.zeros(image.array.shape, dtype=np.float64)
 
        # Define boundary by convolution region.  The region that the
        # correction will be calculated for is one fewer in each dimension
        # because of the second derivative terms.
        # NOTE: these need to use integer math, as we're using start:end as
        # numpy index ranges.
        startX = kLx//2
        endX = -kLx//2
        startY = kLy//2
        endY = -kLy//2
 
        for iteration in range(maxIter):
 
            outArray = scipy.signal.convolve(
                tempImage.array,
                kernelImage.array,
                mode="same",
                method="fft",
            )
            tmpArray = tempImage.getArray()
 
            with np.errstate(invalid="ignore", over="ignore"):
                # First derivative term
                gradTmp = np.gradient(tmpArray[startY:endY, startX:endX])
                gradOut = np.gradient(outArray[startY:endY, startX:endX])
                first = (gradTmp[0]*gradOut[0] + gradTmp[1]*gradOut[1])[1:-1, 1:-1]
 
                # Second derivative term
                diffOut20 = np.diff(outArray, 2, 0)[startY:endY, startX + 1:endX - 1]
                diffOut21 = np.diff(outArray, 2, 1)[startY + 1:endY - 1, startX:endX]
                second = tmpArray[startY + 1:endY - 1, startX + 1:endX - 1]*(diffOut20 + diffOut21)
 
                corr[startY + 1:endY - 1, startX + 1:endX - 1] = 0.5*(first + second)
 
                tmpArray[:, :] = image.getArray()[:, :]
                tmpArray[nanIndex] = 0.
                tmpArray[startY:endY, startX:endX] += corr[startY:endY, startX:endX]
 
            if iteration > 0:
                diff = np.sum(np.abs(prev_image - tmpArray), dtype=np.float64)
 
                if diff < threshold:
                    break
                prev_image[:, :] = tmpArray[:, :]
 
        image.getArray()[startY + 1:endY - 1, startX + 1:endX - 1] += \
            corr[startY + 1:endY - 1, startX + 1:endX - 1]
 
    return diff, iteration
 
 

fluxConservingBrighterFatterCorrection()

lsst.ip.isr.brighterFatterKernel.fluxConservingBrighterFatterCorrection	(		exposure,
			kernel,
			maxIter,
			threshold,
			applyGain,
			gains = None,
			correctionMode = True )

Apply brighter fatter correction in place for the image.

This version presents a modified version of the algorithm
found in ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``
which conserves the image flux, resulting in improved
correction of the cores of stars. The convolution has also been
modified to mitigate edge effects.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure to have brighter-fatter correction applied.  Modified
    by this method.
kernel : `np.ndarray`
    Brighter-fatter kernel to apply.
maxIter : scalar
    Number of correction iterations to run.
threshold : scalar
    Convergence threshold in terms of the sum of absolute
    deviations between an iteration and the previous one.
applyGain : `Bool`
    If True, then the exposure values are scaled by the gain prior
    to correction.
gains : `dict` [`str`, `float`]
    A dictionary, keyed by amplifier name, of the gains to use.
    If gains is None, the nominal gains in the amplifier object are used.
correctionMode : `Bool`
    If True (default) the function applies correction for BFE.  If False,
    the code can instead be used to generate a simulation of BFE (sign
    change in the direction of the effect)

Returns
-------
diff : `float`
    Final difference between iterations achieved in correction.
iteration : `int`
    Number of iterations used to calculate correction.

Notes
-----
Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.

This correction takes a kernel that has been derived from flat
field images to redistribute the charge.  The gradient of the
kernel is the deflection field due to the accumulated charge.

Given the original image I(x) and the kernel K(x) we can compute
the corrected image Ic(x) using the following equation:

Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))

Improved algorithm at this step applies the divergence theorem to
obtain a pixelised correction.

Because we use the measured counts instead of the incident counts
we apply the correction iteratively to reconstruct the original
counts and the correction.  We stop iterating when the summed
difference between the current corrected image and the one from
the previous iteration is below the threshold.  We do not require
convergence because the number of iterations is too large a
computational cost.  How we define the threshold still needs to be
evaluated, the current default was shown to work reasonably well
on a small set of images.

Edges are handled in the convolution by padding.  This is still not
a physical model for the edge, but avoids discontinuity in the correction.

Author of modified version: Lance.Miller@physics.ox.ac.uk
(see DM-38555).

Definition at line 799 of file brighterFatterKernel.py.

                                           gains=None, correctionMode=True):
    """Apply brighter fatter correction in place for the image.
 
    This version presents a modified version of the algorithm
    found in ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``
    which conserves the image flux, resulting in improved
    correction of the cores of stars. The convolution has also been
    modified to mitigate edge effects.
 
    Parameters
    ----------
    exposure : `lsst.afw.image.Exposure`
        Exposure to have brighter-fatter correction applied.  Modified
        by this method.
    kernel : `np.ndarray`
        Brighter-fatter kernel to apply.
    maxIter : scalar
        Number of correction iterations to run.
    threshold : scalar
        Convergence threshold in terms of the sum of absolute
        deviations between an iteration and the previous one.
    applyGain : `Bool`
        If True, then the exposure values are scaled by the gain prior
        to correction.
    gains : `dict` [`str`, `float`]
        A dictionary, keyed by amplifier name, of the gains to use.
        If gains is None, the nominal gains in the amplifier object are used.
    correctionMode : `Bool`
        If True (default) the function applies correction for BFE.  If False,
        the code can instead be used to generate a simulation of BFE (sign
        change in the direction of the effect)
 
    Returns
    -------
    diff : `float`
        Final difference between iterations achieved in correction.
    iteration : `int`
        Number of iterations used to calculate correction.
 
    Notes
    -----
    Modified version of ``lsst.ip.isr.isrFunctions.brighterFatterCorrection``.
 
    This correction takes a kernel that has been derived from flat
    field images to redistribute the charge.  The gradient of the
    kernel is the deflection field due to the accumulated charge.
 
    Given the original image I(x) and the kernel K(x) we can compute
    the corrected image Ic(x) using the following equation:
 
    Ic(x) = I(x) + 0.5*d/dx(I(x)*d/dx(int( dy*K(x-y)*I(y))))
 
    Improved algorithm at this step applies the divergence theorem to
    obtain a pixelised correction.
 
    Because we use the measured counts instead of the incident counts
    we apply the correction iteratively to reconstruct the original
    counts and the correction.  We stop iterating when the summed
    difference between the current corrected image and the one from
    the previous iteration is below the threshold.  We do not require
    convergence because the number of iterations is too large a
    computational cost.  How we define the threshold still needs to be
    evaluated, the current default was shown to work reasonably well
    on a small set of images.
 
    Edges are handled in the convolution by padding.  This is still not
    a physical model for the edge, but avoids discontinuity in the correction.
 
    Author of modified version: Lance.Miller@physics.ox.ac.uk
    (see DM-38555).
    """
    image = exposure.getMaskedImage().getImage()
 
    # The image needs to be units of electrons/holes
    with gainContext(exposure, image, applyGain, gains):
 
        # get kernel and its shape
        kLy, kLx = kernel.shape
        kernelImage = afwImage.ImageD(kLx, kLy)
        kernelImage.getArray()[:, :] = kernel
        tempImage = afwImage.ImageD(image, deep=True)
 
        nanIndex = np.isnan(tempImage.getArray())
        tempImage.getArray()[nanIndex] = 0.
 
        outImage = afwImage.ImageD(image.getDimensions())
        corr = np.zeros(image.array.shape, dtype=np.float64)
        prevImage = np.zeros(image.array.shape, dtype=np.float64)
        convCntrl = afwMath.ConvolutionControl(False, False, 1)
        fixedKernel = afwMath.FixedKernel(kernelImage)
 
        # set the padding amount
        # ensure we pad by an even amount larger than the kernel
        kLy = 2 * ((1+kLy)//2)
        kLx = 2 * ((1+kLx)//2)
 
        # The deflection potential only depends on the gradient of
        # the convolution, so we can subtract the mean, which then
        # allows us to pad the image with zeros and avoid wrap-around effects
        # (although still not handling the image edges with a physical model)
        # This wouldn't be great if there were a strong image gradient.
        imYdimension, imXdimension = tempImage.array.shape
        imean = np.mean(tempImage.getArray()[~nanIndex], dtype=np.float64)
        # subtract mean from image
        tempImage -= imean
        tempImage.array[nanIndex] = 0.0
        padArray = np.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
        outImage = afwImage.ImageD(np.pad(outImage.getArray(), ((0, kLy), (0, kLx))))
        # Convert array to afw image so afwMath.convolve works
        padImage = afwImage.ImageD(padArray.shape[1], padArray.shape[0])
        padImage.array[:] = padArray
 
        for iteration in range(maxIter):
 
            # create deflection potential, convolution of flux with kernel
            # using padded counts array
            afwMath.convolve(outImage, padImage, fixedKernel, convCntrl)
            tmpArray = tempImage.getArray()
            outArray = outImage.getArray()
 
            # trim convolution output back to original shape
            outArray = outArray[:imYdimension, :imXdimension]
 
            # generate the correction array, with correctionMode set as input
            corr[...] = transferFlux(outArray, tmpArray, correctionMode=correctionMode)
 
            # update the arrays for the next iteration
            tmpArray[:, :] = image.getArray()[:, :]
            tmpArray += corr
            tmpArray[nanIndex] = 0.
            # update padded array
            # subtract mean
            tmpArray -= imean
            tempImage.array[nanIndex] = 0.
            padArray = np.pad(tempImage.getArray(), ((0, kLy), (0, kLx)))
 
            if iteration > 0:
                diff = np.sum(np.abs(prevImage - tmpArray), dtype=np.float64)
 
                if diff < threshold:
                    break
                prevImage[:, :] = tmpArray[:, :]
 
        image.getArray()[:] += corr[:]
 
    return diff, iteration

transferFlux()

lsst.ip.isr.brighterFatterKernel.transferFlux	(		cFunc,
			fStep,
			correctionMode = True )

Take the input convolved deflection potential and the flux array
to compute and apply the flux transfer into the correction array.

Parameters
----------
cFunc: `np.array`
    Deflection potential, being the convolution of the flux F with the
    kernel K.
fStep: `np.array`
    The array of flux values which act as the source of the flux transfer.
correctionMode: `bool`
    Defines if applying correction (True) or generating sims (False).

Returns
-------
corr:
    BFE correction array

Definition at line 715 of file brighterFatterKernel.py.

def transferFlux(cFunc, fStep, correctionMode=True):
    """Take the input convolved deflection potential and the flux array
    to compute and apply the flux transfer into the correction array.
 
    Parameters
    ----------
    cFunc: `np.array`
        Deflection potential, being the convolution of the flux F with the
        kernel K.
    fStep: `np.array`
        The array of flux values which act as the source of the flux transfer.
    correctionMode: `bool`
        Defines if applying correction (True) or generating sims (False).
 
    Returns
    -------
    corr:
        BFE correction array
    """
 
    if cFunc.shape != fStep.shape:
        raise RuntimeError(f'transferFlux: array shapes do not match: {cFunc.shape}, {fStep.shape}')
 
    # set the sign of the correction and set its value for the
    # time averaged solution
    if correctionMode:
        # negative sign if applying BFE correction
        factor = -0.5
    else:
        # positive sign if generating BFE simulations
        factor = 0.5
 
    # initialise the BFE correction image to zero
    corr = np.zeros(cFunc.shape, dtype=np.float64)
 
    # Generate a 2D mesh of x,y coordinates
    yDim, xDim = cFunc.shape
    y = np.arange(yDim, dtype=int)
    x = np.arange(xDim, dtype=int)
    xc, yc = np.meshgrid(x, y)
 
    # process each axis in turn
    for ax in [0, 1]:
 
        # gradient of phi on right/upper edge of pixel
        diff = np.diff(cFunc, axis=ax)
 
        # expand array back to full size with zero gradient at the end
        gx = np.zeros(cFunc.shape, dtype=np.float64)
        yDiff, xDiff = diff.shape
        gx[:yDiff, :xDiff] += diff
 
        # select pixels with either positive gradients on the right edge,
        # flux flowing to the right/up
        # or negative gradients, flux flowing to the left/down
        for i, sel in enumerate([gx > 0, gx < 0]):
            xSelPixels = xc[sel]
            ySelPixels = yc[sel]
            # and add the flux into the pixel to the right or top
            # depending on which axis we are handling
            if ax == 0:
                xPix = xSelPixels
                yPix = ySelPixels+1
            else:
                xPix = xSelPixels+1
                yPix = ySelPixels
            # define flux as the either current pixel value or pixel
            # above/right
            # depending on whether positive or negative gradient
            if i == 0:
                # positive gradients, flux flowing to higher coordinate values
                flux = factor * fStep[sel]*gx[sel]
            else:
                # negative gradients, flux flowing to lower coordinate values
                flux = factor * fStep[yPix, xPix]*gx[sel]
            # change the fluxes of the donor and receiving pixels
            # such that flux is conserved
            corr[sel] -= flux
            corr[yPix, xPix] += flux
 
    # return correction array
    return corr