ncempy.algo package

Submodules

ncempy.algo.align module

ncempy.algo.distortion module

Module to handle distortions in diffraction patterns.

ncempy.algo.distortion.filter_ring(points, center, rminmax)[source]

Filter points to be in a certain radial distance range from center.

Parameters:
  • points (np.ndarray) – Candidate points.
  • center (np.ndarray or tuple) – Center position.
  • rminmax (tuple) – Tuple of min and max radial distance.
Returns:

List of filtered points, two column array.

Return type:

np.ndarray

ncempy.algo.distortion.optimize_center(points, center, maxfev=1000, verbose=None)[source]

Optimize the center by minimizing the sum of square deviations from the mean radial distance.

Parameters:
  • points (np.ndarray) – The points to which the optimization is done (x,y coords in org image).
  • center (np.ndarray or tuple) – Initial center guess.
  • maxfev (int) – Max number of iterations forwarded to scipy.optimize.leastsq().
  • verbose (bool) – Set to get verbose output.
Returns:

The optimized center.

Return type:

np.ndarray

ncempy.algo.distortion.optimize_distortion(points, ns, maxfev=1000, verbose=False)[source]

Optimize distortions.

The orders in the list ns are first fitted subsequently and the result is refined in a final fit simultaneously fitting all orders.

Parameters:
  • points (np.ndarray) – Points to optimize to (in polar coords).
  • ns (tuple) – List of orders to correct for.
  • maxfev (int) – Max number of iterations forwarded to scipy.optimize.leastsq().
  • verbose (bool) – Set for verbose output.
Returns:

Optimized parameters according to ns.

Return type:

np.ndarray

ncempy.algo.distortion.points_topolar(points, center)[source]

Convert points to polar coordinate system.

Can be either in pixel or real dim, but should be the same for points and center.

Parameters:
  • points (np.ndarray) – Positions as two column array.
  • center (np.ndarray or tuple) – Origin of the polar coordinate system.
Returns:

Positions in polar coordinate system as two column array (r, theta).

Return type:

np.ndarray

ncempy.algo.distortion.rad_dis(theta, alpha, beta, order=2)[source]

Radial distortion due to ellipticity or higher order distortion.

Relative distortion, to be multiplied with radial distance.

Parameters:
  • theta (np.ndarray) – Angles at which to evaluate. Must be float.
  • alpha (float) – Orientation of major axis.
  • beta (float) – Strength of distortion (beta = (1-r_min/r_max)/(1+r_min/r_max).
  • order (int) – Order of distortion.
Returns:

Distortion factor.

Return type:

np.ndarray

ncempy.algo.distortion.residuals_center(param, data)[source]

Residual function for minimizing the deviations from the mean radial distance.

Parameters:
  • param (np.ndarray) – The center to optimize.
  • data (np.ndarray) – The points in x,y coordinates of the original image.
Returns:

Residuals.

Return type:

np.ndarray

ncempy.algo.distortion.residuals_dis(param, points, ns)[source]

Residual function for distortions.

Parameters:
  • param (np.ndarray) – Parameters for distortion.
  • points (np.ndarray) – Points to fit to.
  • ns (tuple) – List of orders to account for.
Returns:

Residuals.

Return type:

np.ndarray

ncempy.algo.fourier_operations module

ncempy.algo.fourier_operations.bandpass_filter(im, inner, outer, sigma)[source]

Apply a bandpass filter with smoothed edges in Fourier space.

Parameters:
  • im (ndarray, 2D) – The image to apply the bandpass filter to.
  • inner (float) – The lowest spatial frequency allowed in terms of the Nyquist frequency. Should be between 0 and 0.5.
  • outer (float) – The highest spatial frequency allowed in terms of the Nyquist frequency. Should be between 0 and 0.5.
  • sigma (float) – The standard deviation of the Gaussian smoothing to remove hard edges.
Returns:

The image with the bandpass filter applied.

Return type:

ndarray, 2D

ncempy.algo.fourier_operations.rotateImage(im, theta, pad=False)[source]

Use three shears in Fourier space to exactly (and reversibly) rotate an image.

Currently only works for square images.

Parameters:
  • im (ndarray, 2D) – The image to rotate
  • theta (float) – The angle to rotate by in radians
  • pad (bool) – Add padding to the image before rotating
ncempy.algo.fourier_operations.shearImage(im, dim, shear_factor)[source]

Exact shear of an image using the frequency space of the FFT of the image.

Currently only works for square images.

Parameters:
  • im (ndarray, 2D) – The image to shear.
  • dim (int) – The axis to shear along
  • shear_factor (float) – The amount of shear.
Returns:

The sheared image.

Return type:

ndarray

ncempy.algo.fourier_operations.shiftImage(im, shift)[source]

Exact shear of an image using the frequency space of the FFT of the image.

Parameters:
  • im (ndarray, 2D) – The image to shear.
  • shift (tuple of floats) – The number of pixels to shift in each direction
Returns:

The shifted image.

Return type:

ndarray

ncempy.algo.local_max module

Module to find local maxima in an image.

ncempy.algo.local_max.local_max(img, r, thresh)[source]

Find local maxima from comparing dilated and eroded images.

Calculates images with maximum and minimum within given radius. If the difference is larger than the threshold, the original pixel position with max value is detected as local maximum.

Parameters:
  • img (np.ndarray) – Input image.
  • r (int) – Radius for locality.
  • thresh (int or float) – Intensity difference threshold.
Returns:

Array of points.

Return type:

np.ndarray

ncempy.algo.local_max.points_todim(points, dims)[source]

Convert points from px coordinates to real dim.

Points are expected to be array indices for the first two dimensions in dims.

Parameters:
  • points (np.ndarray) – Points to convert.
  • dims (tuple) – Tuple of dimensions.
Returns:

Converted points.

Return type:

np.ndarray

ncempy.algo.math module

Module containing definitions of basic math functions, which can be used for fitting.

If you add functions here, do not forget to update the lookup table at the bottom as well.

ncempy.algo.math.const(x, param)[source]

Constant function.

f = param[0] .

Parameters:
  • x (np.ndarray) – Positions at which to evaluate function.
  • param (np.ndarray) – Necessary parameters.
Returns:

Values of function at x.

Return type:

(np.ndarray)

ncempy.algo.math.linear(x, param)[source]

Linear function.

f = param[0] * x + param[1] .

Parameters:
  • x (np.ndarray) – Positions at which to evaluate function.
  • param (np.ndarray) – Necessary parameters.
Returns:

Values of function at x.

Return type:

(np.ndarray)

ncempy.algo.math.lkp_funcs = {'const': (<function const>, 1), 'linear': (<function linear>, 2), 'powlaw': (<function powlaw>, 2), 'voigt': (<function voigt>, 4)}

Look-up table for functions implemented in this module. Functions are identified by strings, the entries give handles to the functions as well as the number of arguments necessary.

Type:(dict)
ncempy.algo.math.powlaw(x, param)[source]

Power law.

f = param[0] * x^param[1] .

Parameters:
  • x (np.ndarray) – Positions at which to evaluate function.
  • param (np.ndarray) – Necessary parameters.
Returns:

Values of function at x.

Return type:

(np.ndarray)

ncempy.algo.math.sum_functions(x, funcs, param)[source]

Sum functions in funcs over range x.

Parameters:
  • x (np.ndarray) – Positions at which to evaluate the functions.
  • funcs (list) – List of strings identifying function implemented in ncempy.algo.math.
  • param (np.ndarray) – Concatenated parameters for functions in funcs.
Returns:

Values of sum of functions at x.

Return type:

(np.ndarray)

ncempy.algo.math.voigt(x, param)[source]

Voigt peak function.

f = param[0]/(param[2] * sqrt(2*pi)) * Re( Faddeeva( ((x - param[1]) + i*param[3])/(param[2] * sqrt(2)) ) )

Parameters:
  • x (np.ndarray) – Positions at which to evaluate function.
  • param (np.ndarray) – Necessary parameters.
Returns:

Values of function at x.

Return type:

(np.ndarray)

ncempy.algo.moments module

Calculate the image raw moments.

param im:A 2D array of the image to calculate the moments for
type im:ndarray
param order:The order to calculate the moments to.
type order:int, default=3
returns:An ndarray of (order+1, order+1)
rtype:ndarray

ncempy.algo.multicorr_funcs module

Module to correlate two images, functionally written.

Todo

  • Cant use rfft2 currently. This gives one shift as 1/2 the value. How can this be improved to improve speed?
ncempy.algo.multicorr_funcs.dftUpsample(image_corr, upsample_factor, xy_shift)[source]

This performs a matrix multiply DFT around a small neighboring region of the initial correlation peak. By using the matrix multiply DFT to do the Fourier upsampling, the efficiency is greatly improved. This is adapted from the subfunction dftups found in the dftregistration function on the Matlab File Exchange.

https://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation

The matrix multiplication DFT is from Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, “Efficient subpixel image registration algorithms,” Opt. Lett. 33, 156-158 (2008). http://www.sciencedirect.com/science/article/pii/S0045790612000778

Parameters:
  • image_corr (ndarray) – Correlation image between two images in Fourier space.
  • upsample_factor (int) – Scalar integer of how much to upsample.
  • xy_shift (list of 2 floats) – Single pixel shift between images previously computed. Used to center the matrix multiplication on the correlation peak.
Returns:

image_upsample – Upsampled image from region around correlation peak.

Return type:

ndarray

ncempy.algo.multicorr_funcs.imageShifter(g1, xy_shift)[source]

Multiply im by a plane wave that has the real space effect of shifting ifft2(G2) by [x, y] pixels.

Parameters

g1 : complex ndarray
The Fourier transform of an image.
xy_shift : list
A two element list of the shifts along each axis.
G2shift : complex ndarray
Fourier shifted image FFT
>> shiftIm0 = np.real(np.fft.ifft2(multicorr.imageShifter(np.fft.fft2(im0),[11.1,22.2]))) >> plt.imshow(shiftIm0)
ncempy.algo.multicorr_funcs.initial_correlation_image(g1, g2, method='cross', verbose=False)[source]

Generate correlation image at initial resolution using the method specified.

Parameters:
  • g1 (complex ndarray) – Fourier transform of reference image.
  • g2 (complex ndarray) – Fourier transform of the image to register (the kernel).
  • method (str, optional) – The correlation method to use. Must be ‘phase’ or ‘cross’ or ‘hybrid’ (default = ‘cross’)
  • verbose (bool, default is False) – Print output.
Returns:

imageCorr – Correlation array which has not yet been inverse Fourier transformed.

Return type:

ndarray complex

ncempy.algo.multicorr_funcs.upsampleFFT(image_init, upsample_factor)[source]

This does a Fourier upsample of the imageInit. imageInit is the Fourier transform of the correlation image. The function returns the real space correlation image that has been Fourier upsampled by the upsample_factor. An upsample factor of 2 is generally sufficient.

The way it works is that it embeds imageInit in a larger array of zeros, then does the inverse Fourier transform to return the Fourier upsampled image in real space.

Parameters:
  • image_init (ndarray complex) – The image to be Fourier upsampled. This should be in the Fourier domain.
  • upsample_factor (int) – THe upsample factor (usually 2).
Returns:

imageUpsampleReal – The inverse Fourier transform of imageInit upsampled by the upsampleFactor.

Return type:

ndarray complex

imageSize = imageInit.shape imageUpsample = np.zeros(tuple((i*upsampleFactor for i in imageSize))) + 0j imageUpsample[:imageSize[0], :imageSize[1]] = imageInit imageUpsample = np.roll(np.roll(imageUpsample, -int(imageSize[0]/2), 0), -int(imageSize[1]/2),1) imageUpsampleReal = np.real(np.fft.ifft2(imageUpsample)) return imageUpsampleReal

ncempy.algo.multicorr_funcs.upsampled_correlation(image_corr, upsample_factor, verbose=False)[source]

Upsamples the correlation image by a set integer factor upsample_factor. If upsample_factor == 2, then it is naively Fourier upsampled. If the upsample_factor is higher than 2, then it uses dftUpsample, which is a more efficient way to Fourier upsample the image.

Parameters:
  • image_corr (ndarray complex) – Fourier transformed correlation image returned by initial_correlation_image.
  • upsample_factor (int) – Upsampling factor.
  • verbose (bool) – Provide output for debugging
Returns:

xyShift – Shift in x and y of G2 with respect to G1.

Return type:

list

ncempy.algo.radial_profile module

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.
ncempy.algo.radial_profile.calc_polarcoords(center, dims, ns=None, dists=None)[source]

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:

Polar coordinates (r, theta) as two np.ndarrays with same dimensions as original image.

Return type:

tuple

ncempy.algo.radial_profile.calc_radialprofile(img, rs, rMax, dr, rsigma, mask=None)[source]

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 of radial and intensity axes.

Return type:

tuple

ncempy.algo.radial_profile.correct_distortion(img, dims, center, ns, dists)[source]

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:

Corrected diffraction pattern.

Return type:

np.ndarray

ncempy.algo.radial_profile.fit_radialprofile(r, intens, funcs, init_guess, maxfev=None)[source]

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:

Optimized parameters.

Return type:

np.ndarray

ncempy.algo.radial_profile.residuals_fit(param, r, intens, funcs)[source]

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:

Residuals.

Return type:

np.ndarray

ncempy.algo.radial_profile.run_singleImage(img, dims, settings, show=False)[source]

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)

