拟合二维高斯函数，使用scipy。优化。curve_fit - ValueError和minpack.error。

I intend to fit a 2D Gaussian function to images showing a laser beam to get its parameters like FWHM and position. So far I tried to understand how to define a 2D Gaussian function in Python and how to pass x and y variables to it.

我打算将二维高斯函数应用到显示激光束的图像上，以获得像FWHM和位置这样的参数。到目前为止，我试图理解如何在Python中定义二维高斯函数，以及如何将x和y变量传递给它。

I've written a little script which defines that function, plots it, adds some noise to it and then tries to fit it using curve_fit. Everything seems to work except the last step in which I try to fit my model function to the noisy data. Here is my code:

我写了一个小脚本，它定义了这个函数，绘制了它，给它添加了一些噪声，然后尝试使用curve_fit来匹配它。除了最后一步，我试着将模型函数与噪声数据相匹配，一切似乎都有效。这是我的代码:

import scipy.optimize as opt
import numpy as np
import pylab as plt


#define model function and pass independant variables x and y as a list
def twoD_Gaussian((x,y), amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
    xo = float(xo)
    yo = float(yo)    
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
    return offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo) + c*((y-yo)**2)))

# Create x and y indices
x = np.linspace(0, 200, 201)
y = np.linspace(0, 200, 201)
x,y = np.meshgrid(x, y)

#create data
data = twoD_Gaussian((x, y), 3, 100, 100, 20, 40, 0, 10)

# plot twoD_Gaussian data generated above
plt.figure()
plt.imshow(data)
plt.colorbar()

# add some noise to the data and try to fit the data generated beforehand
initial_guess = (3,100,100,20,40,0,10)

data_noisy = data + 0.2*np.random.normal(size=len(x))

popt, pcov = opt.curve_fit(twoD_Gaussian, (x,y), data_noisy, p0 = initial_guess)

Here is the error message I get when running the script using winpython 64-bit Python 2.7:

这里是使用winpython 64位Python 2.7运行脚本时得到的错误消息:

ValueError: object too deep for desired array
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Python\WinPython-64bit-2.7.6.2\python-2.7.6.amd64\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 540, in runfile
    execfile(filename, namespace)
  File "E:/Work Computer/Software/Python/Fitting scripts/2D Gaussian function fit/2D_Gaussian_LevMarq_v2.py", line 39, in <module>
    popt, pcov = opt.curve_fit(twoD_Gaussian, (x,y), data_noisy, p0 = initial_guess)
  File "C:\Python\WinPython-64bit-2.7.6.2\python-2.7.6.amd64\lib\site-packages\scipy\optimize\minpack.py", line 533, in curve_fit
    res = leastsq(func, p0, args=args, full_output=1, **kw)
  File "C:\Python\WinPython-64bit-2.7.6.2\python-2.7.6.amd64\lib\site-packages\scipy\optimize\minpack.py", line 378, in leastsq
    gtol, maxfev, epsfcn, factor, diag)
minpack.error: Result from function call is not a proper array of floats.

What is it that am I doing wrong? Is it how I pass the independent variables to the model function/curve_fit?

我做错了什么?它是如何将自变量传递给模型函数/curve_fit的?

3 个解决方案

#1

The output of twoD_Gaussian needs to be 1D. What you can do is add a .ravel() onto the end of the last line, like this:

twoD_Gaussian的输出必须是1D。你可以做的是在最后一行的末尾添加一个。ravel()，像这样:

def twoD_Gaussian((x, y), amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
    xo = float(xo)
    yo = float(yo)    
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
    g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo) 
                            + c*((y-yo)**2)))
    return g.ravel()

You'll obviously need to reshape the output for plotting, e.g:

很明显，你需要重新设计绘图的输出，例如:

# Create x and y indices
x = np.linspace(0, 200, 201)
y = np.linspace(0, 200, 201)
x, y = np.meshgrid(x, y)

#create data
data = twoD_Gaussian((x, y), 3, 100, 100, 20, 40, 0, 10)

# plot twoD_Gaussian data generated above
plt.figure()
plt.imshow(data.reshape(201, 201))
plt.colorbar()

Do the fitting as before:

和以前一样:

# add some noise to the data and try to fit the data generated beforehand
initial_guess = (3,100,100,20,40,0,10)

data_noisy = data + 0.2*np.random.normal(size=data.shape)

popt, pcov = opt.curve_fit(twoD_Gaussian, (x, y), data_noisy, p0=initial_guess)

And plot the results:

和情节的结果:

data_fitted = twoD_Gaussian((x, y), *popt)

fig, ax = plt.subplots(1, 1)
ax.hold(True)
ax.imshow(data_noisy.reshape(201, 201), cmap=plt.cm.jet, origin='bottom',
    extent=(x.min(), x.max(), y.min(), y.max()))
ax.contour(x, y, data_fitted.reshape(201, 201), 8, colors='w')
plt.show()

#2

curve_fit() wants to the dimension of xdata to be (2,n*m) and not (2,n,m). ydata should have shape (n*m) not (n,m) respectively. So you use ravel() to flatten your 2D arrays:

curve_fit()希望xdata的维度为(2,n*m)，而不是(2,n,m)。ydata应该有形状(n*m)，而不是(n,m)。所以你使用ravel()来压平你的二维数组:

xdata = np.vstack((xx.ravel(),yy.ravel()))
ydata = data_noisy.ravel()
popt, pcov = opt.curve_fit(twoD_Gaussian, xdata, ydata, p0=initial_guess)

By the way: I'm not sure if the parametrization with the trigonometric terms is the best one. E.g., taking the one described here might be a bit more robust under numerical aspects and large deviations.

顺便说一下，我不确定三角函数的参数化是否是最好的。举例来说，在数字方面和大的偏差下，这里所描述的可能会更健壮一些。

#3

To expand on Dietrich's answer a bit, I got the following error when running the suggested solution with Python 3.4 (on Ubuntu 14.04):

为了扩展Dietrich的答案，在使用Python 3.4(在Ubuntu 14.04上)运行建议的解决方案时，我得到了以下错误:

def twoD_Gaussian((x, y), amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
                  ^
SyntaxError: invalid syntax

Running 2to3 suggested the following simple fix:

运行2to3提出了以下简单的解决方案:

def twoD_Gaussian(xdata_tuple, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
    (x, y) = xdata_tuple                                                        
    xo = float(xo)                                                              
    yo = float(yo)                                                              
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)   
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)    
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)   
    g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)         
                        + c*((y-yo)**2)))                                   
    return g.ravel()

The reason for this is that automatic tuple unpacking when it is passed to a function as a parameter has been removed as of Python 3. For more information see here: PEP 3113

这样做的原因是，当将其传递给函数时，自动元组解压缩已经被删除为Python 3。更多信息请看这里:PEP 3113。