使用curve_fit进行多元高斯拟合的问题

问题描述

我一直在尝试拟合 3D 高斯分布，我希望有人能看到我犯的一些愚蠢的错误，因为我花了几个小时调试无济于事。

我有一个 3d 图像存储在一个名为“data”的数组中，其中 data[x,y,z] 给出点 (x,z) 处的灰度强度。我知道这个 3d 图像遵循 3D 高斯分布，在图像中心附近有一个峰值，但我对振幅和传播感兴趣。我正在尝试将此 3d 数组拟合为

形式的高斯我在 Python 中的函数是：

def gaussian_3d(X,A,x0,y0,z0,sigx,sigy,sigz,offset):
    x,z = X
    return offset + A*np.exp(-(x-x0)**2/(2*sigx**2) - \
            (y-y0)**2/(2*sigy**2) - (z-z0)**2/(2*sigz**2))

我这样做的方式如下：如果我的图像大小为 3 x 4 x 5，那么我创建一个网格 (0...2) x (0...3) x (0. ..4），然后尝试将强度值拟合到上面的函数中。我的代码如下所示：

def fit_gauss_3d(data):
    dim = data.shape
    
    # Step 1: set up meshgrid
    x,z = np.arange(0,dim[0]),np.arange(0,dim[1]),dim[2])
    X,Y,Z = np.meshgrid(x,z)
    data_in = np.vstack((X.ravel(),Y.ravel(),Z.ravel()))
    data_out = data.ravel()

    # Step 2: make good guess of the center "peak" point of the gaussian (x0,z0)
    #         by using slices along the middle and finding the position of the maxes
    mid1,mid2,mid3 = dim[0]//2,dim[1]//2,dim[2]//2
    x0,z0 = np.argmax(data[:,mid3]),np.argmax(data[mid1,:,:])

    # Step 3: Set lower/upper bounds for parameter search
    delta = 0.5  # I am saying that the fit peak must be within +/- 0.5 of the initial guess
    p0 = (data_max + 0.05,0.9,0.05)
    #         Note: I kNow that sigmas are between 0.7 and 2.5,and offset is between 0 and 5
    lower_bound = [data_max * 0.9,x0 - delta,y0 - delta,z0 - delta,0.7,0]
    upper_bound = [data_max*1.1 + 0.1,x0 + delta,y0 + delta,z0 + delta,2.5,5]

    # Step 4: Fit
    p_param,p_cov = opt.curve_fit(gauss_3d,data_in,data_out,p0=p0,maxfev=50000,bounds=(lower_bound,upper_bound))

    return p_param


def predict_gauss_3d(params,dims):
    x = np.arange(0,dims[0])
    y = np.arange(0,dims[1])
    z = np.arange(0,dims[2])
    XX,YY,ZZ = np.meshgrid(x,z)
    X = np.vstack((XX.ravel(),YY.ravel(),ZZ.ravel()))
    return gaussian_3d(X,*params).reshape(dims)

def plot_results(orig,sec):
    ''' Plot original and second fitted image'''
    mid1,dim[2]//2

    fig = plt.figure()
    ax1 = fig.add_subplot(3,1,1)
    ax1.plot(orig[:,mid3],label='orig')
    ax1.plot(sec[:,label='fitted')
    ax1.legend(loc="upper left")

    ax2 = fig.add_subplot(3,2)
    ax2.plot(orig[mid1,label='orig')
    ax2.plot(sec[mid1,label='fitted')
    ax2.legend(loc="upper left")

    ax3 = fig.add_subplot(3,3)
    ax3.plot(orig[mid1,mid2],label='orig')
    ax3.plot(sec[mid1,label='fitted')
    ax3.legend(loc="upper left")
    plt.tight_layout()
    plt.show()

我绘制了沿中间轴拟合的投影。第一张图片是改变 x 并将 y、z 保持在它们的中点，第二张是改变 y 并将 x、z 保持在它们的中点，依此类推。我的一些搭配是合理的，像这样：

虽然大多数都非常糟糕，甚至看起来都不是高斯！对于下图，它选择了以下参数：

。显然，我要么绘制错误，要么拟合错误。有人可以帮我吗？我的网格是否以某种方式搞砸了？

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

curve-fitting gaussian python scientific-computing scipy scipy