2D Fit eines Modells zu einem Bild in Python

Question

Ich möchte ein Modell (hier ein 2D Gaussian, aber es könnte etwas anderes sein) mit einem Bild in Python passen.

Ich versuche zu verwenden scipy.optimize.curve_fit Ich habe ein paar Fragen. Siehe unten.

Beginnen wir mit einigen Funktionen:

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) 

Ausgang:

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

Nun, wenn wir die gleiche Passform mit Grenzen (oder Einschränkungen) tun.

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

Ausgang:

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

Warum die Passform nicht gut durchführen, wenn ich Grenzen hinzufügen?

---

Jetzt eine andere Sache, die ich nicht verstehe. Warum ist diese Anpassung nicht robust, wenn sie auf echte Daten angewendet wird? Siehe unten.

# 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) 

Ausgang:

--------------------------------------------------------------------------- 
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. 

Und

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

Ausgang:

Answer 1

Paar Dinge, erstens, Ihre Anfangsparameter x0 und y0 sind falsch, sie sind nicht in der Mitte des Bildes, aber an der Grenze, sie sein sollten

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

Wenn Sie sie an der Grenze des Bildes haben, kann dies im eingeschränkten Fall zu Problemen führen, da sie in einigen Richtungen nicht frei beweglich ist. Das ist meine Spekulation und hängt von der Anpassungsmethode ab.

auch von Methoden gesprochen, curve_fit kann entweder von drei verwenden: lm, trf und dogbox vom scipy least_squares documentation:

• 'trf': Vertrauen Region Reflective Algorithmus, besonders geeignet für große spärliche Probleme mit Grenzen . In der Regel robuste Methode.
• 'Dogbox': Dogleg-Algorithmus mit rechteckigen Vertrauensbereiche, typischer Anwendungsfall ist kleine Probleme mit Grenzen. Nicht empfohlen für Probleme mit rank-defizienten Jacobi.
• 'lm': Levenberg-Marquardt-Algorithmus wie in MINPACK implementiert. Behandelt nicht Grenzen und spärliche Jacobians. In der Regel die effizienteste Methode für kleine unbeschränkte Probleme.

curve_fitwill use different methods for bounded and unbounded cases

Standard ist 'lm' für unbeschränkte Probleme und 'trf', wenn Grenzen vorgesehen sind

So empfehle definiere ich eine Methode zu verwenden, ich habe gut Ergebnisse mit entweder trf und dogbox mit Ihrem Beispiel nach Korrektur der Anfangsparameter, aber Sie sollten überprüfen, welche Methode besser mit Ihren realen Daten funktioniert.

Answer 2

Ich schrieb eine lightweight class um genau dies zu tun. Grenzen werden nicht sehr gut implementiert, aber dies könnte geändert werden, um Ihren Bedürfnissen zu entsprechen.

ist es drei Hauptprobleme, Sie hier zu haben sind:

1. Dir Modellierung keine Offset in Ihren Daten. Wenn Sie Ihr Bild betrachten, beträgt der Hintergrund 1300 bis 1400 Zählungen.
2. Die ersten Schätzwerte sind keine guten Schätzungen (wie von @Noel gezeigt)
3. Sie passen viel mehr Daten als nötig, würde ich empfehlen, die Region um den Gipfel zu wählen. Wenn Sie die Anfangsparameter besser schätzen, können Sie Ihre Anpassung auf ein Fenster beschränken, das auf Ihrer Schätzung basiert. x0 und y0.

Es gibt zwei Möglichkeiten der Problemlösung 1:

• Schätzen Sie den Hintergrund und subtrahieren sie von Ihren Daten (die median oder mode der Daten ist eine gute Wette, my class verwendet den RANSAC Algorithmus zur Verfügung gestellt von sklearn, um dies auf eine anspruchsvollere Weise zu schätzen)
• Oder Sie könnten den Hintergrund direkt modellieren.

Für Problem 2. Sie können skimage 's blob detection algorithms verwenden. Ich habe auch eine andere class geschrieben, die skimage 's DOG-Algorithmus umschließt, um dies zu erleichtern. Sobald Sie das Problem 2 gelöst haben, wird auch das Problem 3 gelöst.