问题描述
我想将二维高斯总和拟合到此数据中:
最初未能将总和拟合到此之后,我改为分别对每个峰进行采样(图像) 并通过 find it's moment 返回合适值(主要使用 此代码)..>
不幸的是,由于相邻峰的重叠信号,这会导致峰位置测量不正确.下面是单独拟合的总和图.显然他们的峰都向中心倾斜.为了返回正确的峰值位置,我需要考虑到这一点.
我有绘制二维高斯包络函数 (twoD_Gaussian()) 的工作代码,我通过 optimize.leastsq 将其解析为使用 numpy.ravel 和适当的误差函数的一维数组,但这会导致废话输出.
我尝试在总和中拟合一个峰值并得到以下错误输出:
对于我可以尝试进行的工作或其他方法(如果这不合适)的任何建议,我将不胜感激.当然欢迎所有输入!
代码如下:
from scipy.optimize import leastsq将 numpy 导入为 np导入 matplotlib.pyplot 作为 pltdef twoD_Gaussian(amp0, x0, y0, amp1=13721, x1=356, y1=247, amp2=14753, x2=291, y2=339, sigma=40):x0 = 浮动(x0)y0 = 浮动(y0)x1 = 浮点数(x1)y1 = 浮动(y1)x2 = 浮点数(x2)y2 = 浮动(y2)返回 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):"""返回(高度,x,y,width_x,width_y)通过拟合找到的二维分布的高斯参数"""errorfunction = lambda p: np.ravel(twoD_Gaussian(*p)(*np.indices(np.shape(data))) - 数据)p, 成功 = optimize.leastsq(errorfunction, params)返回 p# 创建数据索引I = image # 扫描图像的红色通道,相当于这篇文章中显示的第一张图像.p = np.asarray(I).astype('float')w,h = np.shape(I)x, y = np.mgrid[0:h, 0:w]xy = (x,y)# 以 150 dpi = 5.91 点/毫米扫描dpmm = 5.905511811绘图宽度 = 40*dpmm# 创建函数索引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)#裁剪图像以供显示crp_data = 图像[xdims[0]:xdims[1], ydims[0]:ydims[1]]z = crp_data# 从单独拟合获得的参数振幅 = (13245, 13721, 15374)像素 = (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]) # 尝试在总和中拟合单峰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',范围=(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 函数有了更多的了解.有效.为了后代,下面的示例输出和代码(仍然有冗余).
def twoD_Gaussian(amp0, x0, y0, amp1=13721, x1=356, y1=247, amp2=14753, x2=291, y2=339, sigma=40):x0 = 浮动(x0)y0 = 浮动(y0)x1 = 浮点数(x1)y1 = 浮动(y1)x2 = 浮点数(x2)y2 = 浮动(y2)返回 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 = 浮动(x0)y0 = 浮动(y0)x1 = 浮点数(x1)y1 = 浮动(y1)x2 = 浮点数(x2)y2 = 浮动(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))返回 g.ravel()# 创建数据索引I = image # 扫描图像的红色通道,相当于这篇文章中显示的第一张图像.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显示宽度 = 80initial_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)峰值 = [(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',范围=(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)对于范围内的 k (0,N_points):plt.plot(peaks[k][0],peaks[k][1],'bo',markersize=7)plt.show()
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.
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!
Code below:
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()
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!)的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!