Source code for qpimage.bg_estimate

import lmfit
import numpy as np


#: valid values for keyword argument `fit_offset` in :func:`estimate`
VALID_FIT_OFFSETS = ["fit", "gauss", "mean", "mode"]
#: valid values for keyword argument `fit_profile` in :func:`estimate`
VALID_FIT_PROFILES = ["offset", "poly2o", "tilt"]


[docs]def estimate(data, fit_offset="mean", fit_profile="tilt", border_px=0, from_mask=None, ret_mask=False): """Estimate the background value of an image Parameters ---------- data: np.ndarray Data from which to compute the background value fit_profile: str The type of background profile to fit: - "offset": offset only - "poly2o": 2D 2nd order polynomial with mixed terms - "tilt": 2D linear tilt with offset (default) fit_offset: str The method for computing the profile offset - "fit": offset as fitting parameter - "gauss": center of a gaussian fit - "mean": simple average - "mode": mode (see `qpimage.bg_estimate.mode`) border_px: float Assume that a frame of `border_px` pixels around the image is background. from_mask: boolean np.ndarray or None Use a boolean array to define the background area. The boolean mask must have the same shape as the input data. `True` elements are used for background estimation. ret_mask: bool Return the boolean mask used to compute the background. Notes ----- If both `border_px` and `from_mask` are given, the intersection of the two is used, i.e. the positions where both, the frame mask and `from_mask`, are `True`. """ if fit_profile not in VALID_FIT_PROFILES: msg = "`fit_profile` must be one of {}, got '{}'".format( VALID_FIT_PROFILES, fit_profile) raise ValueError(msg) if fit_offset not in VALID_FIT_OFFSETS: msg = "`fit_offset` must be one of {}, got '{}'".format( VALID_FIT_OFFSETS, fit_offset) raise ValueError(msg) # initial mask image if from_mask is not None: assert isinstance(from_mask, np.ndarray) mask = from_mask.copy() else: mask = np.ones_like(data, dtype=bool) # multiply with border mask image (intersection) if border_px > 0: border_px = int(np.round(border_px)) mask_px = np.zeros_like(mask) mask_px[:border_px, :] = True mask_px[-border_px:, :] = True mask_px[:, :border_px] = True mask_px[:, -border_px:] = True # intersection np.logical_and(mask, mask_px, out=mask) # compute background image if fit_profile == "tilt": bgimg = profile_tilt(data, mask) elif fit_profile == "poly2o": bgimg = profile_poly2o(data, mask) else: bgimg = np.zeros_like(data, dtype=float) # add offsets if fit_offset == "fit": if fit_profile == "offset": msg = "`fit_offset=='fit'` only valid when `fit_profile!='offset`" raise ValueError(msg) # nothing else to do here, using offset from fit elif fit_offset == "gauss": bgimg += offset_gaussian((data - bgimg)[mask]) elif fit_offset == "mean": bgimg += np.mean((data - bgimg)[mask]) elif fit_offset == "mode": bgimg += offset_mode((data - bgimg)[mask]) if ret_mask: ret = (bgimg, mask) else: ret = bgimg return ret
[docs]def offset_gaussian(data): """Fit a gaussian model to `data` and return its center""" nbins = 2 * int(np.ceil(np.sqrt(data.size))) mind, maxd = data.min(), data.max() drange = (mind - (maxd - mind) / 2, maxd + (maxd - mind) / 2) histo = np.histogram(data, nbins, density=True, range=drange) dx = abs(histo[1][1] - histo[1][2]) / 2 hx = histo[1][1:] - dx hy = histo[0] # fit gaussian gauss = lmfit.models.GaussianModel() pars = gauss.guess(hy, x=hx) out = gauss.fit(hy, pars, x=hx) return out.params["center"]
[docs]def offset_mode(data): """Compute Mode using a histogram with `sqrt(data.size)` bins""" nbins = int(np.ceil(np.sqrt(data.size))) mind, maxd = data.min(), data.max() histo = np.histogram(data, nbins, density=True, range=(mind, maxd)) dx = abs(histo[1][1] - histo[1][2]) / 2 hx = histo[1][1:] - dx hy = histo[0] idmax = np.argmax(hy) return hx[idmax]
[docs]def profile_tilt(data, mask): """Fit a 2D tilt to `data[mask]`""" params = lmfit.Parameters() params.add(name="mx", value=0) params.add(name="my", value=0) params.add(name="off", value=np.average(data[mask])) fr = lmfit.minimize(tilt_residual, params, args=(data, mask)) bg = tilt_model(fr.params, data.shape) return bg
[docs]def profile_poly2o(data, mask): """Fit a 2D 2nd order polynomial to `data[mask]`""" # lmfit params = lmfit.Parameters() params.add(name="mx", value=0) params.add(name="my", value=0) params.add(name="mxy", value=0) params.add(name="ax", value=0) params.add(name="ay", value=0) params.add(name="off", value=np.average(data[mask])) fr = lmfit.minimize(poly2o_residual, params, args=(data, mask)) bg = poly2o_model(fr.params, data.shape) return bg
[docs]def poly2o_model(params, shape): """lmfit 2nd order polynomial model""" mx = params["mx"].value my = params["my"].value mxy = params["mxy"].value ax = params["ax"].value ay = params["ay"].value off = params["off"].value bg = np.zeros(shape, dtype=float) + off x = np.arange(bg.shape[0]) - bg.shape[0] // 2 y = np.arange(bg.shape[1]) - bg.shape[1] // 2 x = x.reshape(-1, 1) y = y.reshape(1, -1) bg += ax * x**2 + ay * y**2 + mx * x + my * y + mxy * x * y return bg
[docs]def poly2o_residual(params, data, mask): """lmfit 2nd order polynomial residuals""" bg = poly2o_model(params, shape=data.shape) res = (data - bg)[mask] return res.flatten()
[docs]def tilt_model(params, shape): """lmfit tilt model""" mx = params["mx"].value my = params["my"].value off = params["off"].value bg = np.zeros(shape, dtype=float) + off x = np.arange(bg.shape[0]) - bg.shape[0] // 2 y = np.arange(bg.shape[1]) - bg.shape[1] // 2 x = x.reshape(-1, 1) y = y.reshape(1, -1) bg += mx * x + my * y return bg
[docs]def tilt_residual(params, data, mask): """lmfit tilt residuals""" bg = tilt_model(params, shape=data.shape) res = (data - bg)[mask] return res.flatten()