Source code for ncempy.algo.radial_profile

Module to calculate radial profiles.

The **settings** dict used for the evaluation has the following entries:

    * lmax_r (int):    Radius used for local maximum detection in [px].
    * lmax_thresh (int):    Threshold for local maximum detection in [intensity].
    * lmax_cinit (int, int):    Initial guess for center in [px].
    * lmax_range (float, float):    r range to cinit used to filter local maxima [dims].
    * ns (int, ...):    Distortion orders to correct.
    * fit_rrange (float, float):    r range used to fit radial profile [dims].
    * back_xs (float, ...):    Fixpoints for fitting background [dims].
    * back_xswidth (float):    Range around fixpoints to take into account [dims].
    * back_init (float):    Initial guess for subtracting background.
    * fit_funcs (str, ...):    List of functions to model radial profile.
    * fit_init (float, ...):    Initial guess for fitting.

    optional (set to None to use defaults):

        * plt_imgminmax (float, float):    Relative range to plot the img.
        * rad_rmax (float):    Max r to be used in radial profile [dims].
        * rad_dr (float): Stepsize for r-axis in radial profile [dims].
        * rad_sigma (float): Sigma for Gaussian used as kernel density estimator [dims].
        * mask (np.ndarray): Binary image as img, 0 for pixels to exclude.
        * fit_maxfev (int): Maxfev forwarded to scipy optimize.


import copy
import numpy as np
import scipy.ndimage.filters
import scipy.interpolate

import ncempy.algo.math
import ncempy.algo.distortion
import ncempy.viz

