拟合2D高斯和，scipy.optimise.leastsq（回答：使用curve_fit！） [英] Fitting 2D sum of gaussians, scipy.optimise.leastsq (Ans: Use curve_fit!)
问题描述
我想对这些数据拟合2维高斯：
I want to fit an 2D sum of gaussians to this data:
在最初不为此求和后，我改为分别对每个峰进行采样（图片），并通过查找当前时刻返回拟合值（基本上使用此代码）。
After failing at fitting a sum to this initially I instead sampled each peak separately (image) and returned a fit by find it's moments (essentially using this code).
不幸的是，由于信号重叠，导致峰位置测量不正确相邻的山峰。下面是各个拟合值总和的图。显然，他们的高峰都向中心倾斜。为了返回正确的峰值位置，我需要考虑这一点。
Unfortunately, this results in an incorrect peak position measurement, due to the overlapping signal of the neighbouring peaks. Below is a plot of the sum of the separate fits. Obviously their peak all lean toward the centre. I need to account for this in order to return the correct peak position.
3个独立的高斯拟合之和>
我有绘制二维高斯包络函数（twoD_Gaussian（））的工作代码，通过使用numpy.ravel和适当的误差函数，通过optimize.leastsq将其解析为一维数组，但是这会导致无意义的输出。
I've got working code which plots a 2D gaussian envelope function (twoD_Gaussian()), and I parse this through optimize.leastsq as a 1D array using numpy.ravel and an appropriate error function, however this results in a nonsense output.
我尝试在其中拟合单个峰总和并获得以下错误输出：
I tried fitting a single peak within the sum and get the following erroneous output:
对于我可以尝试执行此操作的任何建议，或者如果不是这样的替代方法，我将不胜感激。不适合。当然欢迎所有输入！
I'd appreciate any advice on what i could try to make this work, or alternative approaches if this isn't appropriate. All input welcomed of course!
以下代码：
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()
推荐答案
在放弃并再次尝试curve_fit之前，我尝试了许多其他方法，尽管掌握了更多有关解析lambda函数的知识。有效。为了后代，下面的示例输出和代码（仍带有冗余）。
I tried any number of other things before giving up and trying curve_fit again, albeit with more knowledge of parsing lambda functions. It worked. Example output and code below (still with redundancies) for the sake of posterity.
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()
这篇关于拟合2D高斯和，scipy.optimise.leastsq（回答：使用curve_fit！）的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
