"""Generic functions to generate and fit gaussians
based on 'CodedMaskImage.py' python program of Cyril Lachaud
modfied from http://scipy-cookbook.readthedocs.io/items/FittingData.html
"""
import logging
import numpy as np
from scipy import optimize
logger = logging.getLogger(__name__)
[docs]def gaussian_1d(height, center, sigma):
"""Returns a gaussian function + offset with the given parameters
:param height: height of the gaussian function
:type height: float
:param center: center of the gaussian function
:type center: float
:param sigma: sigma of the gaussian function
:type sigma: float
"""
return lambda x: height * np.exp(-((center - x) / sigma) ** 2 / 2.)
[docs]def gaussian_2d_fixedwidth_offset(height, center_x, center_y, offset):
"""Returns a gaussian function + offset with the given parameters
sigma is unique and fixed.
FWHM_fit = FWHM_PSF = sqrt(m**2+d**2) (=sqrt(m**2+1) en pixel)
m = mask element
d = detector pitch
sigma = sqrt(m**2+1)/2.3548 (with m=2.53)
sigma = 1.155 pixel
:param height: height of the gaussian function
:type height: float
:param center_x: center in x of the gaussian function
:type center_x: float
:param center_y: center in y of the gaussian function
:type center_y: float
:param offset: offset of the gaussian function
:type offset float
:return: gaussian function with the given parameters
:rtype: lambda function
"""
return lambda x, y, sigma: height * np.exp(-(((center_x - x) / sigma) ** 2 +
((center_y - y) / sigma) ** 2) / 2.) + offset
[docs]def gaussian_2d_1width_offset(height, center_x, center_y, sigma, offset):
"""Returns a gaussian function + offset with the given parameters
sigma is the same for x and y direction.
:param height: height of the gaussian function
:type height: float
:param center_x: center in x of the gaussian function
:type center_x: float
:param center_y: center in y of the gaussian function
:type center_y: float
:param sigma: sigma of the gaussian function
:type sigma: float
:param offset: offset of the gaussian function
:type offset float
:return: gaussian function with the given parameters
:rtype: lambda function
"""
sigma = float(sigma)
return lambda x, y, _: height * np.exp(-(((center_x - x) / sigma) ** 2 +
((center_y - y) / sigma) ** 2) / 2.) + offset
[docs]def gaussian_2d_2width_offset(height, center_x, center_y, sigma_x, sigma_y, offset):
"""Returns a gaussian function + offset with the given parameters
:param height: height of the gaussian function
:type height: float
:param center_x: center in x of the gaussian function
:type center_x: float
:param center_y: center in y of the gaussian function
:type center_y: float
:param sigma_x: sigma in x of the gaussian function
:type sigma_x: float
:param sigma_x: sigma in y of the gaussian function
:type sigma_x: float
:param offset: offset of the gaussian function
:type offset float
:return: gaussian function with the given parameters
:rtype: lambda function
"""
sigma_x = float(sigma_x)
sigma_y = float(sigma_y)
return lambda x, y, _: height * np.exp(-(((center_x - x) / sigma_x) ** 2 + \
((center_y - y) / sigma_y) ** 2) / 2.) + offset
[docs]def gaussian_2d_fixedwidth(height, center_x, center_y):
"""Returns a gaussian function with the given parameters
sigma is unique and fixed.
FWHM_fit = FWHM_PSF = sqrt(m**2+d**2) (=sqrt(m**2+1) en pixel)
m = mask element
d = detector pitch
sigma = sqrt(m**2+1)/2.3548 (with m=2.53)
sigma = 1.155 pixel
:param height: height of the gaussian function
:type height: float
:param center_x: center in x of the gaussian function
:type center_x: float
:param center_y: center in y of the gaussian function
:type center_y: float
:return: gaussian function with the given parameters
:rtype: lambda function
"""
return lambda x, y, sigma: height * np.exp(-(((center_x - x) / sigma) ** 2 + \
((center_y - y) / sigma) ** 2) / 2.)
[docs]def gaussian_2d_1width(height, center_x, center_y, sigma):
"""Returns a gaussian function with the given parameters
sigma is the same for x and y direction.
:param height: height of the gaussian function
:type height: float
:param center_x: center in x of the gaussian function
:type center_x: float
:param center_y: center in y of the gaussian function
:type center_y: float
:param sigma: sigma of the gaussian function
:type sigma: float
:return: gaussian function with the given parameters
:rtype: lambda function
"""
sigma = float(sigma)
return lambda x, y, _: height * np.exp(-(((center_x - x) / sigma) ** 2 + \
((center_y - y) / sigma) ** 2) / 2.)
[docs]def gaussian_2d_2width(height, center_x, center_y, sigma_x, sigma_y):
"""Returns a gaussian function with the given parameters
:param height: height of the gaussian function
:type height: float
:param center_x: center in x of the gaussian function
:type center_x: float
:param center_y: center in y of the gaussian function
:type center_y: float
:param sigma_x: sigma in x of the gaussian function
:type sigma_x: float
:param sigma_x: sigma in y of the gaussian function
:type sigma_x: float
:return: gaussian function with the given parameters
:rtype: lambda function
"""
sigma_x = float(sigma_x)
sigma_y = float(sigma_y)
return lambda x, y, _: height * np.exp(-(((center_x - x) / sigma_x) ** 2 + \
((center_y - y) / sigma_y) ** 2) / 2.)
[docs]def gaussian_2d_moments(data, fit_method, fit_sigma):
"""Guess the gaussian parameters of a distribution
Use for initial fitting parameters
the gaussian parameters of a 2D distribution by calculating its moments
Modified from : http://scipy-cookbook.readthedocs.io/items/FittingData.html
:param data: 2D distribution to evaluate
:type data: float 2D-array
:param fit_method: fit method that will be used: 'gaussian_2d_fixedwidth',
'gaussian_2d_1width' or 'gaussian_2d_2width'
:type fit_method: string
:param fit_sigma: standard deviation of the gaussian fit_function. Default is 1.555.
:type fit_sigma: float>0
:return: (height, x, y), (height, x, y, sigma_x, sigma_y) or (height, x, y, sigma) depending
on fit_method
:rtype: floats (number of which depend on the fit_method)
"""
height = np.nanmax(data)
x, y = np.where(data == height)
x = x[0]
y = y[0]
offset = 0 # TODO: should not we do the same for offset we did for sigma ???
if fit_method.__name__ == 'gaussian_2d_fixedwidth':
fit_params = height, x, y
elif fit_method.__name__ == 'gaussian_2d_1width':
fit_params = height, x, y, fit_sigma
elif fit_method.__name__ == 'gaussian_2d_2width':
fit_params = height, x, y, fit_sigma, fit_sigma
elif fit_method.__name__ == 'gaussian_2d_fixedwidth_offset':
fit_params = height, x, y, offset
elif fit_method.__name__ == 'gaussian_2d_1width_offset':
fit_params = height, x, y, fit_sigma, offset
elif fit_method.__name__ == 'gaussian_2d_2width_offset':
fit_params = height, x, y, fit_sigma, fit_sigma, offset
return fit_params
[docs]def fit_gaussian_2d(data, fit_method, fit_sigma, verbose=False):
"""fit the 2D distribution with 2D gaussian function
The fit success parameter is the parameter returned by the optimize.leastsq function.
Returns (height, x, y, sigma_x, sigma_y),
the gaussian parameters of a 2D distribution found by a fit.
error is computed based on the covariance matrix.
:param data: 2D distribution to fit
:type data: 2D array
:param fit_method: fit method to be used: 'gaussian_2d_fixedwidth', 'gaussian_2d_1width'
or 'gaussian_2d_2width'
:type fit_method: string
:param fit_sigma: standard deviation of the gaussian fit_function. Default is 1.555.
:type fit_sigma: float>0
:param verbose: verbosity parameter. Default=False
:type verbose: bool
:return: success_flag, fitparameters (height, x, y, sigma_x, sigma_y), error (=0 if success=0)
:rtype: bool, (floats), float
"""
# define set of parameters to be fitted.
params = gaussian_2d_moments(data, fit_method, fit_sigma)
X, Y = np.indices(data.shape)
# remove NaN in input data.
mask = ~np.isnan(data)
x = X[mask]
y = Y[mask]
data_without_nan = data[mask]
# define cost function to be minimized when fitting
# and run the least square minimization.
errorfunction = lambda p: np.ravel((fit_method(*p)(x, y, fit_sigma) - data_without_nan))
p, cov, infodict, mesg, success = optimize.leastsq(errorfunction, x0=params, full_output=1)
# check convergence. If so estimate fitting error.
if success and cov is not None:
# https://stackoverflow.com/questions/14581358/getting-standard-errors-on-fitted-parameters-using-the-optimize-leastsq-method-i
reduced_chi2 = np.sum(errorfunction(p) ** 2) / (data_without_nan.size - len(p))
errors = np.sqrt((cov * reduced_chi2).diagonal())
else:
errors = 0
if verbose == True:
logger.info(f'fit status: {success}')
logger.debug(infodict)
logger.debug(mesg)
logger.debug(p)
logger.debug(cov)
return success, p, errors