# Source code for ncempy.algo.gaussND

""" 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.

"""

import numpy as np

[docs]def gauss1D(x, x0, sigma):
""" 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: ndarray
A vector of size (N,) of the Gaussian distribution evaluated
at the input x values.

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

"""
return np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

[docs]def lorentz1D(x, x0, w):
""" 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: ndarray
A vector of size (N,) of the Lorentzian distribution evaluated
at the input x values.

"""
return w ** 2 / ((x - x0) ** 2 + w ** 2)

[docs]def gaussLorentz1D(x, x0, w):
""" 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: ndarray
A vector of size (N,) of the Lorentzian Gaussian distribution
evaluated at the input x values.

"""
gSigma = w / (2 * np.sqrt(2 * np.log(2)))  # change the width to sigma
return gauss1D(x, x0, gSigma) + lorentz1D(x, x0, w)

[docs]def gauss2D(x, y, x0, y0, sigma_x, sigma_y):
""" 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
----------
x, y : ndarray
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 : ndarray
A ndarray of size (N, M) of the Gaussian distribution evaluated
at the input x values.

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

"""
x0 = float(x0)
y0 = float(y0)
g2 = np.exp(-((x - x0) ** 2 / (2 * sigma_x ** 2) + (y - y0) ** 2 / (2 * sigma_y ** 2)))
g2_norm = g2 / np.max(g2.flatten())
return g2_norm

[docs]def gauss2D_FIT(xy, x0, y0, sigma_x, sigma_y):
""" 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: ndarray
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.

"""
x0 = float(x0)
y0 = float(y0)
x = xy
y = xy
g2 = np.exp(-((x - x0) ** 2 / (2 * sigma_x ** 2) + (y - y0) ** 2 / (2 * sigma_y ** 2)))
g2_norm = g2 / np.max(g2.flatten())
return g2_norm.reshape(-1)

[docs]def gauss2D_theta(x, y, x0, y0, sigma_x, sigma_y, theta):
""" 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: ndarray
A (N, M) sized 2D ndarray with a Gaussian distribution rotated.

"""

x0 = float(x0)
y0 = float(y0)
a = (np.cos(theta) ** 2) / (2 * sigma_x ** 2) + (np.sin(theta) ** 2) / (2 * sigma_y ** 2)
b = -(np.sin(2 * theta)) / (4 * sigma_x ** 2) + (np.sin(2 * theta)) / (4 * sigma_y ** 2)
c = (np.sin(theta) ** 2) / (2 * sigma_x ** 2) + (np.cos(theta) ** 2) / (2 * sigma_y ** 2)
g2 = np.exp(- (a * ((x - x0) ** 2) + 2 * b * (x - x0) * (y - y0) + c * ((y - y0) ** 2)))
g2_norm = g2 / np.max(g2.flatten())
return g2_norm  # return a 2D array

[docs]def gauss2D_theta_FIT(xy, x0, y0, sigma_x, sigma_y, theta):
""" 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: ndarray
A (N, M) sized 2D ndarray with a Gaussian distribution rotated.

"""
x0 = float(x0)
y0 = float(y0)
x = xy
y = xy
a = (np.cos(theta) ** 2) / (2 * sigma_x ** 2) + (np.sin(theta) ** 2) / (2 * sigma_y ** 2)
b = -(np.sin(2 * theta)) / (4 * sigma_x ** 2) + (np.sin(2 * theta)) / (4 * sigma_y ** 2)
c = (np.sin(theta) ** 2) / (2 * sigma_x ** 2) + (np.cos(theta) ** 2) / (2 * sigma_y ** 2)
g2 = np.exp(- (a * ((x - x0) ** 2) + 2 * b * (x - x0) * (y - y0) + c * ((y - y0) ** 2)))
g2_norm = g2 / np.max(g2.flatten())
return g2_norm.ravel()  # return a 1D vector

[docs]def gauss2D_poly_FIT(xy, x0, y0, A, B, C):
""" 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: ndarray
A (N, M) sized 2D ndarray with a Gaussaian distribution rotated.
"""

x0 = float(x0)
y0 = float(y0)
x = xy  # retrieve the array from the tuple
y = xy
g2 = np.exp(-(A * (x - x0) ** 2 + 2 * B * (x - x0) * (y - y0) + C * (y - y0) ** 2))
g2_norm = g2 / np.max(g2.flatten())
return g2_norm.ravel()

[docs]def gauss3D(x, y, z, x0, y0, z0, sigma_x, sigma_y, sigma_z):
"""
gauss3D(x,y,z,x0,y0,z0,sigma_x,sigma_y,sigma_z)
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
----------
x, y, z: ndarray, from numpy.meshgrid
2D arrays of points (from meshgrid)

x0, y0, z0: float
The x, y, z centers of the Gaussian

sigma_x, sigma_y, sigma_z: float
The standard deviations of the Gaussian.
Returns
-------
g3_norm: ndarray
A 3D ndarray

"""
x0 = float(x0)
y0 = float(y0)
z0 = float(z0)
g3 = np.exp(-((x - x0) ** 2 / (2 * sigma_x ** 2) + (y - y0) ** 2 / (2 * sigma_y ** 2) + (z - z0) ** 2 / (
2 * sigma_z ** 2)))
g3_norm = g3 / np.max(g3.flatten())
return g3_norm

