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

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可以使用三种：lm，trf和dogbox从scipy least_squares documentation：

• “TRF”：信赖域反射算法，特别适合于大型稀疏问题界限。一般健壮的方法。
• '狗箱'：具有矩形信任区域的狗腿算法，典型的用例是带边界的小问题。不建议用于排位不足的雅可比行列式问题。
• '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算法sklearn以更复杂的方式进行评估）
• 或者您可以直接对背景进行建模。

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