Return type:

Optimized parameters of fitting the radial profile according to the settings.

ncempy.algo.rebin module

Rebin a 2D array. From stackoverflow: https://stackoverflow.com/questions/4624112/grouping-2d-numpy-array-in-average

This should be used carefully. There is not antialiasing applied which could produce odd effects for some data sets (such as high resolution data)

param im:2D array to reduce
type im:ndarray
param f:The factor to rebin by. Must be integer
type f:int
param funcType:The type of reduction. mean or sum are implemented.
type funcType:str
returns:The 2D array with the new size.
rtype:ndarray

ncempy.algo.gaussND module

A module with Gaussian functions for creating Gaussian distributions in 1D, 2D and 3D. There are also functions with the same name and _FIT which can be used to fit an intensity distribution in a ndarray with the corresponding Gaussian distribution using scipy’s curve_fit.

A 1D Lorentz and a 1D gaussLorentz function is also included for fitting zero loss peaks in EELS spectra.

These functions were written assuming you use np.meshgrid() with indexing = ‘xy’. If you use ‘ij’ indexing then be careful about how you pass in the x and y coordinates.

ncempy.algo.gaussND.MatrixQuaternionRot(vector, theta)[source]

Quaternion tp rotate a given theta angle around the given vector.

adapted from code by Yongsoo Yang, yongsoo.ysyang@gmail.com

Parameters:
  • vector (ndarray, 3-element) – A non-zero 3-element numpy array representing rotation axis
  • theta (float, degrees) – Rotation angle in degrees
Returns:

Authoryongsoo.ysyang@gmail.com

Return type:

Yongsoo Yang, Dept. of Physics and Astronomy, UCLA

ncempy.algo.gaussND.gauss1D(x, x0, sigma)[source]

Returns the value of a gaussian at a 2D set of points for the given standard deviation with maximum normalized to 1.

Parameters:
  • x (ndarray) – A vector of size (N,) of points for the x values
  • x0 (float) – The center of the Gaussian.
  • sigma (float) – Standard deviation of the Gaussian.
Returns:

g – A vector of size (N,) of the Gaussian distribution evaluated at the input x values.

Return type:

ndarray

Note

Calculate the half width at half maximum (HWHM) as >> HWHM = sqrt(2*log(2))*stDev ~ 1.18*stDev or if x goes from -1 to 1 >> 0.5*size(x)*stDev

ncempy.algo.gaussND.gauss2D(x, y, x0, y0, sigma_x, sigma_y)[source]

Calculates the value of a Gaussian at a 2D set of points for the given standard deviations with maximum normalized to 1. The Gaussian axes are assumed to be 90 degrees from each other.

Parameters:
  • y (x,) – A 2D array of size (M,N) of points for the x values. Use x, y = np.meshgrid(range(M),range(N),indexing=’xy’).
  • x0 (float) – The center of the Gaussian along the x direction.
  • y0 (float) – The center of the Gaussian along the y direction.
  • sigma_x (float) – Standard deviation of the Gaussian along x.
  • sigma_y (float) – Standard deviation of the Gaussian along y.
Returns:

g – A ndarray of size (N, M) of the Gaussian distribution evaluated at the input x values.

