Fitting a 2D Gaussian function using scipy.optimize.curve_fit - ValueError and minpack.error

The output of twoD_Gaussian needs to be 1D. What you can do is add a .ravel() onto the end of the last line, like this:

def twoD_Gaussian((x, y), amplitude, xo, yo, sigma_x, sigma_y, theta, offset):    xo = float(xo)    yo = float(yo)        a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)    g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)                             + c*((y-yo)**2)))    return g.ravel()

You'll obviously need to reshape the output for plotting, e.g:

# Create x and y indicesx = np.linspace(0, 200, 201)y = np.linspace(0, 200, 201)x, y = np.meshgrid(x, y)#create datadata = twoD_Gaussian((x, y), 3, 100, 100, 20, 40, 0, 10)# plot twoD_Gaussian data generated aboveplt.figure()plt.imshow(data.reshape(201, 201))plt.colorbar()

Do the fitting as before:

# add some noise to the data and try to fit the data generated beforehandinitial_guess = (3,100,100,20,40,0,10)data_noisy = data + 0.2*np.random.normal(size=data.shape)popt, pcov = opt.curve_fit(twoD_Gaussian, (x, y), data_noisy, p0=initial_guess)

And plot the results:

data_fitted = twoD_Gaussian((x, y), *popt)fig, ax = plt.subplots(1, 1)ax.hold(True)ax.imshow(data_noisy.reshape(201, 201), cmap=plt.cm.jet, origin='bottom',    extent=(x.min(), x.max(), y.min(), y.max()))ax.contour(x, y, data_fitted.reshape(201, 201), 8, colors='w')plt.show()

To expand on Dietrich's answer a bit, I got the following error when running the suggested solution with Python 3.4 (on Ubuntu 14.04):

def twoD_Gaussian((x, y), amplitude, xo, yo, sigma_x, sigma_y, theta, offset):                  ^SyntaxError: invalid syntax

Running 2to3 suggested the following simple fix:

def twoD_Gaussian(xdata_tuple, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):    (x, y) = xdata_tuple                                                            xo = float(xo)                                                                  yo = float(yo)                                                                  a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)       b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)        c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)       g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)                                 + c*((y-yo)**2)))                                       return g.ravel()

The reason for this is that automatic tuple unpacking when it is passed to a function as a parameter has been removed as of Python 3. For more information see here: PEP 3113

curve_fit() wants to the dimension of xdata to be (2,n*m) and not (2,n,m). ydata should have shape (n*m) not (n,m) respectively. So you use ravel() to flatten your 2D arrays:

xdata = np.vstack((xx.ravel(),yy.ravel()))ydata = data_noisy.ravel()popt, pcov = opt.curve_fit(twoD_Gaussian, xdata, ydata, p0=initial_guess)

By the way: I'm not sure if the parametrization with the trigonometric terms is the best one. E.g., taking the one described here might be a bit more robust under numerical aspects and large deviations.