[docs]def gauss3D_FIT(xyz, x0, y0, z0, sigma_x, sigma_y, sigma_z):
"""
gauss3D_FIT((x,y,z),x0,y0,z0,sigma_x,sigma_y,sigma_z)
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)
x0, y0, z0: float
The x, y, z centers of the Gaussian
sigma_x, sigma_y, sigma_z: float
The standard deviations of the Gaussian.
Returns
-------
g3_norm: ndarray
A flattened array for fitting.

"""

x0 = float(x0)
y0 = float(y0)
z0 = float(z0)
x = xyz
y = xyz
z = xyz
g3 = np.exp(-((x - x0) ** 2 / (2 * sigma_x ** 2) + (y - y0) ** 2 / (2 * sigma_y ** 2) + (z - z0) ** 2 / (
2 * sigma_z ** 2)))
g3_norm = g3 / np.max(g3.flatten())
return g3_norm.ravel()

[docs]def gauss3D_poly(x, y, z, x0, y0, z0, A, B, C, D, E, F):
"""
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
----------
x, y, z : ndarray, 3D
3D arrays of points (from meshgrid)

x0, y0, z0 : float
The x, y, z centers of the Gaussian

A, B, C, D, E, F : float
The polynomial values for the fit

Returns
-------
: ndarray, 3D
The 3D Gaussian
"""
x0 = float(x0)
y0 = float(y0)
z0 = float(z0)
g3 = np.exp(-(A * (x - x0) ** 2 + B * (y - y0) ** 2 + C * (z - z0) ** 2 + 2 * D * (x - x0) * (y - y0) + 2 * E * (
x - x0) * (z - z0) + 2 * F * (y - y0) * (z - z0)))
g3_norm = g3 / np.max(g3.flatten())
return g3_norm

[docs]def gauss3D_poly_FIT(xyz, x0, y0, z0, A, B, C, D, E, F):
"""
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

x0, y0, z0 : float
The x, y, z centers of the Gaussian

A, B, C, D, E, F : float
The polynomial values for the fit

Returns
-------
: ndarray, 3D
The 3D Gaussian
"""
x0 = float(x0)
y0 = float(y0)
z0 = float(z0)
x = xyz
y = xyz
z = xyz
g3 = np.exp(-(A * (x - x0) ** 2 + B * (y - y0) ** 2 + C * (z - z0) ** 2 + 2 * D * (x - x0) * (y - y0) + 2 * E * (
x - x0) * (z - z0) + 2 * F * (y - y0) * (z - z0)))
g3_norm = g3 / np.max(g3.flatten())
return g3_norm.ravel()

[docs]def gauss3DGEN_FIT(xyz, x0, y0, z0, sigma_x, sigma_y, sigma_z, Angle1, Angle2, Angle3, BG, Height):
"""
gauss3DGEN_FIT((x,y,z),x0,y0,z0,sigma_x,sigma_y,sigma_z,Angle1,Angle2,Angle3,BG,Height)
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

x0, y0, z0 : float
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
Angle1, Angle2, Angle3 : float, degrees
Tait-Bryan angles in ZYX convention for 3D rotation in degrees
BG : float
Background
Height: float
The peak height of the Gaussian function

Returns
-------
: ndarray, 3D
The 3D Gaussian
"""
# 3D vectors for each sampled positions
print('Warning: this function is using Fortran ordered array for use in tomviz. Need to test this')
# Todo: Remove Fortran ordering
v = np.array([xyz.reshape(-1, order='F') - x0,
xyz.reshape(-1, order='F') - y0,
xyz.reshape(-1, order='F') - z0])

# rotation axes for Tait-Bryan angles
vector1 = np.array([0, 0, 1])
rotmat1 = MatrixQuaternionRot(vector1, Angle1)

vector2 = np.array([0, 1, 0])
rotmat2 = MatrixQuaternionRot(vector2, Angle2)

vector3 = np.array([1, 0, 0])
rotmat3 = MatrixQuaternionRot(vector3, Angle3)

# full rotation matrix
# Todo : Remove dependency on np.matrix()
rotMAT = np.matrix(rotmat3) * np.matrix(rotmat2) * np.matrix(rotmat1)

# 3x3 matrix for applying sigmas
D = np.matrix(np.array([[1. / (2 * sigma_x ** 2), 0, 0, ],
[0, 1. / (2 * sigma_y ** 2), 0],
[0, 0, 1. / (2 * sigma_z ** 2)]]))

# apply 3D rotation to the sigma matrix
WidthMat = np.transpose(rotMAT) * D * rotMAT

# calculate 3D Gaussian
RHS_calc = WidthMat * np.matrix(v)
Result = Height * np.exp(-1 * np.sum(v * RHS_calc.A, axis=0)) + BG

return Result

[docs]def MatrixQuaternionRot(vector, theta):
"""
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
-------

Author: Yongsoo Yang, Dept. of Physics and Astronomy, UCLA
yongsoo.ysyang@gmail.com
"""
theta = theta * np.pi / 180
vector = vector / np.float(np.sqrt(np.dot(vector, vector)))
w = np.cos(theta / 2)
x = -np.sin(theta / 2) * vector
y = -np.sin(theta / 2) * vector
z = -np.sin(theta / 2) * vector
RotM = np.array([[1. - 2 * y ** 2. - 2 * z ** 2, 2. * x * y + 2 * w * z, 2. * x * z - 2. * w * y],
[2. * x * y - 2. * w * z, 1. - 2. * x ** 2 - 2. * z ** 2, 2. * y * z + 2. * w * x],
[2 * x * z + 2 * w * y, 2 * y * z - 2. * w * x, 1 - 2. * x ** 2 - 2. * y ** 2]])

return RotM