Return type:

ndarray

Note

The Gaussian is normalized such that the peak == 1. To normalize the integral divide by 2 * np.pi * sigma_x * sigma_y

ncempy.algo.gaussND.gauss2D_FIT(xy, x0, y0, sigma_x, sigma_y)[source]

Version of gauss2D used for fitting (1 x_data input (xy) and flattened output). Returns the value of a gaussian at a 2D set of points for the given standard deviation with maximum normalized to 1. The Gaussian axes are assumed to be 90 degrees from each other.

Parameters:
  • xy (tuple) – A (N,M) shaped tuple containing the vectors of points for the evaluation points.
  • x0 (float) – The x center of the Gaussian.
  • y0 (float) – The y center of the Gaussian.
  • sigma_x (float) – The standard deviation of the Gaussian along x.
  • sigma_y (float) – The standard deviation of the Gaussian along y.
Returns:

g2_norm – The Gaussian distribution with maximum value normalized to 1. The 2D ndarray is reshaped into a (N*M,) array for use in fitting functions in numpy and scipy.

Return type:

ndarray

ncempy.algo.gaussND.gauss2D_poly_FIT(xy, x0, y0, A, B, C)[source]

Returns the flattened values of an elliptical gaussian at a 2D set of points for the given polynomial pre-factors with maximum normalized to 1. The Gaussian axes are assumed to be 90 degrees from each other.

The matrix looks like [[A,B],[B,C]]. See https://en.wikipedia.org/wiki/Gaussian_function

Parameters:
  • xy (tuple) – Evaluation points along x and y of shape (N,M)
  • x0 (float) – Center of the maximum along x.
  • y0 (float) – Center of the maximum along y.
  • A (float) – The A pre-factor for the polynomial expansion in the exponent.
  • B (float) – The B pre-factor for the polynomial expansion in the exponent.
  • C (float) – The C pre-factor for the polynomial expansion in the exponent.
Returns:

g2_norm – A (N, M) sized 2D ndarray with a Gaussaian distribution rotated.

Return type:

ndarray

ncempy.algo.gaussND.gauss2D_theta(x, y, x0, y0, sigma_x, sigma_y, theta)[source]

Returns the value of a generalized Gaussian at a 2D set of points for the given standard deviation with maximum normalized to 1. The Gaussian axes (assumed to be 90 degrees from each other) can be oriented at different angles to the output array axes using theta.

Parameters:
  • x (ndarray or list) – Evaluation points along x of size (N,)
  • y (ndarray) – Evaluation points along y of size (M,)
  • x0 (float) – Center of the maximum along x.
  • y0 (float) – Center of the maximum along y.
  • sigma_x (float) – Standard deviation along x.
  • sigma_y (float) – Standard deviation along y.
  • theta (float) – The rotation of the Gaussian principle axes from the array horizontal and vertical. In radians.
Returns:

g2_norm – A (N, M) sized 2D ndarray with a Gaussian distribution rotated.

Return type:

ndarray

ncempy.algo.gaussND.gauss2D_theta_FIT(xy, x0, y0, sigma_x, sigma_y, theta)[source]

Version of gauss2D_theta used for fitting (1 x_data input and flattened output). Returns the value of a gaussian at a 2D set of points for the given standard deviation with maximum normalized to 1. The Gaussian axes can be oriented at different angles (theta) in radians.

Parameters:
  • xy (tuple) – Evaluation points along x and y of shape (N,M)
  • x0 (float) – Center of the maximum along x.
  • y0 (float) – Center of the maximum along y.
  • sigma_x (float) – Standard deviation along x.
  • sigma_y (float) – Standard deviation along y.
  • theta (float) – The rotation of the Gaussian principle axes from the array horizontal and vertical. In radians.
Returns:

g2_norm – A (N, M) sized 2D ndarray with a Gaussian distribution rotated.

Return type:

ndarray

ncempy.algo.gaussND.gauss3D(x, y, z, x0, y0, z0, sigma_x, sigma_y, sigma_z)[source]

Returns the value of a gaussian at a 2D set of points for the given standard deviations with maximum normalized to 1. The Gaussian axes are assumed to be 90 degrees from each other.

Note

Be careful about the indexing used in meshgrid and the order in which you pass the x, y, z variables in.

Parameters:
  • y, z (x,) – 2D arrays of points (from meshgrid)
  • y0, z0 (x0,) – The x, y, z centers of the Gaussian
  • sigma_y, sigma_z (sigma_x,) – The standard deviations of the Gaussian.
Returns:

g3_norm – A 3D ndarray

Return type:

ndarray

ncempy.algo.gaussND.gauss3DGEN_FIT((x, y, z), x0, y0, z0, sigma_x, sigma_y, sigma_z, Angle1, Angle2, Angle3, BG, Height)[source]

Returns the value of a gaussian at a 3D set of points for the given sub-pixel positions with standard deviations, 3D Eular rotation angles, constant Background value, and Gaussian peak height.

adapted from code by Yongsoo Yang, yongsoo.ysyang@gmail.com

Note

This is a work in progress. Needs more testing.

Parameters:
  • xyz (tuple of 3 np.ndarray) – 3D arrays of points (from meshgrid) combined in a tuple
  • y0, z0 (x0,) – The x, y, z centers of the Gaussian
  • sigma_x,sigma_y,sigma_z (float) – standard deviations along x,y,z direction before 3D angular rotation
  • Angle2, Angle3 (Angle1,) – Tait-Bryan angles in ZYX convention for 3D rotation in degrees
  • BG (float) – Background
  • Height (float) – The peak height of the Gaussian function
