python - 将模型二维拟合到 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’ : Trust Region Reflective algorithm, particularly suitable for large sparse problems with bounds. Generally robust method.
• ‘dogbox’ : dogleg algorithm with rectangular trust regions, typical use case is small problems with bounds. Not recommended for problems with rank-deficient Jacobian.
• ‘lm’ : Levenberg-Marquardt algorithm as implemented in MINPACK. Doesn’t handle bounds and sparse Jacobians. Usually the most efficient method for small unconstrained problems.

curve_fit will use different methods for bounded and unbounded cases

Default is ‘lm’ for unconstrained problems and ‘trf’ if bounds are provided

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