Fiting a sum of 2D gaussians to 2d data in python?

Question

I have some 2D data, specifically scanned X-ray films. These have measurements of overlapping point source exposures. Example of data: https://i.sstatic.net/oawlU.png

I want to find the peak positions, by fitting a sum of 2D gaussians to the data. I've tried several other methods with some success, including a "star search" method that locates the global maximum, fits a Gaussian and subtracts it. Looping this finds all of the peaks reasonably well, but it's unstable and not very accurate. I would like to use the star search output as a first guess for a fit, but I'm having trouble, implementing scipy.optimize.curve_fit.

I've made a function twoD_envelope that creates the 2D envelope of all the Gaussians found from the star search. This produces this output: https://i.sstatic.net/4KnpG.png

I was under the impression that I could use this to as an initial guess in curve_fit, however I get the following TypeError:

TypeError: twoD_envelope() takes 4 positional arguments but 358802 were given

358802 is 1 more than size of the data, which is a huge clue, but I can't work out what the problem is! I'm a physicist with "pragmatic" coding knowledge, so any input is very much appreciated.

The code is below.

def twoD_envelope(npoints, xls, yls, pars):

    envl = copy.copy(sq_image)
    envl[:] = 0


    for n in range(0,npoints):
        height, cx, cy, width_x, width_y = pars[n]
        FWHM = 0.5*(width_x+width_y)
        g=makeGaussian(shape(sq_image)[0],fwhm=FWHM,center=[cx+xls[n],cy+yls[n]])
        envl = envl + g

    return envl.ravel()

# Create x and y indices
x = np.linspace(0, np.size(sq_image[0]), np.size(sq_image[0])+1)
y = np.linspace(0, np.size(sq_image[1]), np.size(sq_image[1])+1)
x, y = np.meshgrid(x, y)
coords = [x,y]

data = sq_image

initial_guess = twoD_envelope(9,xls,yls,pars)

pars_opt, pars_cov = opt.curve_fit(twoD_envelope, coords, data, p0=initial_guess)

(sq_image is the data, an ndarray with shape (599,599))

(pars, xls, yxl = lists of Gaussian fit parameters from star search)

(makeGaussian is a function defined elsewhere)

Answer 1

I think your error message is a result of the function you pass to curve_fit which doesn't look like it is of the correct form. This may result in the data array interpreted as vars which results in the TypeError you see (difficult to tell without being able to run the code). According to the documentation, the function passed to curve_fit should,

take the independent variable as the first argument and the parameters to fit as separate remaining arguments.

You may also need to use a star before initial_guess to pass it as a tuple, e.g.

pars_opt, pars_cov = opt.curve_fit(twoD_envelope, coords, data, p0=*initial_guess)

To fit a single 2D Gaussian, where p0 will be about 7 parameters, there is a very good answer which may help: Fitting a 2D Gaussian function using scipy.optimize.curve_fit - ValueError and minpack.error.

A possible solution to your problem may be to loop through your 2D sq_image and use a single 2D Gaussian fit locally with the initial parameters for each from the star search...

EDIT: Code to fit Gaussian.

import scipy.optimize as opt
import numpy as np
import pylab as plt
import matplotlib.cm as cm
import Image


def twoD_Gaussian((x, y), amplitude, xo, yo, sigma_x, sigma_y, theta):
    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 = amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo) 
                            + c*((y-yo)**2)))
    return g.ravel()


# Create x and y indices
I = Image.open('./30155885.png')
p = np.asarray(I).astype('float')
w,h = I.size
x, y = np.mgrid[0:h, 0:w]

#Use only one channel of image
p = p[:,:,0]

#Fit 2D Gaussian
initial_guess = (3,10,10,20,40,0)
popt, pcov = opt.curve_fit(twoD_Gaussian, (x, y), p.ravel(), p0=initial_guess)
data_fitted = twoD_Gaussian((x, y), *popt)


fig, ax = plt.subplots(1, 1)
cb = ax.imshow(p, cmap=plt.cm.jet, origin='bottom',
    extent=(x.min(), x.max(), y.min(), y.max()))
ax.contour(x, y, data_fitted.reshape(x.shape[0], y.shape[1]), 8, colors='w')


plt.colorbar(cb)
plt.show()

where image 30155885 is,

Note that only one channel of the image data is used (you should replace the data array with the one from sq_image). This results in,