Returns:

The 3D Gaussian

Return type:

ndarray, 3D

ncempy.algo.gaussND.gauss3D_FIT((x, y, z), x0, y0, z0, sigma_x, sigma_y, sigma_z)[source]

Returns the value of a gaussian at a 2D set of points for the given standard deviations with maximum normalized to 1. The Gaussian axes are assumed to be 90 degrees from each other.

xyz - x0, y0, z0 = the x, y, z centers of the Gaussian sigma_x, sigma_y, sigma_z = The std. deviations of the Gaussian.

Note

Be careful about the indexing used in meshgrid and the order in which you pass the x, y, z variables in.

Parameters:
  • xyz (tuple of ndarrays) – A tuple containing the 3D arrays of points (from meshgrid)
  • y0, z0 (x0,) – The x, y, z centers of the Gaussian
  • sigma_y, sigma_z (sigma_x,) – The standard deviations of the Gaussian.
Returns:

g3_norm – A flattened array for fitting.

Return type:

ndarray

ncempy.algo.gaussND.gauss3D_poly(x, y, z, x0, y0, z0, A, B, C, D, E, F)[source]

gauss3Dpoly_FIT((x,y,z),x0,y0,z0,A,B,C,D,E,F) Returns the value of a gaussian at a 2D set of points for the given standard deviations with maximum normalized to 1. The Gaussian axes are not locked to be 90 degrees.

Parameters:
  • y, z (x,) – 3D arrays of points (from meshgrid)
  • y0, z0 (x0,) – The x, y, z centers of the Gaussian
  • B, C, D, E, F (A,) – The polynomial values for the fit
Returns:

The 3D Gaussian

Return type:

ndarray, 3D

ncempy.algo.gaussND.gauss3D_poly_FIT(xyz, x0, y0, z0, A, B, C, D, E, F)[source]

gauss3Dpoly_FIT((x,y,z),x0,y0,z0,A,B,C,D,E,F) Returns the value of a gaussian at a 2D set of points for the given standard deviations with maximum normalized to 1. The Guassian axes are not locked to be 90 degrees.

Parameters:
  • xyz (tuple of ndarrays) – 3D arrays of points (from meshgrid) combined in a tuple
  • y0, z0 (x0,) – The x, y, z centers of the Gaussian
  • B, C, D, E, F (A,) – The polynomial values for the fit
Returns:

The 3D Gaussian

Return type:

ndarray, 3D

ncempy.algo.gaussND.gaussLorentz1D(x, x0, w)[source]

A Gaussian-Lorentzian function in one dimension.

Parameters:
  • x (ndarray or list) – A 1D vector of points for the dependent variable
  • x0 (float) – The center of the peak maximum
  • w (float) – The parameter that modifies the width of the distribution.
Returns:

lg – A vector of size (N,) of the Lorentzian Gaussian distribution evaluated at the input x values.

Return type:

ndarray

ncempy.algo.gaussND.lorentz1D(x, x0, w)[source]

Returns the probability density function of the Lorentzian (aka Cauchy) function with maximum at 1.

Parameters:
  • x (ndarray or list) – A 1D vector of points for the dependent variable
  • x0 (float) – The center of the peak maximum
  • w (float) – The parameter that modifies the width of the distribution.
Returns:

l – A vector of size (N,) of the Lorentzian distribution evaluated at the input x values.

Return type:

ndarray

ncempy.algo.peak_find module

Module to find local maxima in 2D and 3D data such as images and density maps.

Many functions are based on work by Colin Ophus, cophus@lbl.gov and his excellent RealSpaceLattice01.m code.

author: Peter Ercius, percius@lbl.gov

ncempy.algo.peak_find.applyLatticeLimit(lattice, bounds)[source]

Remove lattice points outside the data bounds. For 2D and 3D data.

Parameters:
  • lattice (ndarray; (N, 2) or (N, 3)) – From lattice2D
  • bounds (tuple,) – Minimum and maximum for axes 0 and 1 (min0, max0, min1, max1) or axes 0, 1 and 2 (min0, max0, min1, max1, min2, max2)
Returns:

Same as lattice input except only containing points within the bounds specified. M <= N

Return type:

ndarray; (M, 2) or (M, 3)

ncempy.algo.peak_find.calculate_unit_cell(image, lattice, u, v, unit_cell_size)[source]

IN PROGRESS Calculate a unit cell using the input lattice and lattice parameters. Any unit cell at the edge of the image is not used as part of the unit cell will be invalid.

Parameters:
  • image (ndarray) – The image containing the unit cell in a regular pattern (i.e. a lattice)
  • lattice (ndarray) – The starting position of each unit cell in the image.
  • v (u,) – The u and v vectors for the fitted lattice
  • unit_cell_size (int or tuple) – If int then this is the number of pixels in the unit cell along the u vector. The size of the v vector will be automatically calculated to make the unit cell pixels square. If tuple then the number of pixels along each u and v dimension are used directly.
Returns:

An array of size unit_cell_size which contains the average value of each position in the unit cell for each peak position and u,v, lattice point.

Return type:

ndarray

ncempy.algo.peak_find.doubleRoll(image, vec)[source]

Roll an 2D ndarray (e.g. an image) along each dimension.

Parameters:
  • image (np.ndarray) – The 2D array to roll.
  • vec (iterable) – The number of time to apply the roll to the rows and columns where the tuple has the structure (row, col). Values in tuple must be int.
Returns:

The rolled image.

Return type:

np.ndarray

ncempy.algo.peak_find.enforceMinDist(positions, intensities, minDistance)[source]

