integrate elliptical gaussian.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import circulant
from scipy.optimize import curve_fit, brentq
from scipy.interpolate import interp1d
from scipy.ndimage import gaussian_filter
import scipy.ndimage as nd

def gaussian(x, a, mu, sigma, c):
    """
    Gaussian function
    
    a * exp(-((x - mu) ** 2) / 2 / sigma ** 2) + c
    """
    return a * np.exp(-((x - mu) ** 2) / 2 / sigma ** 2) + c


def guss_gaussian(x):
    """
    Find a set of better starting parameters for Gaussian function fitting
    
    """
    c_guess = (x[0] + x[-1]) / 2
    a_guess = x.max() - c_guess
    mu_guess = x.argmax()
    x_inter = interp1d(range(len(x)), x)

    def _(i):
        return x_inter(i) - a_guess / 2 - c_guess

    try:
        sigma_l_guess = brentq(_, 0, mu_guess)
    except:
        sigma_l_guess = len(x) / 4
    try:
        sigma_r_guess = brentq(_, mu_guess, len(x) - 1)
    except:
        sigma_r_guess = 3 * len(x) / 4
    return a_guess, mu_guess, (sigma_r_guess - sigma_l_guess) / 2.35482, c_guess


def fit_gaussian(x):
    """
    Fit Gaussian function and return its parameter 
    """
    p, q = curve_fit(gaussian, list(range(x.size)), x, p0=guss_gaussian(x))
    return p

def gaussian2D(X, Y, A=1, x0=0, y0=0, sx=1, sy=2):
    return A * np.exp(-(X-x0)**2/(2*sx**2) - (Y-y0)**2/(2*sy**2))
    

x = np.linspace(-5,5,200)
X,Y = np.meshgrid(x,x)

Z = gaussian2D(X,Y)


fig, axs = plt.subplots(2,1,figsize=(6,8))
axs[0].imshow(Z)

integration = np.sum(Z,axis=0)
axs[1].plot(integration, label='integral over y', lw=2)

p = fit_gaussian(integration)
fit = gaussian(np.arange(len(integration)),*p)
axs[1].plot(fit, label='Gaussian fit',color='r', ls='dashed', lw=2)

leg = axs[1].legend(); leg.draw_frame(False)

plt.show()

output: