Python lmfit: Fitting a 2D Model

Question

I'm trying to fit a 2D-Gaussian to some greyscale image data, which is given by one 2D array. The lmfit library implements a easy-to-use Model class, that should be capable of doing this. Unfortunately the documentation (http://lmfit.github.io/lmfit-py/model.html) does only provide examples for 1D fitting. For my case I simply construct the lmfit Model with 2 independent variables.

The following code seems valid for me, but causes scipy to throw a "minpack.error: Result from function call is not a proper array of floats."

Tom sum it up: How to input 2D (x1,x2)->(y) data to a Model of lmfit.?

Here is my approach: Everything is packed in a GaussianFit2D class, but here are the important parts: That's the Gaussian function. The documentation says about user defined functions

Of course, the model function will have to return an array that will be the same size as the data being modeled. Generally this is handled by also specifying one or more independent variables.

I don't really get what this should mean, since for given values x1,x2 the only reasonable result is a scalar value.

def _function(self, x1, x2, amp, wid, cen1, cen2):
    val = (amp/(np.sqrt(2*np.pi)*wid)) * np.exp(-((x1-cen1)**2+(x2-cen2)**2)/(2*wid**2))
    return val

Here the model is generated:

def _buildModel(self, **kwargs):
    model = lmfit.Model(self._function, independent_vars=["x1", "x2"],
                        param_names=["amp", "wid", "cen1", "cen2"])
    return model

That's the function that takes the data, builds the model and params and calls lmfit fit():

def fit(self, data, freeX, **kwargs):
    freeX = np.asarray(freeX, float)
    model = self._buildModel(**kwargs)
    params = self._generateModelParams(model, **kwargs)

    model.fit(data, x1=freeX[0], x2=freeX[1], params=params)

Anf finally here this fit function gets called:

    data = np.asarray(img, float)
    gaussFit = GaussianFit2D()
    x1 = np.arange(len(img[0, :]))
    x2 = np.arange(len(img[:, 0]))
    fit = gaussFit.fit(data, [x1, x2])

Do you require lmfit, or are other tools (curve_fit or leastsq from scipy) also fine? — user707650
lmfit is charming, because it's easy to add init values for parameters, constraints and so on... If I don't manage do get it working, scipy is probably the next thing to try. — m0e

m0e m0e · Accepted Answer · 2014-10-24T01:41:41

Ok, wrote with the devs and got the answer from them (thanks to Matt here).

The basic idea is to flatten all the input to 1D data, hiding from lmfit the >1 dimensional input. Here's how you do it. Modify your function:

 def function(self, x1, x2):
       return (x1+x2).flatten()

Flatten your 2D input array you want to fit to:

...
data = data.flatten()
...

Modify the two 1D x-variables such that you have any combination of them:

...
x1n = []
x2n = []
    for i in x1:
         for j in x2:
              x1n.append(i)
              x2n.append(j)
x1n = np.asarray(x1n)
x2n = np.asarray(x2n)
...

And throw anything into the fitter:

model.fit(data, x1=x1n, x2=x2n, params=params)

Python lmfit: Fitting a 2D Model

2 Answers