使用 scipy.optimize.curve_fit - ValueError 和 minpack.error 拟合 2D 高斯函数

Posted

技术标签:

【中文标题】使用 scipy.optimize.curve_fit - ValueError 和 minpack.error 拟合 2D 高斯函数【英文标题】:Fitting a 2D Gaussian function using scipy.optimize.curve_fit - ValueError and minpack.error 【发布时间】:2014-03-01 05:49:01 【问题描述】:

我打算将二维高斯函数拟合到显示激光束的图像中,以获取其参数,例如 FWHM 和位置。到目前为止,我试图了解如何在 Python 中定义 2D 高斯函数以及如何将 x 和 y 变量传递给它。

我编写了一个小脚本,它定义了该函数,绘制它,为其添加一些噪音,然后尝试使用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)

这是我在使用 winpython 64-bit 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.

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

【问题讨论】:

由于高斯完全由其均值和协方差矩阵指定,您实际上根本不需要curve_fit。这是一个实现:code.google.com/p/agpy/source/browse/trunk/agpy/gaussfitter.py 【参考方案1】:

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)

顺便说一句:我不确定使用三角函数的参数化是否是最好的。例如,采用here 描述的那个可能在数值方面和大偏差下更加稳健。

【讨论】:

【参考方案2】:

twoD_Gaussian 的输出需要是一维的。您可以在最后一行的末尾添加一个.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()

您显然需要重塑输出以进行绘图,例如:

# 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()

像以前一样进行拟合:

# 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)

并绘制结果:

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()

【讨论】:

太好了,行得通!有趣的是,将数据实际拟合到高斯模型的方法比Zhenya 提出的code.google.com/p/agpy/source/browse/trunk/agpy/gaussfitter.py 更快。它实际上是我首先尝试的。现在拟合例程可以工作了,现在很容易将其扩展到其他功能。 虽然上面的代码确实提供了正确的输出参数,但对于非零 theta 值,结果将不正确地显示(轮廓会消失),除非将 origin='lower' 添加到 imshow 命令。这是因为 contour 和 imshow 具有不同的默认坐标系。 @Mike 非常正确 - 由于 OP 对示例参数的特殊选择,我没有发现这一点。我将在我的示例中添加额外的标志。 @Kokomoking 很好,是的,如果您使用完全正确的参数进行初始化,它肯定会收敛得更快:^) def twoD_Gaussian 似乎有一个 SyntaxError【参考方案3】:

为了稍微扩展 Dietrich 的回答,在使用 Python 3.4(在 Ubuntu 14.04 上)运行建议的解决方案时出现以下错误:

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

运行 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()

这样做的原因是,从 Python 3 开始,在将元组作为参数传递给函数时自动解包已被删除。有关更多信息,请参见此处:PEP 3113

【讨论】:

以上是关于使用 scipy.optimize.curve_fit - ValueError 和 minpack.error 拟合 2D 高斯函数的主要内容,如果未能解决你的问题,请参考以下文章

在使用加载数据流步骤的猪中,使用(使用 PigStorage)和不使用它有啥区别?

今目标使用教程 今目标任务使用篇

Qt静态编译时使用OpenSSL有三种方式(不使用,动态使用,静态使用,默认是动态使用)

MySQL db 在按日期排序时使用“使用位置;使用临时;使用文件排序”

使用“使用严格”作为“使用强”的备份

Kettle java脚本组件的使用说明(简单使用升级使用)