"""
Background substraction is a necessary component of reflectometry reduction,
where the background scattering is removed from the reflected intensity.
Herein are some functions to enable that for a two-dimensional detector image,
as well as simple dataclasses in which we can store some information relating to
the background subtraction, and any fitting that we might have carried out.
"""
from dataclasses import dataclass
from typing import Callable, List
import numpy as np
from scipy.stats import norm
from scipy.optimize import curve_fit
from .region import Region
from .image import Image
[docs]@dataclass
class FitInfo:
"""
A simple dataclass in which we can store data relating to the quality of a
fit.
"""
popt: np.ndarray
pcov: np.ndarray
fit_function: Callable
[docs]@dataclass
class BkgSubInfo:
"""
A simple data class in which we can store information relating to a
background subtraction.
"""
bkg: float
bkg_e: float
bkg_sub_function: Callable
fit_info: FitInfo = None
[docs]def roi_subtraction(image, list_of_regions: List[Region]):
"""
Carry out background subtraction by taking a series of rectangular regions
of interested (ROIs) as being fair Poissonian measurements of the
background.
Args:
image:
The islatu.image.Image object from which we should subtract
background from.
list_of_regions:
A list of instances of the Regions class corresponding to background
regions.
"""
# We're going to need to count all intensity in all the background, as well
# as the number of pixels used in our measurement of the background.
sum_of_bkg_areas = 0
total_num_pixels = 0
# Make sure we've been given multiple regions. If not, np: make a list.
if isinstance(list_of_regions, Region):
list_of_regions = [list_of_regions]
# Add up all the intensity in all the pixels.
for region in list_of_regions:
# Now add the total intensity in this particular background region to
# the intensity measured in all the background regions so far.
sum_of_bkg_areas += np.sum(
image.array_original[
int(region.x_start):int(region.x_end),
int(region.y_start):int(region.y_end)
]
)
# Add the number of pixels in this background ROI to the total number of
# pixels used to compute the background measurement overall.
total_num_pixels += region.num_pixels
# Now Poisson stats can be abused to only calculate a single sqrt.
err_of_bkg_areas = np.sqrt(sum_of_bkg_areas)
if err_of_bkg_areas == 0:
err_of_bkg_areas = 1
# Get the per pixel background mean and stddev.
bkg_per_pixel = sum_of_bkg_areas / total_num_pixels
bkg_error_per_pixel = err_of_bkg_areas / total_num_pixels
# Expose the calculated background and background_error per pixel.
return BkgSubInfo(bkg_per_pixel, bkg_error_per_pixel, roi_subtraction)
[docs]def univariate_normal(data, mean, sigma, offset, factor):
"""
Produce a univariate normal distribution.
Args:
data (:py:attr:`array_like`): Abscissa data.
mean (:py:attr:`float`): Mean (horizontal).
sigma (:py:attr:`float`): Variance (horizontal).
offset (:py:attr:`float`): Offset from the 0 for the ordinate, this is
the background level.
factor (:py:attr:`float`): Multiplicative factor for area of normal
distribution.
Returns:
:py:attr:`array_like`: Ordinate data for univariate normal distribution.
"""
# Creation of the bivariate normal distribution
normal = norm(loc=mean, scale=sigma)
return offset + normal.pdf(data).flatten() * factor
[docs]def fit_gaussian_1d(image: Image, params_0=None, bounds=None, axis=0):
"""
Fit a one-dimensional Gaussian function with some ordinate offset to an
image with uncertainty. This is achieved by averaging in a given ``axis``
before performing the fit. Return the results, and index of the offset.
Args:
image:
The islatu image object to fit.
params_0 (:py:attr:`list`, optional):
An initial guess at the parameters. Defaults to values based on the
image.
bounds (:py:attr:`list` of :py:attr:`tuple`, optional):
Bounds for the fitting. Defaults to values based on the image.
axis (:py:attr:`int`):
The dimension along which the averaging will be performed.
Returns:
:py:attr:`tuple`: Containing:
- :py:attr:`array_like`: The results (with uncertainties) for each
of the 6 parameters fit.
- :py:attr:`int`: The index of the offset.
- :py:attr:`None`: As it is not possible to describe the reflected
peak width.
"""
arr, arr_e = image.array, image.array_e
ordinate = arr.mean(axis=axis)
# Now we can generate an array of errors.
ordinate_e = np.sqrt(np.mean(arr_e**2, axis=axis))
# Setting default values.
if params_0 is None:
# Now we generate the initial values for our Gaussian fit.
# These values are crucial – as this is a high dimensional fitting
# problem, it is likely that we'll get stuck in a local minimum if these
# aren't good.
# Guess that the Gaussian mean is at the most intense mean pixel value.
mean0 = np.argmax(ordinate)
# Guess that the standard deviation is a single pixel.
sdev0 = 1
# Guess that the background (offset) is the median pixel value.
offset0 = np.median(ordinate)
# Guess that the scale is equal to the largest recorded value.
scale0 = arr.max()
params_0 = [mean0, sdev0, offset0, scale0]
if bounds is None:
bounds = ([0, 0, 0, 0],
[ordinate.shape[0], ordinate.shape[0], scale0, scale0 * 10])
# Perform the fitting.
fit_popt_pcov = curve_fit(
univariate_normal,
np.arange(0, ordinate.shape[0], 1), ordinate, bounds=bounds,
sigma=ordinate_e, p0=params_0, maxfev=2000 * (len(params_0) + 1))
fit_info = FitInfo(fit_popt_pcov[0], fit_popt_pcov[1], univariate_normal)
# Determine uncertainty from covarience matrix.
# Note: the stddev of the fit Gaussian can be accessed via popt[1].
p_sigma = np.sqrt(np.diag(fit_info.pcov))
return BkgSubInfo(fit_info.popt[2], p_sigma[2], fit_gaussian_1d, fit_info)