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 

する必要がありますそれらを画像の境界に置くと、制約された場合にいくつかの方向に移動する余裕がないため、問題が発生する可能性があります。これは私の推測であり、フィッティング方法に依存します。 scipy least_squares documentationからlm、trfとdogbox：

もcurve_fit三のいずれかを使用することができ、方法と言えば

• 'TRF'：信頼領域Reflectiveアルゴリズム、境界で大きな疎な問題のために特に適して。一般的に頑強な方法。
• 'dogbox'：矩形の信頼領域を持つdoglegアルゴリズムで、典型的な使用例は境界の小さな問題です。ランクが不足しているヤコビ行列の問題にはお勧めできません。
• 'lm'：MINPACKで実装されているLevenberg-Marquardtアルゴリズム。境界線やスパースヤコビアンを扱いません。通常、小さな拘束されていない問題の最も効率的な方法です。

curve_fitwill use different methods for bounded and unbounded cases

境界が

を提供している場合は、デフォルトではだから私は、使用するメソッドを定義し、私はいいまし示唆制約のない問題と 'TRF' の 'LM' です最初のパラメータを修正した後にtrfとdogboxのいずれかの結果が表示されますが、実際のデータでどの方法が効果的かを確認する必要があります。

Answer 2

私はこれを正確に行うためにlightweight classと書いています。境界は非常にうまく実装されていませんが、これはニーズに合わせて変更できます。

あなたがここに持っている三つの主要な問題があります：

1. はあなたがあなたのデータでオフセットをモデリングしていません。あなたのイメージを見て、背景は約1300-1400カウントです。
2. 最初の推測パラメータは良い推定値ではありません（@Noelによって指摘されているように）
3. あなたは必要以上に多くのデータを扱っていますが、ピーク付近の領域を選択することをお勧めします。あなたが最初のパラメータのより良い推測をするならば、あなたの推測を中心にしたウィンドウにあなたのフィットを制限することができますx0とy0。

あり、問題1を解決するための2つの方法があります。

• 推定値が背景と（データのmedianまたはmodeは良い賭けである、あなたのデータからそれを引く、my classはによって提供さRANSACアルゴリズムを使用していますsklearnこれをより洗練された方法で見積もる）
• また、背景を直接モデル化することもできます。

問題2の場合、skimageのblob detection algorithmsを使用できます。私はまたこれを簡単にするためにskimageのDOGアルゴリズムをラップする別のclassを書きました。いったん問題2を解決すれば、問題3も解決されます。