Remove peaks that violate a minimum distance requirement. If two peaks are closer together than the minDistance then the one with the highest intensity is chosen.

Works for 2D and 3D data sets.

Parameters:
  • positions (np.ndarray) – An array of shape Nx3 or Nx2 containing the positions of the peaks to test
  • intensities (np.ndarray) – An array of shape Nx1 with the peak intensities.
  • minDistance (int) – The minimum distance in pixels
Returns:

validPeaks – An ndarray of Mx3 or Mx2 (M <= N) with invalid peaks removed.

Return type:

np.ndarray

ncempy.algo.peak_find.fit_peaks_gauss2d(image, peaks, cutOut, init, bounds, remove_edge_peaks=True)[source]

Fit peaks to a 2D Gaussian. The Gaussian function is gaussND.gauss2D()

Todo: Write tests. This function is copied from test_peakFind.ipynb Todo: Add example test_peakFind.ipynb from my jupyter notebook folder

Parameters:
  • image (np.ndarray) – The 2D image with the intensities to fit to
  • peaks (np.ndarray) – The peaks in a ndarray of shape (N,2) where N is the number of peaks. Should match output of peakFind.peakFind2D.
  • cutOut (int) – The size (+/-) of the region around each peak to fit.
  • init (tuple) – A tuple of initial sigma values in x and y. (sigma_x, sigma_y)
  • bounds (tuple) – A tuple of 2 tuples indicating the upper and lower bounds for the fitting. See scipy.optimize.curve_fit for more information. Ordered as ( (center_x_low, center_y_low, sigma_x_low, sigma_y_low), (center_x_high, center_y_high, sigma_x_high, sigma_y_high))
  • remove_edge_peaks (bool, default = True) – Peak positions at the edge of the image are set to NaN and are removed. If you want to have those positions returned as nan set this to False.

Note

This function uses np.meshgrid and indexing=’ij’ internally.

Returns:
Returns a tuple of 3 arrays:
[0] optimized peak positions as a (M, 2) ndarray. M <= N. Peaks on the edge of the image are removed. [1] peak intensity value interpolated at the optimized peak position. [2] the fitting values for each peak
Return type:tuple
ncempy.algo.peak_find.fit_peaks_gauss3d(volume, peaks, cutOut, init, bounds, remove_edge_peaks=True)[source]

Fit peaks to a 2D Gaussian. The Gaussian function is gaussND.gauss2D()

Todo: Write tests. This function is copied from test_peakFind_3d.ipynb

Parameters:
  • volume (np.ndarray) – The 3D image with the intensities to fit to.
  • peaks (np.ndarray) – The peaks in a ndarray of shape (N, 3) where N is the number of peaks. Should match output of peakFind.peakFind3D.
  • cutOut (int) – The size (+/-) of the region around each peak to fit.
  • init (tuple) – A 3-tuple of initial sigma values in x and y. (sigma_x, sigma_y, sigma_z)
  • bounds (tuple) – A 2-tuple of 6-tuples indicating the upper and lower bounds for the fitting. See scipy.optimize.curve_fit for more information. Ordered as ( (center_x_low, center_y_low, center_z_low, sigma_x_low, sigma_y_low, sigma_z_low), (center_x_high, center_y_high, center_z_high, sigma_x_high, sigma_y_high, sigma_z_high))
  • remove_edge_peaks (bool, default = True) – Peak positions at the edge of the image are set to NaN and are removed. If you want to have those positions returned as nan set this to False.

Note

This function uses np.meshgrid and indexing=’ij’ internally.

Returns:
Returns a tuple of 3 arrays:
[0] optimized peak positions as a (M, 3) ndarray. M <= N. Peaks on the edge of the image are removed unless otherwise specified. Then the edge peaks are returned as NAN values. [1] peak intensity value interpolated at the optimized peak position. [2] the local fitting values for each peak.
Return type:tuple
ncempy.algo.peak_find.generateLatticeFromRefinement(origin, u, v, ab, fraction=(1, 1))[source]

Generate lattice positions from refined output of refineLattice2D. This can be used to generate all positions when fraction is used in refineLattice2D

Parameters:
  • origin (tuple) – Origin. 2-tuple
  • u (tuple) – Initial u vector. 2-tuple
  • v (tuple) – Initial v vector. 2-tuple
  • fraction (tuple) – The fraction used in refineLattice2D. 2-tuple
  • ab (np.ndarray) – The set of positions in terms of u and v lattice vectors. (number, position) (M, 2)
  • fraction – The fractional coordinates for the unit cell. 2-tuple
Returns:

The lattice site positions of the given lattice vectors of shape (M, 2)

Return type:

np.ndarray

ncempy.algo.peak_find.lattice2D(u, v, a, b, origin, num_points)[source]

A modified version of peakFind.lattice2D_norm which can use non normalized u,v,vectors

Parameters:
  • v (u,) – 2 element vectors defining the lattice directions
  • b (a,) – values to multiply each vector by (if u,v,w are not normalized then set these to 1)
  • origin (tuple) – The origin in the format (x0,y0)
  • num_points (tuple) – 2 element tuple of the number of repeats of the lattice along each direction
Returns:

The set of points in the lattice. Size (num_points[0] * num_points[1], 2)

Return type:

ndarray

ncempy.algo.peak_find.lattice2D_2(u, v, a, b, xy0, numPoints)[source]

A modified version of lattice2D which uses non-normalized u, v vectors

ncempy.algo.peak_find.lattice2D_norm(u, v, a, b, origin, num_points)[source]

Returns a set of points in a lattice according to the u, v unit vectors (vectors are normalized internally) and lengths a,b centered at origin. The lattice has num_points along each u,v vector.

