cristallina_analysis_package/examples/2d_gaussian_fitting_example.ipynb

In [1]:

%matplotlib widget

from matplotlib import pyplot as plt
import numpy as np

import cristallina as cr

"Automatic p-group extraction from the current working directory didn't work."

Example 2D Gaussian fitting

To fit arbitrary Gaussian-shaped peaks on 2D detector images we use the lmfit-package with some additions.

Here is a short example:

In [2]:

def gauss2d(x=0, y=0, mx=0, my=0, sx=1, sy=1):
    """ Defines normalized 2D gaussian function.
    """ 
    return (
        1 / (2 * np.pi * sx * sy) * np.exp(-((x - mx) ** 2 / (2 * sx**2.0) + (y - my) ** 2 / (2 * sy**2)))
    )

x = np.arange(0, 150, 1)
y = np.arange(0, 100, 1)
x, y = np.meshgrid(x, y)  

z = gauss2d(x, y, mx=40, my=50, sx=20, sy=40)

In [3]:

fig, ax = plt.subplots()
ax.imshow(z, origin="lower")

Out [3]:

<matplotlib.image.AxesImage at 0x7fbe3f915d90>

In [4]:

center_x, center_y, result = cr.analysis.fit_2d_gaussian(z, plot=True)

In [5]:

result

Out [5]:

Fit Result

Model: Model(gaussian2d)

Fit Statistics

fitting method	leastsq
# function evals	52
# data points	15000
# variables	5
chi-square	9.3716e-36
reduced chi-square	6.2498e-40
Akaike info crit.	-1354057.78
Bayesian info crit.	-1354019.70
R-squared	1.00000000

Parameters

name	value	standard error	relative error	initial value	min	max	vary	expression
amplitude	1.00000000	5.6014e-18	(0.00%)	0.08151715916717606	-inf	inf	True
centerx	40.0000000	7.4667e-17	(0.00%)	40.0	-inf	inf	True
centery	50.0000000	1.7928e-16	(0.00%)	50.0	-inf	inf	True
sigmax	20.0000000	7.7413e-17	(0.00%)	24.833333333333332	0.00000000	inf	True
sigmay	40.0000000	2.7538e-16	(0.00%)	16.5	0.00000000	inf	True
fwhmx	47.0964000	1.8229e-16	(0.00%)	58.478030000000004	-inf	inf	False	2.3548200*sigmax
fwhmy	94.1927986	6.4847e-16	(0.00%)	38.854530000000004	-inf	inf	False	2.3548200*sigmay
height	1.9894e-04	7.9319e-22	(0.00%)	3.166285616970345e-05	-inf	inf	False	0.1591549*amplitude/(max(1e-15, sigmax)*max(1e-15, sigmay))

Correlations (unreported values are < 0.100)

Parameter1	Parameter 2	Correlation
amplitude	sigmay	+0.8115
amplitude	sigmax	+0.3522

Example 2D Gaussian with rotation

In [6]:

import lmfit
from scipy.interpolate import griddata

def gaussian2d_rot(
    x, y=0.0, amplitude=1.0, center_x=0.0, center_y=0.0, sigma_x=1.0, sigma_y=1.0, rotation=0, background=0
):
    """Returns a two-dimensional Gaussian model from lmfit with a rotation in radians around the center."""
    sr = np.sin(rotation)
    cr = np.cos(rotation)

    center_x_rot = center_x * cr - center_y * sr
    center_y_rot = center_x * sr + center_y * cr

    x_rot = x * cr - y * sr
    y_rot = x * sr + y * cr

    return (
        lmfit.models.gaussian2d(
            x_rot,
            y=y_rot,
            amplitude=amplitude,
            centerx=center_x_rot,
            centery=center_y_rot,
            sigmax=sigma_x,
            sigmay=sigma_y,
        )
        + background
    )


npoints = 5000
np.random.seed(2021)

x = np.random.rand(npoints) * 100
y = np.random.rand(npoints) * 50


height = 30  # /(2*np.pi*0.6*0.8)
z = gaussian2d_rot(x, y, height, 40, 30, 6, 20, 1.5)
#z += 0.2 * (np.random.rand(*z.shape) - 0.5)


# define normalized 2D gaussian
def gauss2d_rotated(x=0, y=0, center_x=0, center_y=0, sx=1, sy=1, rotation=0.5):

    sr = np.sin(rotation)
    cr = np.cos(rotation)

    center_x_rot = center_x * cr - center_y * sr
    center_y_rot = center_x * sr + center_y * cr

    x_rot = x * cr - y * sr
    y_rot = x * sr + y * cr

    return (1 / (2 * np.pi * sx * sy) * np.exp(-((x_rot - center_x_rot) ** 2 / (2 * sx**2.0) + (y_rot - center_y_rot) ** 2 / (2 * sy**2))))

x = np.arange(0, 150, 1)
y = np.arange(0, 100, 1)
x, y = np.meshgrid(x, y)  

z = 100*gauss2d_rotated(x, y, center_x=40, center_y=50, sx=10, sy=20, rotation=0.5)
z += 1E-2 * (np.random.rand(*z.shape) - 0.5)
        
# for general x and y that have a floating     
#X, Y = np.meshgrid(np.linspace(0, x.max(), 100), np.linspace(0, y.max(), 100))
#X, Y = np.meshgrid(np.arange(0, 100), np.arange(0, 50))
#Z = griddata((x, y), z, (X, Y), method="linear", fill_value=0)

fig, ax = plt.subplots()
art = ax.pcolor(x, y, z, shading="auto")

In [7]:

fig, ax = plt.subplots()
ax.imshow(z, origin='lower')

Out [7]:

<matplotlib.image.AxesImage at 0x7fbe3ebe4310>

In [8]:

center_x, center_y, result = cr.analysis.fit_2d_gaussian_rotated(z, vary_rotation=True, plot=True)

In [9]:

#fig, ax = plt.subplots()

#ax.plot(Z[30,:])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 3
      1 fig, ax = plt.subplots()
----> 3 ax.plot(Z[30,:])

NameError: name 'Z' is not defined