I want to fit an intensity distribution function to 2D image data using scipy.optimize.curve_fit and can't locate the error in my code:
# Define doughnut beam intensity distribution function
def doughnut(x, y, x0, y0, A, FWHM):
'''2D intensity distribution function of doughnut beams (DOI: 10.1126/science.aak9913,
https://science.sciencemag.org/content/sci/suppl/2016/12/21/science.aak9913.DC1/Balzarotti_SM.pdf).
Parameters
----------
x, y : float
X and Y coordinates, orthogonal to beam axis
x0 : float
X offset
y0 : float
Y offset
A : float
Peak intensity
FWHM : float
Full width at half maximum
'''
return A*np.exp(1)*4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2*np.exp(-4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2)
# Read image file names
pathname = '/home/user/doughnut_beam/'
filenameList = [filename for filename in os.listdir(pathname)
if filename.endswith('.tif')]
# Open image files, fit doughnut beam intensity distribution function
for filename in filenameList:
img = Image.open(pathname + filename)
X, Y = img.size
xRange = np.arange(1, X+1)
yRange = np.arange(1, Y+1)
xGrid, yGrid = np.meshgrid(xRange, yRange)
xyGrid = np.vstack((xGrid.ravel(), yGrid.ravel())) # scipy.optimize.curve_fit requires 2xN-array
imgArray = np.array(img)
imgArrayFlat = imgArray.ravel() # Flatten 2D pixel data into 1D array for scipy.optimize.curve_fit
params_opt, params_cov = curve_fit(doughnut, xyGrid, imgArrayFlat)
This is the output from Jupyter Notebook:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-44-eaa3ebdb6469> in <module>()
17 imgArrayFlat = imgArray.ravel() # Flatten 2D pixel data into 1D array for scipy.optimize.curve_fit
18
---> 19 params_opt, params_cov = curve_fit(doughnut, xyGrid, imgArrayFlat)
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
749 # Remove full_output from kwargs, otherwise we're passing it in twice.
750 return_full = kwargs.pop('full_output', False)
--> 751 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
752 popt, pcov, infodict, errmsg, ier = res
753 cost = np.sum(infodict['fvec'] ** 2)
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
381 if not isinstance(args, tuple):
382 args = (args,)
--> 383 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
384 m = shape[0]
385 if n > m:
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
25 def _check_func(checker, argname, thefunc, x0, args, numinputs,
26 output_shape=None):
---> 27 res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
28 if (output_shape is not None) and (shape(res) != output_shape):
29 if (output_shape[0] != 1):
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in func_wrapped(params)
461 if transform is None:
462 def func_wrapped(params):
--> 463 return func(xdata, *params) - ydata
464 elif transform.ndim == 1:
465 def func_wrapped(params):
<ipython-input-43-3e0adae6fbe0> in doughnut(x, y, x0, y0, A, FWHM)
17 Full width at half maximum
18 '''
---> 19 return A*np.exp(1)*4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2*np.exp(-4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2)
ValueError: shapes (2,210) and (2,210) not aligned: 210 (dim 1) != 2 (dim 0)
UPDATE: For some reason, using numpy.dot to square the (offset) variables x+x0 and y+y0 in the function definition does not work. Simply changing to the ** operator results in the correct plot:
# UPDATED: Define doughnut beam intensity distribution function
def doughnut(x, y, x0, y0, A, FWHM):
'''2D intensity distribution function of doughnut beams (DOI: 10.1126/science.aak9913,
https://science.sciencemag.org/content/sci/suppl/2016/12/21/science.aak9913.DC1/Balzarotti_SM.pdf).
Parameters
----------
x, y : float
X and Y coordinates, orthogonal to beam axis
x0 : float
X offset
y0 : float
Y offset
A : float
Peak intensity
FWHM : float
Full width at half maximum
'''
return A*np.exp(1)*4*np.log(2)*((x+x0)**2 + (y+y0)**2)/FWHM**2*np.exp(-4*np.log(2)*((x+x0)**2 + (y+y0)**2)/FWHM**2)
fig = plt.figure()
ax = fig.gca(projection='3d')
# Make data
X = np.arange(-10, 10, 0.25)
Y = np.arange(-10, 10, 0.25)
X, Y = np.meshgrid(X, Y)
Z = doughnut(X, Y, x0=0, y0=0, A=1.5, FWHM=7)
# Plot the surface
surf = ax.plot_surface(X, Y, Z)
plt.show()
=> Plot
BUT: Now I'm getting a new error when trying to fit the data:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-61-eaa3ebdb6469> in <module>()
17 imgArrayFlat = imgArray.ravel() # Flatten 2D pixel data into 1D array for scipy.optimize.curve_fit
18
---> 19 params_opt, params_cov = curve_fit(doughnut, xyGrid, imgArrayFlat)
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
749 # Remove full_output from kwargs, otherwise we're passing it in twice.
750 return_full = kwargs.pop('full_output', False)
--> 751 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
752 popt, pcov, infodict, errmsg, ier = res
753 cost = np.sum(infodict['fvec'] ** 2)
/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
384 m = shape[0]
385 if n > m:
--> 386 raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
387 if epsfcn is None:
388 epsfcn = finfo(dtype).eps
TypeError: Improper input: N=5 must not exceed M=2
This should do the trick. Take a look at the for_fitting function to see how you can package everything in a way curve_fit will accept.
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
def doughnut(y, x, y0, x0, A, FWHM):
"""2D intensity distribution function of doughnut beams (DOI: 10.1126/science.aak9913,
https://science.sciencemag.org/content/sci/suppl/2016/12/21/science.aak9913.DC1/Balzarotti_SM.pdf).
Parameters
----------
y, x : float
X and Y coordinates, orthogonal to beam axis
y0 : float
Y offset
x0 : float
X offset
A : float
Peak intensity
FWHM : float
Full width at half maximum
"""
return (
A
* np.e
* 4
* np.log(2)
* ((x + x0) ** 2 + (y + y0) ** 2)
/ FWHM ** 2
* np.exp(-4 * np.log(2) * ((x + x0) ** 2 + (y + y0) ** 2) / FWHM ** 2)
)
fig0, (ax0, ax1, ax2) = plt.subplots(1, 3, sharex=True, sharey=True)
# Make data
X = np.arange(-10, 10, 0.25)
Y = np.arange(-10, 10, 0.25)
X, Y = np.meshgrid(X, Y)
true_params = (0, 0, 100, 7)
Z = doughnut(Y, X, *true_params)
# Plot the surface
ax0.matshow(Z, extent=(-10, 10, 10, -10))
ax0.set_title("Ground Truth")
def for_fitting(xdata, y0, x0, A, FWHM):
yy, xx = xdata
return doughnut(yy, xx, y0, x0, A, FWHM).ravel()
noisy_data = np.random.poisson(Z) + np.random.randn(*Z.shape)
ax1.matshow(noisy_data, extent=(-10, 10, 10, -10))
ax1.set_title("Noisy Data")
opt_params, cov = curve_fit(for_fitting, (Y, X), noisy_data.ravel(), p0=(0, 0, 10, 1))
print(opt_params)
fit_Z = doughnut(Y, X, *opt_params)
ax2.matshow(fit_Z, extent=(-10, 10, 10, -10))
ax2.set_title("Fit")
fig0.tight_layout()
plt.show()