Parameters:
  • v (u,) – 2 element tuples defining the lattice directions as vectors.
  • b (a,) – values to multiply each vector by (if u,v are not normalized then set these to 1)
  • origin (2-tuple) – The origin in the format (x0, y0)
  • num_points (2-tuple) – 2 element tuple of the number of repeats of the lattice along each direction
Returns:

The set of points in the lattice. Size (num_points[0] * num_points[1], 2)

Return type:

ndarray

ncempy.algo.peak_find.lattice3D(u, v, w, a, b, c, origin, num_points)[source]

Returns a set of points in a lattice according to the u, v, w vectors and lengths a,b centered at origin. The lattice has num_points along each u,v,w vector.

Parameters:
  • v, w (u,) – 3 element vectors defining the lattice directions
  • b, c (a,) – Values to multiply each vector by (if u,v,w are not normalized then set these to 1)
  • origin (tuple) – The origin in the format (x0,y0,z0)
  • num_points (tuple) – 3 element tuple of the number of repeats of the lattice along each direction
Returns:

xyz – A (N,3) shaped set of coordinates

Return type:

ndarray

ncempy.algo.peak_find.lattice3D_2(u, v, w, a, b, c, xyz0, numPoints)[source]

A modified version of lattice3D which uses non-normalized u,v vectors

ncempy.algo.peak_find.latticeDisplacements(peaks, u, v, origin)[source]

Find the displacements of the experimental peaks to the fitted lattice parameters

Parameters:
  • peaks (ndarray) – The experimental peaks with shape (num_peaks, 2). For example, the output of fit_peaks_gauss2d.
  • v (u,) – The u and v vectors for the fitted lattice
  • origin (tuple) – The origin of the fitted lattice for u, v vectors
Returns:

The displacement of each peak from the expected lattice position.

Return type:

ndarray

ncempy.algo.peak_find.match_lattice_peaks(peaks, u, v, origin)[source]

Find the matching lattice points to the experimental peaks. This is useful to generate all fitted lattice points for a set of experimental peaks. This can then be used directly to calculate displacements for example.

Parameters:
  • peaks (ndarray) – The experimental peaks as a (num_peaks, 2) ndarray
  • v (u,) – The u and v vectors for the fitted lattice
  • origin (tuple) – The origin of the fitted lattice for u, v vectors
Returns:

An array of shape (num_peaks, 2) where each coordinate is the closest lattice point for each input peak.

Return type:

ndarray

ncempy.algo.peak_find.peak3View(fg, vol, peak_positions)[source]

Plot 3 orthogonal slices and the corresponding peaks in each slice.

Not fully test. Unsure whether the slices and peaks are exactly the same.

todo: move this to viz

Parameters:
  • fg (matplotlib figure) – The figure to use. This allows the user to reuse figures
  • vol (ndarray (3D)) – The volume to slice along each orthogonal direction. Set up for [X,Z,Y] orientation
  • peak_positions (ndarray, list, tuple (1D)) – This selects the point shown in all three slices
Returns:

Return type:

None

ncempy.algo.peak_find.peakFind2D(image, threshold)[source]

Find peaks in a 2D image based on the roll technique. A threshold is first applied to remove low intensity noisy peaks.

Parameters:
  • image (ndarray) – An 2D ndarray of image intensities with peaks to find.
  • threshold (float) – Threshold [0,1) the image to remove noisy peaks based on max intensity in the image
Returns:

positions – A Nx2 array with the index locations of the n number of peaks

Return type:

array

ncempy.algo.peak_find.peakFind3D(vol, threshold)[source]

Find peaks in a 3D volume based on the roll technique. A threshold is first applied to remove low intensity noisy peaks.

Parameters:
  • vol (ndarray) – An 3D ndarray of image intensities with peaks to find.
  • threshold (float) – Threshold [0,1) the image to remove noisy peaks based on max intensity in the image
Returns:

positions – A Nx3 array with the index locations of the n number of peaks

Return type:

ndarray

ncempy.algo.peak_find.peakPlot3D(X, Y, Z, mkr, myAxes3D)[source]

Plot a set of peaks in a 3D plot using matplotlib. See the example below for how to set up a figure as input to this function using matplotlib Axes3D.

todo: move this to viz

Parameters:
  • X (np.ndarray) – The x positions as a 1D array.
  • Y (np.ndarray) – The y positions as a 1D array.
  • Z (np.ndarray) – The Z positions as a 1D array.
  • mkr (string) – String indicating the marks type and color. Passed directly to pyplot.plt command (ex. ‘go’ for green circles)
  • myAxes3D (mpl_toolkits.mplot3d.Axes3d) – An Axes3D object. (see example below)
Returns:

Return type:

None

Example

This function requires the input of a Axes3D to plot a set of peaks. See below how to import and set up a figure for use with this function. >> import matplotlib.pyplot as plt >> import mpl_toolkits.mplot3d >> from ncempy.algo import peak_find >> fg1 = plt.figure() >> ax1 = mpl_toolkits.mplot3d.Axes3D(fg1) >> peak_find.peakPlot3D(peakList[:,2], peakList[:,1], peakList[:,0], ‘go’, ax1)

ncempy.algo.peak_find.peaksToImage(peakList, imShape, gaussSigma, gaussSize, indexing='ij')[source]

Place 2D Gaussian at a set of peak positions.

