"""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
Created on 20 mars 2018
@author: Catalano Camille, APC
"""
import numpy as np
from scipy import optimize
[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
"""
sigma = 1.155
return lambda x,y: 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
"""
sigma = 1.155
return lambda x,y: 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,fitmethod):
"""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 distrbution to evaluate
:type data: float 2D-array
:param fitmethod: fit method that will be used: 'gaussian_2d_fixedwidth', 'gaussian_2d_1width' or 'gaussian_2d_2width'
:type fitmethod: string
:return: (height, x, y), (height, x, y, sigma_x, sigma_y) or (height, x, y, sigma) depending on fitmethod
:rtype: floats
"""
#general case
#total = data.sum()
#X, Y = np.indices(data.shape)
#x = (X*data).sum()/total
#y = (Y*data).sum()/total
#col = data[:, int(y)]
#width_x = np.sqrt(np.abs((np.arange(col.size)-y)**2*col).sum()/np.abs(col.sum()))
#row = data[int(x), :]
#width_y = np.sqrt(np.abs((np.arange(row.size)-x)**2*row).sum()/np.abs(row.sum()))
#
# simplification in our case
height = np.nanmax(data)
x,y=np.where(data==np.nanmax(data))
x = x[0]
y = y[0]
# largeur de départ fixée à 1 et symétrique (source ponctuelle = 1.155 pix)
sigma = 1.155
offset = 0
if fitmethod.__name__ == 'gaussian_2d_fixedwidth' :
return height, x, y
elif fitmethod.__name__ == 'gaussian_2d_1width' :
return height, x, y, sigma
elif fitmethod.__name__ == 'gaussian_2d_2width' :
return height, x, y, sigma, sigma
elif fitmethod.__name__ == 'gaussian_2d_fixedwidth_offset' :
return height, x, y, offset
elif fitmethod.__name__ == 'gaussian_2d_1width_offset' :
return height, x, y, sigma, offset
elif fitmethod.__name__ == 'gaussian_2d_2width_offset' :
return height, x, y, sigma, sigma, offset
[docs]def fit_gaussian_2d(data,fitmethod,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 fitmethod: fit method to be used: 'gaussian_2d_fixedwidth', 'gaussian_2d_1width' or 'gaussian_2d_2width'
:type fitmethod: string
: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
"""
params = gaussian_2d_moments(data,fitmethod)
X, Y = np.indices(data.shape)
#remove nan
mask = ~np.isnan(data)
x = X[mask]
y = Y[mask]
data_without_nan = data[mask]
errorfunction = lambda p: np.ravel((fitmethod(*p)(x, y) - data_without_nan))
p, cov, infodict, mesg, success = optimize.leastsq(errorfunction, params,full_output=1)
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 :
print('Success = %d'%(success))
print(p)
print(cov)
return success,p,errors