[docs]def calc_polarcoords(center, dims, ns=None, dists=None): """Calculate the polar coordinates for an image of given shape. Center is assumed to be in real coordinates (if not just fake the dims). Distortions are corrected if ns and corresponding dists are given. Parameters ---------- center : np.ndarray/tuple Center of polar coordinate system. dims : tuple Tuple of dimensions. ns : tuple List of distortion orders to correct for. dists : np.ndarray Parameters for distortions. Returns -------- : tuple Polar coordinates (r, theta) as two np.ndarrays with same dimensions as original image. """ # check input try: # check center center = np.array(center) center = np.reshape(center, 2) # check if enough dims available assert (len(dims) >= 2) assert (len(dims[0]) == 3) # check orders if ns is not None: assert (len(ns) >= 1) # check dists assert (dists.shape[0] == len(ns) * 2 + 1) except TypeError: raise TypeError('Something wrong with the input!') # create coordinates xx, yy = np.meshgrid(dims[0][0], dims[1][0], indexing='ij') # calculate polar coordinate system rs = np.sqrt(np.square(xx - center[0]) + np.square(yy - center[1])) thes = np.arctan2(yy - center[1], xx - center[0]) # correct for distortions if ns is not None: for i in range(len(ns)): rs /= ncempy.algo.distortion.rad_dis(thes, dists[i * 2 + 1], dists[i * 2 + 2], ns[i]) return rs, thes
[docs]def correct_distortion(img, dims, center, ns, dists): """Give corrected version of diffraction pattern with respect to distortions. Parameters ----------- img : np.ndarray Image (2D). dims : tuple Dimensions tuple. center : np.ndarray or tuple Center to be used. ns : tuple List of distortion orders. dists : np.ndarray Distortion parameters. Returns -------- : np.ndarray Corrected diffraction pattern. """ # check input try: assert (isinstance(img, np.ndarray)) # check center center = np.array(center) center = np.reshape(center, 2) # check if enough dims available assert (len(dims) >= 2) assert (len(dims[0]) == 3) # check orders assert (len(ns) >= 1) # check dists assert (dists.shape[0] == len(ns) * 2 + 1) except TypeError: raise TypeError('Something wrong with the input!') rs, thes = calc_polarcoords(center, dims) # anti distort for i in range(len(ns)): rs *= ncempy.algo.distortion.rad_dis(thes, dists[i * 2 + 1], dists[i * 2 + 2], ns[i]) dis_xx = rs * np.cos(thes) + center[0] dis_yy = rs * np.sin(thes) + center[1] f = scipy.interpolate.RectBivariateSpline(dims[0][0], dims[1][0], img) return f.ev(dis_xx, dis_yy)
[docs]def calc_radialprofile(img, rs, rMax, dr, rsigma, mask=None): """Calculate the radial profile using Gaussian kernel density estimator. It is suggested to use an rMax such that all directions are still in the image, otherwise outer areas will contribute differently and lead to different signal-to-noise. A value of dr corresponding to 1/10 px and rsigma corresponding to 1 px are good parameters. Parameters ----------- img : np.ndarray Image to take radial profile of intensity from. rs : np.ndarray Array containing the radial distance, same shape as img. rMax : float Maximum radial distance used for the radial profile. dr : float Step size for r-axis of radial distance. rsigma : float Sigma for Gaussian used as kernel density estimator. mask : np.ndarray Mask image, values at 0 are excluded from radial profile. Returns -------- : tuple Tuple of radial and intensity axes. """ # check input try: assert (isinstance(img, np.ndarray)) assert (isinstance(rs, np.ndarray)) assert (np.array_equal(img.shape, rs.shape)) rMax = float(rMax) dr = float(dr) rsigma = float(rsigma) if not mask is None: assert (isinstance(mask, np.ndarray)) assert (np.array_equal(img.shape, mask.shape)) except: raise TypeError('Something wrong with input.') # process mask if mask is None: mask = np.ones(img.shape) mask = mask.astype('float64') mask[np.where(mask > 0)] = 1 mask[np.where(mask == 0)] = np.nan # prepare radial axis for hist rBins = np.arange(0, rMax, dr) radialIndex = np.round(rs / dr) + 1 rBinMax = len(rBins) sel = (radialIndex <= rBinMax) sel = np.logical_and(sel, np.logical_not(np.isnan(mask))) signal = np.histogram(rs[sel], rBins, weights=img[sel]) count = np.histogram(rs[sel], rBins, weights=np.ones(img[sel].shape)) signal_sm = scipy.ndimage.filters.gaussian_filter1d(signal[0], rsigma / dr) count_sm = scipy.ndimage.filters.gaussian_filter1d(count[0], rsigma / dr) # masked regions lead to 0 in count_sm, divide produces nans, just ignore the warning old_err_state = np.seterr(divide='ignore', invalid='ignore') signal_sm = np.divide(signal_sm, count_sm) np.seterr(**old_err_state) return rBins[:-1], signal_sm
[docs]def residuals_fit(param, r, intens, funcs): """Residual function to fit radial profile with flexibility. The arguments for the single functions are automatically cut from the fit parameters. Parameters ----------- param : np.ndarray Fit parameters. r : np.ndarray r-axis. intens : np.ndarray Intensity-axis. funcs : tuple List of functions to include. Returns -------- : np.ndarray Residuals. """ return intens - ncempy.algo.math.sum_functions(r, funcs, param)
[docs]def fit_radialprofile(r, intens, funcs, init_guess, maxfev=None): """Fit the radial profile. Convenience wrapper for fitting. Parameters ----------- r : np.ndarray r-axis of radial profile. intens : np.ndarray Intensity-axis of radial profile. funcs : tuple List of functions. init_guess : np.ndarray Initial guess for parameters of functions in funcs. maxfev : int maximum function calls (for scipy.optimize.leastsq) Returns -------- : np.ndarray Optimized parameters. """ try: # check data assert (isinstance(r, np.ndarray)) assert (isinstance(intens, np.ndarray)) assert (np.array_equal(r.shape, intens.shape)) # funcs and params assert (len(funcs) >= 1) for i in range(len(funcs)): assert (funcs[i] in ncempy.algo.math.lkp_funcs) init_guess = np.array(init_guess) init_guess = np.reshape(init_guess, sum(map(lambda x: ncempy.algo.math.lkp_funcs[x][1], funcs))) except: raise TypeError('Something wrong with the input!') if maxfev is None: maxfev = 1000 popt, flag = scipy.optimize.leastsq(residuals_fit, init_guess, args=(r, intens, funcs), maxfev=maxfev) if flag not in [1, 2, 3, 4]: print('WARNING: fitting of radial profile failed.') return popt
[docs]def run_singleImage(img, dims, settings, show=False): """Evaluate a single ring diffraction pattern with given settings. Parameters ----------- img : np.ndarray Image. dims : tuple Corresponding dim vectors. settings : dict Dict of settings necessary for the evaluation. show : bool Set to directly show plots interactively. Returns -------- (np.ndarray): Optimized parameters of fitting the radial profile according to the settings. """ try: assert (isinstance(img, np.ndarray)) # check if dims availabel assert (len(dims) >= 1) assert (len(dims[0]) == 3) assert (type(settings) is dict) except: raise RuntimeError('Something wrong with the input') # create a copy of settings to decouple mysettings = copy.deepcopy(settings) # get local maxima an turn them into real space coords points = ncempy.algo.local_max.local_max(img, mysettings['lmax_r'], mysettings['lmax_thresh']) points = ncempy.algo.local_max.points_todim(points, dims) # convert center to real space center_init = ncempy.algo.local_max.points_todim(mysettings['lmax_cinit'], dims) # filter to single ring points = ncempy.algo.distortion.filter_ring(points, center_init, mysettings['lmax_range']) if show: if settings['plt_imgminmax'] is None: mysettings['plt_imgminmax'] = (0., 1.) plot = ncempy.viz.plot_points(img, points, vminmax=mysettings['plt_imgminmax'], dims=dims, invert=True, show=show) # optimize center center = ncempy.algo.distortion.optimize_center(points, center_init, verbose=show) # fit distortions points_plr = ncempy.algo.distortion.points_topolar(points, center) dists = ncempy.algo.distortion.optimize_distortion(points_plr, mysettings['ns']) if show: plot = ncempy.viz.plot_distpolar(points_plr, dims, dists, mysettings['ns'], show=show) # calc coordinates in optimized system rs, thes = ncempy.algo.radial_profile.calc_polarcoords(center, dims, mysettings['ns'], dists) if settings['rad_rmax'] is None: mysettings['rad_rmax'] = np.abs(dims[0][0][0] - dims[0][0][1]) * np.min(img.shape) / 2.0 if settings['rad_dr'] is None: mysettings['rad_dr'] = np.abs(dims[0][0][0] - dims[0][0][1]) / 10. if settings['rad_sigma'] is None: mysettings['rad_sigma'] = np.abs(dims[0][0][0] - dims[0][0][1]) # extract radial profile R, I = ncempy.algo.radial_profile.calc_radialprofile(img, rs, mysettings['rad_rmax'], mysettings['rad_dr'], mysettings['rad_sigma'], mask=mysettings['mask']) # cut radial profile sel = (R >= mysettings['fit_rrange'][0]) * (R <= mysettings['fit_rrange'][1]) I = I[sel] R = R[sel] rawRI = np.copy(np.array([R, I]).transpose()) # subtract a power law background fitted to specific points fit_R = np.array([]) fit_I = np.array([]) for xpoint in mysettings['back_xs']: ix = np.where(np.abs(R - xpoint) < mysettings['back_xswidth']) fit_R = np.append(fit_R, R[ix]) fit_I = np.append(fit_I, I[ix]) funcs_back = ['const', 'powlaw'] res_back = ncempy.algo.radial_profile.fit_radialprofile(fit_R, fit_I, funcs_back, mysettings['back_init'], maxfev=1000) if show: plot = ncempy.viz.plot_fit(R, I, dims, funcs_back, res_back, show=show) I = I - ncempy.algo.math.sum_functions(R, funcs_back, res_back) # fit res = ncempy.algo.radial_profile.fit_radialprofile(R, I, mysettings['fit_funcs'], mysettings['fit_init'], maxfev=mysettings['fit_maxfev']) if show: plot = ncempy.viz.plot_fit(R, I, dims, mysettings['fit_funcs'], res, show=show) return np.array([R, I]).transpose(), res, center, dists, rawRI, res_back, mysettings