Parameters:
  • peakList (np.ndarray) – Array of peak positions of size (numPeaks, 2).
  • imShape (tuple) – 2 element tuple of the shape of the image containing the peaks
  • gaussSigma (tuple) – 2 element tuple with that sigma parameters passed to the gaussND.gauss2D() function (sigX,sigY)
  • gaussSize (tuple) – 2 element tuple describing the 2D size of the simulated gaussian box. Must be odd.
  • indexing (str, default = 'ij') – ‘ij’ or ‘xy’ indexing to pass to np.meshgrid.
Returns:

A image containing Gaussian peaks at each position indicated in peakList parameter.

Return type:

np.ndarray

ncempy.algo.peak_find.peaksToVolume(peakList, volShape, gaussSigma, gaussSize, indexing='ij')[source]

Place 3D Gaussian at a set of peak positions.

Parameters:
  • peakList (np.ndarray) – Array of peak positions of size (numPeaks, 3).
  • volShape (tuple) – 3 element tuple of the shape of the volume containing the peaks
  • gaussSigma (tuple) – 3 element tuple with that sigma parameters passed to the gaussND.gauss3D() function (sigX,sigY,sigZ)
  • gaussSize (tuple) – 3 element tuple describing the 3D size of the simulated gaussian box. Must be odd.
  • indexing (str) – String to pass to np.meshgrid to indicate which indexing is used. Either ‘ij’ or ‘xy’ are possible. Default is ‘ij’ for this module.
Returns:

A volume containing gaussian peaks at each position indicated in peakList parameter.

Return type:

np.ndarray

ncempy.algo.peak_find.refineLattice2D(or0, u0, v0, pos, fraction=(1, 1), max_iter=30, refine_locally=True, verbose=False, num_unit_cells=(1, 1))[source]

Refine lattice based on measurements and initial guess at lattice vectors. This code is designed to work only with square lattices. Hexagonal and more complex lattice might not work.

Parameters:
  • or0 (tuple) – Initial guess of origin. 2-tuple
  • u0 (tuple) – Initial u vector. Use num_unit_cells to average over several lengths of the lattice vector in this direction. 2-tuple
  • v0 (tuple) – Initial v vector. Use num_unit_cells to average over several lengths of the lattice vector in this direction. 2-tuple
  • pos (np.ndarray) – The set of positions to fit the lattice to. (number, position) (M, 2)
  • refine_locally (bool) – Refine locally near the origin before using all positions Locally is considered 2 time the larger vector of u0 or v0.
  • fraction (tuple) – Site fraction to take into account peak positions inside the unit cell. FCC imaged along [100] for example would need site_fraction = (2, 2).
  • max_iter (int) – The maximum number of iterations to run to refine. This usually converges in a few iterations. Default is 30.
  • num_unit_cells (tuple) – The number of unit cells the initial guess of u0 and v0 extend over. 2-tuple (num_u0, num_v0) with default
  • verbose (bool) – Print out helpful information.
Returns:

Tuple of of 4 np.ndarray (origin, u, v, ab) where ab is the site positions in fractions of u and v.

Return type:

tuple

ncempy.algo.peak_find.refineLattice3D(or0, u0, v0, w0, pos, fraction=(1, 1, 1), max_iter=30, refine_locally=True)[source]

Refine lattice based on measurements and initial guess at lattice vectors

Warning

Tested code but not production yet. This is a work un progress.

Parameters:
  • or0 (3-element tuple) – Origin.
  • u0 (tuple) – Initial u vector.
  • v0 (tuple) – Initial v vector
  • w0 (tuple) – Initial w vector
  • pos (array (number, position) (M,3)) – The set of positions to fit the lattice to.
  • refine_locally (bool) – Refine locally near the origin before using all positions Locally is considered 2 time the larger vector of u0 or v0.
  • fraction (tuple) – Site fraction to take into account peak positions inside the unit cell. FCC imaged along [100] for example would need site_fraction = (2, 2).
  • max_iter (int) – Maximum number of iterations.
Returns:

Tuple of refined origin, u, v optimized values as np.ndarray and the site positions in fractions of u and v.

Return type:

tuple of 4 np.ndarray (origin, u, v, ab)

ncempy.algo.peak_find.remove_xrays(imageOriginal, threshold, size_median_filter=(3, 3))[source]

Find x-ray pixels in TEM images and replace the xray intensity with a local median value.

Parameters:
  • imageOriginal (ndarray) – The image as a 2D ndarray.
  • threshold (float) – Value above which to consider a pixel intensity an x-ray.
  • size_median_filter (list or tuple) – The footprint of the median filter over which to search for a median value.
Returns:

The image with x-ray values replaced with median local value.

Return type:

ndarray

ncempy.algo.peak_find.tripleRoll(vol, vec)[source]

Roll a 3D ndarray (e.g. an volume) along each dimension.

Parameters:
  • vol (np.ndarray) – The 3D array to roll
  • vec (iterable) – The number of time to apply the roll the tuple has the structure (axis0, axis1, axis2). Values in tuple must be int.
Returns:

The rolled volume.

Return type:

ndarray

ncempy.algo.peak_find.writeXYZ(filename, XYZ, element, comment)[source]

Write out a set of XYZ coordinates that can be read by various crystal viewing software such as Vesta.

todo: Move this to io

Parameters:
  • filename (str) – The name of the file to write to.
  • XYZ (np.ndarray) – Atom coordinates in a 2D ndarray with axes [atom number, position]. Should be shape [N, 3]
  • element (tuple) – Element symbol of each coordinate as a string (e.g. ‘C’). Should be length N.
  • comment (str) – A comment to add to the file.

Example

# Write out the peak positions as Carbon atoms: >> writeXYZ(‘name.xyz’, positions, [‘C’,] * len(positions), ‘writeXYZ example carbon atoms’)