适用于Python中特定坐标的2D高斯拟合 [英] 2D Gaussian Fit for intensities at certain coordinates in Python

查看:156
本文介绍了适用于Python中特定坐标的2D高斯拟合的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我有一组坐标(x,y,z(x,y)),它们描述了在坐标x,y处的强度(z).对于在不同坐标下的一定数量的这些强度,我需要拟合一个使均方误差最小的2D高斯函数. 数据以numpy矩阵表示,对于每个拟合会话,我将具有4、9、16或25个坐标.最终,我只需要获得MSE最小的高斯(x_0,y_0)的中心位置. 我发现的所有示例都使用scipy.optimize.curve_fit,但是它们拥有的输入数据是在整个网格上而不是几个坐标上. 任何帮助,将不胜感激.

解决方案

简介

有多种方法可以解决此问题.您可以使用非线性方法(例如scipy.optimize.curve_fit),但是它们会很慢并且不能保证收敛.您可以线性化问题(快速,独特的解决方案),但是分布尾部"中的任何噪声都会引起问题.实际上,您可以将一些技巧应用于此特定情况,以避免出现后一个问题.我将展示一些示例,但是我现在没有时间展示所有的技巧".

作为一个补充说明,一般的2D波斯语有6个参数,因此您将无法完全用4个点拟合事物.但是,听起来您可能假设x和y之间没有协方差,并且每个方向的方差都相同(即完美的圆"钟形曲线).如果是这样,那么您只需要四个参数.如果您知道高斯的幅度,则只需三个即可.但是,我将从通用解决方案开始,如果需要,您可以稍后对其进行简化.

目前,让我们集中讨论使用非线性方法(例如scipy.optimize.curve_fit)解决此问题.


二维高斯的一般方程为(直接来自维基百科):

其中:

在协方差矩阵上基本上是0.5,A是振幅, 并且(X₀,Y₀)是中心


生成简化的样本数据

让我们将上面的等式写出来:

import numpy as np
import matplotlib.pyplot as plt

def gauss2d(x, y, amp, x0, y0, a, b, c):
    inner = a * (x - x0)**2 
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

然后让我们生成一些示例数据.首先,我们将生成一些易于拟合的数据:

np.random.seed(1977) # For consistency
x, y = np.random.random((2, 10))
x0, y0 = 0.3, 0.7
amp, a, b, c = 1, 2, 3, 4

zobs = gauss2d(x, y, amp, x0, y0, a, b, c)

fig, ax = plt.subplots()
scat = ax.scatter(x, y, c=zobs, s=200)
fig.colorbar(scat)
plt.show()

请注意,我们没有添加任何噪声,并且分布的中心在我们拥有数据的范围内(即中心位于0.3、0.7处,并且x,y观测值的分散程度介于0和1之间).现在,让我们坚持下去,然后看看添加噪声并移动中心时会发生什么.


非线性拟合

首先,让我们使用scpy.optimize.curve_fit对高斯函数进行非线性最小二乘法拟合. (附带说明,您可以使用scipy.optimize中的其他一些功能来尝试使用精确的最小化算法.)

scipy.optimize函数期望的函数签名与我们上面最初写的函数签名略有不同.我们可以编写一个包装来翻译",但让我们改写gauss2d函数:

def gauss2d(xy, amp, x0, y0, a, b, c):
    x, y = xy
    inner = a * (x - x0)**2
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

我们所做的只是让函数期望独立变量(x&y)作为单个2xN数组.

现在,我们需要初步猜测高斯曲线的参数是什么.这是可选的(如果我没记错的话,默认值是全选),但是如果1、1不是特别接近高斯曲线的真实"中心,您可能会遇到收敛问题.因此,我们将使用观察到的最大z值的x和y值作为中心的起点.我将其余参数保留为1,但是如果您知道它们可能始终存在显着不同,请将其更改为更合理的值.

这是完整的独立示例:

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def main():
    x0, y0 = 0.3, 0.7
    amp, a, b, c = 1, 2, 3, 4
    true_params = [amp, x0, y0, a, b, c]
    xy, zobs = generate_example_data(10, true_params)
    x, y = xy

    i = zobs.argmax()
    guess = [1, x[i], y[i], 1, 1, 1]
    pred_params, uncert_cov = opt.curve_fit(gauss2d, xy, zobs, p0=guess)

    zpred = gauss2d(xy, *pred_params)
    print 'True parameters: ', true_params
    print 'Predicted params:', pred_params
    print 'Residual, RMS(obs - pred):', np.sqrt(np.mean((zobs - zpred)**2))

    plot(xy, zobs, pred_params)
    plt.show()

def gauss2d(xy, amp, x0, y0, a, b, c):
    x, y = xy
    inner = a * (x - x0)**2
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

def generate_example_data(num, params):
    np.random.seed(1977) # For consistency
    xy = np.random.random((2, num))

    zobs = gauss2d(xy, *params)
    return xy, zobs

def plot(xy, zobs, pred_params):
    x, y = xy
    yi, xi = np.mgrid[:1:30j, -.2:1.2:30j]
    xyi = np.vstack([xi.ravel(), yi.ravel()])

    zpred = gauss2d(xyi, *pred_params)
    zpred.shape = xi.shape

    fig, ax = plt.subplots()
    ax.scatter(x, y, c=zobs, s=200, vmin=zpred.min(), vmax=zpred.max())
    im = ax.imshow(zpred, extent=[xi.min(), xi.max(), yi.max(), yi.min()],
                   aspect='auto')
    fig.colorbar(im)
    ax.invert_yaxis()
    return fig

main()

在这种情况下,我们完全恢复了原始的"true"参数.

True parameters:  [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [ 1.   0.3  0.7  2.   3.   4. ]
Residual, RMS(obs - pred): 1.01560615193e-16

我们稍后将看到,情况并非总是如此...


添加噪音

让我们的观察结果有些混乱.我在这里所做的只是更改generate_example_data函数:

def generate_example_data(num, params):
    np.random.seed(1977) # For consistency
    xy = np.random.random((2, num))

    noise = np.random.normal(0, 0.3, num)
    zobs = gauss2d(xy, *params) + noise
    return xy, zobs

但是,结果看起来却大不相同:

就参数而言:

True parameters:  [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [  1.129    0.263   0.750   1.280   32.333   10.103  ]
Residual, RMS(obs - pred): 0.152444640098

预测的中心变化不​​大,但是bc参数已经发生了很大变化.

如果我们将函数的中心更改为稍微超出点分散范围的某个地方,则:

x0, y0 = -0.3, 1.1

由于存在噪音,我们将胡说八道! (它仍然可以正常运行而没有噪音.)

True parameters:  [1, -0.3, 1.1, 2, 3, 4]
Predicted params: [  0.546  -0.939   0.857  -0.488  44.069  -4.136]
Residual, RMS(obs - pred): 0.235664449826

当拟合衰减到零的函数时,这是一个常见问题. 尾巴"中的任何噪音都可能导致非常差的结果.有许多策略可以解决这个问题.最简单的方法之一是通过观察到的z值对反演进行加权.这是一维情况的示例:(专注于线性化问题) 解决方案

Introduction

There are multiple ways to approach this. You can use non-linear methods (e.g. scipy.optimize.curve_fit), but they'll be slow and aren't guaranteed to converge. You can linearize the problem (fast, unique solution), but any noise in the "tails" of the distribution will cause issues. There are actually a few tricks you can apply to this particular case to avoid the latter issue. I'll show some examples, but I don't have time to demonstrate all of the "tricks" right now.

Just as a side note, a general 2D guassian has 6 parameters, so you won't be able to fully fit things with 4 points. However, it sounds like you might be assuming that there's no covariance between x and y and that the variances are the same in each direction (i.e. a perfectly "round" bell curve). If that's the case, then you only need four parameters. If you know the amplitude of the guassian, you'll only need three. However, I'm going to start with the general solution, and you can simplify it later on, if you want to.

For the moment, let's focus on solving this problem using non-linear methods (e.g. scipy.optimize.curve_fit).


The general equation for a 2D guassian is (directly from wikipedia):

where:

is essentially 0.5 over the covariance matrix, A is the amplitude, and (X₀, Y₀) is the center


Generate simplified sample data

Let's write the equation above out:

import numpy as np
import matplotlib.pyplot as plt

def gauss2d(x, y, amp, x0, y0, a, b, c):
    inner = a * (x - x0)**2 
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

And then let's generate some example data. To start with, we'll generate some data that will be easy to fit:

np.random.seed(1977) # For consistency
x, y = np.random.random((2, 10))
x0, y0 = 0.3, 0.7
amp, a, b, c = 1, 2, 3, 4

zobs = gauss2d(x, y, amp, x0, y0, a, b, c)

fig, ax = plt.subplots()
scat = ax.scatter(x, y, c=zobs, s=200)
fig.colorbar(scat)
plt.show()

Note that we haven't added any noise, and the center of the distribution is within the range that we have data (i.e. center at 0.3, 0.7 and a scatter of x,y observations between 0 and 1). For the moment, let's stick with this, and then we'll see what happens when we add noise and shift the center.


Non-linear fitting

To start with, let's use scpy.optimize.curve_fit to preform a non-linear least-squares fit to the gaussian function. (On a side note, you can play around with the exact minimization algorithm by using some of the other functions in scipy.optimize.)

The scipy.optimize functions expect a slightly different function signature than the one we originally wrote above. We could write a wrapper to "translate", but let's just re-write the gauss2d function instead:

def gauss2d(xy, amp, x0, y0, a, b, c):
    x, y = xy
    inner = a * (x - x0)**2
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

All we did was have the function expect the independent variables (x & y) as a single 2xN array.

Now we need to make an initial guess at what the guassian curve's parameters actually are. This is optional (the default is all ones, if I recall correctly), but you're likely to have problems converging if 1, 1 is not particularly close to the "true" center of the gaussian curve. For that reason, we'll use the x and y values of our largest observed z-value as a starting point for the center. I'll leave the rest of the parameters as 1, but if you know that they're likely to consistently be significantly different, change them to something more reasonable.

Here's the full, stand-alone example:

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def main():
    x0, y0 = 0.3, 0.7
    amp, a, b, c = 1, 2, 3, 4
    true_params = [amp, x0, y0, a, b, c]
    xy, zobs = generate_example_data(10, true_params)
    x, y = xy

    i = zobs.argmax()
    guess = [1, x[i], y[i], 1, 1, 1]
    pred_params, uncert_cov = opt.curve_fit(gauss2d, xy, zobs, p0=guess)

    zpred = gauss2d(xy, *pred_params)
    print 'True parameters: ', true_params
    print 'Predicted params:', pred_params
    print 'Residual, RMS(obs - pred):', np.sqrt(np.mean((zobs - zpred)**2))

    plot(xy, zobs, pred_params)
    plt.show()

def gauss2d(xy, amp, x0, y0, a, b, c):
    x, y = xy
    inner = a * (x - x0)**2
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

def generate_example_data(num, params):
    np.random.seed(1977) # For consistency
    xy = np.random.random((2, num))

    zobs = gauss2d(xy, *params)
    return xy, zobs

def plot(xy, zobs, pred_params):
    x, y = xy
    yi, xi = np.mgrid[:1:30j, -.2:1.2:30j]
    xyi = np.vstack([xi.ravel(), yi.ravel()])

    zpred = gauss2d(xyi, *pred_params)
    zpred.shape = xi.shape

    fig, ax = plt.subplots()
    ax.scatter(x, y, c=zobs, s=200, vmin=zpred.min(), vmax=zpred.max())
    im = ax.imshow(zpred, extent=[xi.min(), xi.max(), yi.max(), yi.min()],
                   aspect='auto')
    fig.colorbar(im)
    ax.invert_yaxis()
    return fig

main()

In this case, we exactly(ish) recover our original "true" parameters.

True parameters:  [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [ 1.   0.3  0.7  2.   3.   4. ]
Residual, RMS(obs - pred): 1.01560615193e-16

As we'll see in a second, this won't always be the case...


Adding Noise

Let's add some noise to our observations. All I've done here is change the generate_example_data function:

def generate_example_data(num, params):
    np.random.seed(1977) # For consistency
    xy = np.random.random((2, num))

    noise = np.random.normal(0, 0.3, num)
    zobs = gauss2d(xy, *params) + noise
    return xy, zobs

However, the result looks quite different:

And as far as the parameters go:

True parameters:  [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [  1.129    0.263   0.750   1.280   32.333   10.103  ]
Residual, RMS(obs - pred): 0.152444640098

The predicted center hasn't changed much, but the b and c parameters have changed quite a bit.

If we change the center of the function to somewhere slightly outside of our scatter of points:

x0, y0 = -0.3, 1.1

We'll wind up with complete nonsense as a result in the presence of noise! (It still works correctly without noise.)

True parameters:  [1, -0.3, 1.1, 2, 3, 4]
Predicted params: [  0.546  -0.939   0.857  -0.488  44.069  -4.136]
Residual, RMS(obs - pred): 0.235664449826

This is a common problem when fitting a function that decays to zero. Any noise in the "tails" can result in a very poor result. There are a number of strategies to deal with this. One of the easiest is to weight the inversion by the observed z-values. Here's an example for the 1D case: (focusing on linearized the problem) How can I perform a least-squares fitting over multiple data sets fast? If I have time later, I'll add an example of this for the 2D case.

这篇关于适用于Python中特定坐标的2D高斯拟合的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