Quero ajustar uma soma 2D de gaussianos a esses dados:

Depois de falhar em ajustar uma soma a isso, inicialmente eu experimentei cada pico separadamente (imagem) e retornou um ajuste ao encontrar seus momentos (essencialmente usandoesse código)

Infelizmente, isso resulta em uma medição incorreta da posição do pico, devido ao sinal de sobreposição dos picos vizinhos. Abaixo está um gráfico da soma dos ajustes separados. Obviamente, todos os seus picos se inclinam em direção ao centro. Preciso levar isso em conta para retornar a posição correta do pico.

Eu tenho um código de trabalho que plota uma função de envelope gaussiano 2D (twoD_Gaussian ()), e analiso isso através de optimize.leastsq como uma matriz 1D usando numpy.ravel e uma função de erro apropriada, no entanto, isso resulta em uma saída sem sentido.

Tentei ajustar um único pico dentro da soma e obtive a seguinte saída incorreta:

Agradeço qualquer conselho sobre o que eu poderia tentar fazer com que isso funcione, ou abordagens alternativas, se isso não for apropriado. Todas as informações são bem-vindas, é claro!

Código abaixo:

from scipy.optimize import leastsq
import numpy as np
import matplotlib.pyplot as plt

def twoD_Gaussian(amp0, x0, y0, amp1=13721, x1=356, y1=247, amp2=14753, x2=291,  y2=339, sigma=40):

    x0 = float(x0)
    y0 = float(y0)
    x1 = float(x1)
    y1 = float(y1)
    x2 = float(x2)
    y2 = float(y2)

    return lambda x, y:  (amp0*np.exp(-(((x0-x)/sigma)**2+((y0-y)/sigma)**2)/2))+(
                             amp1*np.exp(-(((x1-x)/sigma)**2+((y1-y)/sigma)**2)/2))+(
                             amp2*np.exp(-(((x2-x)/sigma)**2+((y2-y)/sigma)**2)/2))

def fitgaussian2D(x, y, data, params):

    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution found by a fit"""
    errorfunction = lambda p: np.ravel(twoD_Gaussian(*p)(*np.indices(np.shape(data))) - data)

    p, success = optimize.leastsq(errorfunction, params)
    return p     

# Create data indices
I = image   # Red channel of a scanned image, equivalent to the  1st image displayed in this post.
p = np.asarray(I).astype('float')
w,h = np.shape(I)
x, y = np.mgrid[0:h, 0:w]
xy = (x,y)

# scanned at 150 dpi = 5.91 dots per mm
dpmm = 5.905511811
plot_width = 40*dpmm

# create function indices
fdims = np.round(plot_width/2)  
xdims = (RC[0] - fdims, RC[0] + fdims)
ydims = (RC[1] - fdims, RC[1] + fdims)
fx = np.linspace(xdims[0], xdims[1], np.round(plot_width))
fy = np.linspace(ydims[0], ydims[1], np.round(plot_width))
fx,fy = np.meshgrid(fx,fy)

#Crop image for display
crp_data = image[xdims[0]:xdims[1], ydims[0]:ydims[1]]
z = crp_data    

# Parameters obtained from separate fits
Amplitudes = (13245, 13721, 15374)
px = (410, 356, 290)
py = (350, 247, 339)

initial_guess_sum = (Amp[0], px[0], py[0], Amp[1], px[1], py[1], Amp[2], px[2], py[2])
initial_guess_peak3 = (Amp[0], px[0], py[0]) # Try fitting single peak within sum


fitted_pars = fitgaussian2D(x, y, z, initial_guess_sum)
#fitted_pars = fitgaussian2D(x, y, z, initial_guess_peak3)

data_fitted= twoD_Gaussian(*fitted_pars)(fx,fy)
#data_fitted= twoD_Gaussian(*initial_guess_sum)(fx,fy)



fig = plt.figure(figsize=(10, 30))
ax = fig.add_subplot(111, aspect="equal")
#fig, ax = plt.subplots(1)
cb = ax.imshow(p, cmap=plt.cm.jet, origin='bottom',
    extent=(x.min(), x.max(), y.min(), y.max()))
ax.contour(fx, fy, data_fitted.reshape(fx.shape[0], fy.shape[1]), 4, colors='w')

ax.set_xlim(np.int(RC[0])-135, np.int(RC[0])+135)
ax.set_ylim(np.int(RC[1])+135, np.int(RC[1])-135)

#plt.colorbar(cb)
plt.show()