模型与 Python 中图像的 2D 拟合

Question

我想用 Python 中的图像拟合模型（这里是 2D 高斯，但可能是其他东西）。

尝试使用scipy.optimize.curve_fit 我有一些问题。见下文。

让我们从一些函数开始：

import numpy as np
from scipy.optimize import curve_fit
from scipy.signal import argrelmax

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.patches import Circle

from tifffile import TiffFile

# 2D Gaussian model
def func(xy, x0, y0, sigma, H):

    x, y = xy

    A = 1 / (2 * sigma**2)
    I = H * np.exp(-A * ( (x - x0)**2 + (y - y0)**2))
    return I

# Generate 2D gaussian
def generate(x0, y0, sigma, H):

    x = np.arange(0, max(x0, y0) * 2 + sigma, 1)
    y = np.arange(0, max(x0, y0) * 2 + sigma, 1)
    xx, yy = np.meshgrid(x, y)

    I = func((xx, yy), x0=x0, y0=y0, sigma=sigma, H=H)

    return xx, yy, I

def fit(image, with_bounds):

    # Prepare fitting
    x = np.arange(0, image.shape[1], 1)
    y = np.arange(0, image.shape[0], 1)
    xx, yy = np.meshgrid(x, y)

    # Guess intial parameters
    x0 = int(image.shape[0]) # Middle of the image
    y0 = int(image.shape[1]) # Middle of the image
    sigma = max(*image.shape) * 0.1 # 10% of the image
    H = np.max(image) # Maximum value of the image
    initial_guess = [x0, y0, sigma, H]

    # Constraints of the parameters
    if with_bounds:
        lower = [0, 0, 0, 0]
        upper = [image.shape[0], image.shape[1], max(*image.shape), image.max() * 2]
        bounds = [lower, upper]
    else:
        bounds = [-np.inf, np.inf]

    pred_params, uncert_cov = curve_fit(func, (xx.ravel(), yy.ravel()), image.ravel(),
                                        p0=initial_guess, bounds=bounds)

    # Get residual
    predictions = func((xx, yy), *pred_params)
    rms = np.sqrt(np.mean((image.ravel() - predictions.ravel())**2))

    print("True params : ", true_parameters)
    print("Predicted params : ", pred_params)
    print("Residual : ", rms)

    return pred_params

def plot(image, params):

    fig, ax = plt.subplots()
    ax.imshow(image, cmap=plt.cm.BrBG, interpolation='nearest', origin='lower')

    ax.scatter(params[0], params[1], s=100, c="red", marker="x")

    circle = Circle((params[0], params[1]), params[2], facecolor='none',
            edgecolor="red", linewidth=1, alpha=0.8)
    ax.add_patch(circle)

# Simulate and fit model
true_parameters = [50, 60, 10, 500]
xx, yy, image = generate(*true_parameters)

# The fit performs well without bounds
params = fit(image, with_bounds=False)
plot(image, params)

输出：

True params :  [50, 60, 10, 500]
Predicted params :  [  50.   60.   10.  500.]
Residual :  0.0

现在，如果我们对边界（或约束）进行同样的拟合。

# The fit is really bad with bounds
params = fit(image, with_bounds=True)
plot(image, params)

输出：

True params :  [50, 60, 10, 500]
Predicted params :  [ 130.          130.            0.72018729    1.44948159]
Residual :  68.1713019773

为什么当我添加边界时拟合效果不佳？

---

现在还有一件事我不明白。为什么这种拟合在应用于真实数据时并不稳健？见下文。

# Load some real data
image = TiffFile("../data/spot.tif").asarray()
plt.imshow(image, aspect='equal', origin='lower', interpolation="none", cmap=plt.cm.BrBG)
plt.colorbar()

# Fit is not possible without bounds
params = fit(image, with_bounds=False)
plot(image, params)

输出：

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-14-3187b53d622d> in <module>()
      1 # Fit is not possible without bounds
----> 2 params = fit(image, with_bounds=False)
      3 plot(image, params)

<ipython-input-11-f14c9dec72f2> in fit(image, with_bounds)
     54 
     55     pred_params, uncert_cov = curve_fit(func, (xx.ravel(), yy.ravel()), image.ravel(),
---> 56                                         p0=initial_guess, bounds=bounds)
     57 
     58     # Get residual

/home/hadim/local/conda/envs/ws/lib/python3.5/site-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, **kwargs)
    653         cost = np.sum(infodict['fvec'] ** 2)
    654         if ier not in [1, 2, 3, 4]:
--> 655             raise RuntimeError("Optimal parameters not found: " + errmsg)
    656     else:
    657         res = least_squares(func, p0, args=args, bounds=bounds, method=method,

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 1000.

和

# Fit works but is not accurate at all with bounds
params = fit(image, with_bounds=True)
plot(image, params)

输出：

True params :  [50, 60, 10, 500]
Predicted params :  [   19.31770886    10.52153346    37.          1296.22524248]
Residual :  83.1944464761

Answer 1

几件事，首先，你的初始参数x0和y0是错误的，它们不在图像的中间，而是在边界，应该是

x0 = int(image.shape[0])/2 # Middle of the image
y0 = int(image.shape[1])/2 # Middle of the image

将它们放在图像的边界可能会在受限情况下产生一些问题，因为它不会给它在某些方向上移动的空间。这是我的猜测，取决于拟合方法。

同样说到方法，curve_fit 可以使用三个中的任何一个：scipy least_squares documentation 中的lm、trf 和dogbox：

• ‘trf’：信任区域反射算法，特别适用于有界的大稀疏问题。一般稳健的方法。
• ‘dogbox’：具有矩形信任区域的 dogleg 算法，典型用例是有边界的小问题。不推荐用于秩不足的雅可比行列式问题。
• ‘lm’：在MINPACK 中实现的Levenberg-Marquardt 算法。不处理边界和稀疏雅可比行列式。通常是解决小型无约束问题的最有效方法。

curve_fitwill use different methods for bounded and unbounded cases

对于不受约束的问题，默认为“lm”，如果提供了边界，则默认为“trf”

所以我建议定义一个要使用的方法，在纠正初始参数后，我使用您的示例使用trf 和dogbox 得到了很好的结果，但是您应该检查哪种方法更适合您的真实数据。

Answer 2

我写了一个lightweight class 来做到这一点。边界没有很好地实现，但可以根据您的需要进行更改。

您在这里遇到三个主要问题：

1. 您没有在数据中建模偏移。从您的图片来看，背景大约是 1300-1400 计数。
2. 最初的猜测参数不是很好的估计（@Noel 指出）
3. 您拟合的数据比需要的多得多，我建议选择峰值附近的区域。如果您对初始参数做出更好的猜测，则可以将拟合限制在以猜测 x0 和 y0 为中心的窗口内。

解决问题1有两种方法：

• 估计背景并从数据中减去它（数据的median 或mode 是一个不错的选择，my class 使用RANSAC 提供的RANSAC 算法在更复杂的方式）
• 或者您可以直接为背景建模。

对于问题2，你可以使用skimage的blob detection algorithms。我还编写了另一个 class 包装 skimage 的 DOG 算法，以使这更容易。一旦你解决了问题 2，问题 3 也解决了。