拟合二维高斯和，scipy.optimise.leastsq公司（回答：使用曲线拟合！）

Question

我想把高斯数的二维和拟合到这个数据上：

在一开始无法拟合和之后，我改为分别对每个峰值进行采样（image），并返回一个fit by find it's moments（基本上使用this code）。在

不幸的是，由于相邻峰的信号重叠，这会导致峰值位置测量不正确。下面是单独拟合的总和图。很明显，它们的峰都向中心倾斜。为了返回正确的峰值位置，我需要对此进行解释。在

我有一个绘制二维高斯包络函数（twoD_gaussian（））的工作代码，我通过优化.leatsq作为1D阵列使用努比·拉威尔以及一个适当的错误函数，但是这会导致无意义的输出。在

我尝试在求和中拟合一个峰值，得到以下错误输出：

如果您能给我提供什么样的建议，我将不胜感激，或者如果这不合适的话，可以选择其他方法。当然欢迎所有的意见！在

代码如下：

from scipy.optimize import leastsq
import numpy as np
import matplotlib.pyplot as plt

def twoD_Gaussian(amp0, x0, y0, amp1=13721, x1=356, y1=247, amp2=14753, x2=291,  y2=339, sigma=40):

    x0 = float(x0)
    y0 = float(y0)
    x1 = float(x1)
    y1 = float(y1)
    x2 = float(x2)
    y2 = float(y2)

    return lambda x, y:  (amp0*np.exp(-(((x0-x)/sigma)**2+((y0-y)/sigma)**2)/2))+(
                             amp1*np.exp(-(((x1-x)/sigma)**2+((y1-y)/sigma)**2)/2))+(
                             amp2*np.exp(-(((x2-x)/sigma)**2+((y2-y)/sigma)**2)/2))

def fitgaussian2D(x, y, data, params):

    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution found by a fit"""
    errorfunction = lambda p: np.ravel(twoD_Gaussian(*p)(*np.indices(np.shape(data))) - data)

    p, success = optimize.leastsq(errorfunction, params)
    return p     

# Create data indices
I = image   # Red channel of a scanned image, equivalent to the  1st image displayed in this post.
p = np.asarray(I).astype('float')
w,h = np.shape(I)
x, y = np.mgrid[0:h, 0:w]
xy = (x,y)

# scanned at 150 dpi = 5.91 dots per mm
dpmm = 5.905511811
plot_width = 40*dpmm

# create function indices
fdims = np.round(plot_width/2)  
xdims = (RC[0] - fdims, RC[0] + fdims)
ydims = (RC[1] - fdims, RC[1] + fdims)
fx = np.linspace(xdims[0], xdims[1], np.round(plot_width))
fy = np.linspace(ydims[0], ydims[1], np.round(plot_width))
fx,fy = np.meshgrid(fx,fy)

#Crop image for display
crp_data = image[xdims[0]:xdims[1], ydims[0]:ydims[1]]
z = crp_data    

# Parameters obtained from separate fits
Amplitudes = (13245, 13721, 15374)
px = (410, 356, 290)
py = (350, 247, 339)

initial_guess_sum = (Amp[0], px[0], py[0], Amp[1], px[1], py[1], Amp[2], px[2], py[2])
initial_guess_peak3 = (Amp[0], px[0], py[0]) # Try fitting single peak within sum


fitted_pars = fitgaussian2D(x, y, z, initial_guess_sum)
#fitted_pars = fitgaussian2D(x, y, z, initial_guess_peak3)

data_fitted= twoD_Gaussian(*fitted_pars)(fx,fy)
#data_fitted= twoD_Gaussian(*initial_guess_sum)(fx,fy)



fig = plt.figure(figsize=(10, 30))
ax = fig.add_subplot(111, aspect="equal")
#fig, ax = plt.subplots(1)
cb = ax.imshow(p, cmap=plt.cm.jet, origin='bottom',
    extent=(x.min(), x.max(), y.min(), y.max()))
ax.contour(fx, fy, data_fitted.reshape(fx.shape[0], fy.shape[1]), 4, colors='w')

ax.set_xlim(np.int(RC[0])-135, np.int(RC[0])+135)
ax.set_ylim(np.int(RC[1])+135, np.int(RC[1])-135)

#plt.colorbar(cb)
plt.show()

Answer 1

在放弃之前，我尝试了很多其他的方法，并再次尝试曲线拟合，尽管我对解析lambda函数有更多的了解。成功了。下面的示例输出和代码（仍有冗余）是为了子孙后代。在

def twoD_Gaussian(amp0, x0, y0, amp1=13721, x1=356, y1=247, amp2=14753, x2=291,  y2=339, sigma=40):

    x0 = float(x0)
    y0 = float(y0)
    x1 = float(x1)
    y1 = float(y1)
    x2 = float(x2)
    y2 = float(y2)

    return lambda x, y:  (amp0*np.exp(-(((x0-x)/sigma)**2+((y0-y)/sigma)**2)/2))+(
                             amp1*np.exp(-(((x1-x)/sigma)**2+((y1-y)/sigma)**2)/2))+(
                             amp2*np.exp(-(((x2-x)/sigma)**2+((y2-y)/sigma)**2)/2))

def twoD_GaussianCF(xy, amp0, x0, y0, amp1=13721, amp2=14753, x1=356, y1=247, x2=291,  y2=339, sigma_x=12, sigma_y=12):

    x0 = float(x0)
    y0 = float(y0)
    x1 = float(x1)
    y1 = float(y1)
    x2 = float(x2)
    y2 = float(y2)

    g = (amp0*np.exp(-(((x0-x)/sigma_x)**2+((y0-y)/sigma_y)**2)/2))+(
        amp1*np.exp(-(((x1-x)/sigma_x)**2+((y1-y)/sigma_y)**2)/2))+(
        amp2*np.exp(-(((x2-x)/sigma_x)**2+((y2-y)/sigma_y)**2)/2))

    return g.ravel()

# Create data indices
I = image   # Red channel of a scanned image, equivalent to the  1st image displayed in this post.
p = np.asarray(I).astype('float')
w,h = np.shape(I)
x, y = np.mgrid[0:h, 0:w]
xy = (x,y)

N_points = 3
display_width = 80

initial_guess_sum = (Amp[0], px[0], py[0], Amp[1], px[1], py[1], Amp[2], px[2], py[2])

popt, pcov = opt.curve_fit(twoD_GaussianCF, xy, np.ravel(p), p0=initial_guess_sum)

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

peaks = [(popt[1],popt[2]), (popt[5],popt[6]), (popt[7],popt[8])]

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, aspect="equal")
cb = ax.imshow(p, cmap=plt.cm.jet, origin='bottom',
    extent=(x.min(), x.max(), y.min(), y.max()))
ax.contour(x, y, data_fitted.reshape(x.shape[0], y.shape[1]), 20, colors='w')

ax.set_xlim(np.int(RC[0])-135, np.int(RC[0])+135)
ax.set_ylim(np.int(RC[1])+135, np.int(RC[1])-135)

for k in range(0,N_points):
    plt.plot(peaks[k][0],peaks[k][1],'bo',markersize=7)
plt.show()