2014-10-22 184 views
2

我试图拟合二维高斯到一些灰度图像数据,这是由一个二维数组给出的。 lmfit库实现了一个易于使用的Model类,它应该能够做到这一点。 不幸的是,文档(http://lmfit.github.io/lmfit-py/model.html)只提供了一维拟合的例子。对于我的情况,我简单地用两个独立变量构造lmfit模型。Python lmfit:拟合2D模型

下面的代码似乎对我有效,但会导致scipy抛出“minpack.error:函数调用的结果不是一个适当的浮点数组。

Tom总结:如何将2D(x1,x2) - >(y)数据输入到lmfit模型。

以下是我的方法: 所有内容都包含在GaussianFit2D类中,但以下是重要部分: 这就是高斯函数。文档说有关用户定义的函数

Of course, the model function will have to return an array that will be the same size as the data being modeled. Generally this is handled by also specifying one or more independent variables.

我真的不明白这是什么意思应该,因为对于给定的值X1,X2的唯一合理的结果是标量值。

def _function(self, x1, x2, amp, wid, cen1, cen2): 
    val = (amp/(np.sqrt(2*np.pi)*wid)) * np.exp(-((x1-cen1)**2+(x2-cen2)**2)/(2*wid**2)) 
    return val 

这里生成的模型:

def _buildModel(self, **kwargs): 
    model = lmfit.Model(self._function, independent_vars=["x1", "x2"], 
         param_names=["amp", "wid", "cen1", "cen2"]) 
    return model 

这是获取数据的功能,构建模型,而params,并呼吁lmfit符合():终于在这里

def fit(self, data, freeX, **kwargs): 
    freeX = np.asarray(freeX, float) 
    model = self._buildModel(**kwargs) 
    params = self._generateModelParams(model, **kwargs) 

    model.fit(data, x1=freeX[0], x2=freeX[1], params=params) 

ANF这个适合的功能被称为:

data = np.asarray(img, float) 
    gaussFit = GaussianFit2D() 
    x1 = np.arange(len(img[0, :])) 
    x2 = np.arange(len(img[:, 0])) 
    fit = gaussFit.fit(data, [x1, x2]) 
+0

你需要'lmfit',或其他工具('curve_fit'或SciPy的'leastsq')也未尝不可? – Evert 2014-10-22 07:57:01

+0

lmfit很迷人,因为它很容易为参数,约束等添加初始值......如果我不管理它,那么scipy可能是第二次尝试。 – m0e 2014-10-23 06:52:19

回答

0

好吧,与开发人员一起写下来,并从他们那里得到答案(在此感谢马特)。

其基本思想是将所有输入平坦化为1D数据,从lmfit隐藏> 1维输入。 以下是你如何做到的。 修改你的函数:

def function(self, x1, x2): 
     return (x1+x2).flatten() 

拼合你的二维输入数组要适合:

... 
data = data.flatten() 
... 

修改两个1D X变量,例如,你有他们的任意组合:

... 
x1n = [] 
x2n = [] 
    for i in x1: 
     for j in x2: 
       x1n.append(i) 
       x2n.append(j) 
x1n = np.asarray(x1n) 
x2n = np.asarray(x2n) 
... 

并将任何东西扔进钳工:

model.fit(data, x1=x1n, x2=x2n, params=params) 
0

这是一个供您参考的例子,希望它能帮助您。

import numpy 
from lmfit import Model 

def gaussian(x, cenu, cenv, wid): 

    u = x[:, 0] 
    v = x[:, 1] 
    return (1/(2*numpy.pi*wid**2)) * numpy.exp(-(u-cenu)**2/(2*wid**2)-(v-cenv)**2/(2*wid**2)) 

data = numpy.empty((25,3)) 
x = numpy.arange(-2,3,1) 
y = numpy.arange(-2,3,1) 
xx, yy = numpy.meshgrid(x, y) 
data[:,0] = xx.flatten() 
data[:,1] = yy.flatten() 


data[:, 2]= gaussian(data[:,0:2],0,0,0.5) 

print 'xx\n', xx 
print 'yy\n',yy 
print 'data to be fit\n', data[:, 2] 

cu = 0.9 
cv = 0.5 
wid = 1 
gmod = Model(gaussian) 

gmod.set_param_hint('cenu', value=cu, min=cu-2, max=cu+2) 
gmod.set_param_hint('cenv', value=cv, min=cv -2, max=cv+2) 
gmod.set_param_hint('wid', value=wid, min=0.1, max=5) 
params = gmod.make_params() 
result = gmod.fit(data[:, 2], x=data[:, 0:2], params=params) 
print result.fit_report(min_correl=0.25) 
print result.best_values 
print result.best_fit 
+0

解释你的答案也会有帮助。 – Fred 2015-04-16 16:11:35