Matlab code for 2D Gaussian surface fitting

Updated 10/21/2011

I have some code on Matlab Central to automatically fit a 1D Gaussian to a curve and a 2D Gaussian or Gabor to a surface. The 2D Gaussian code can optionally  fit a tilted Gaussian.

In its basic form curve/surface fitting is straightforward (a call to lsqcurvefit will do the trick), but the problem is that the error function to be minimized has no guarantees of convexity. This brings two issues:

  1. The minimization can get stuck in local minima
  2. Error bars based on the second order derivatives of the error function are meaningless

I solved the first issue by an initial exhaustive search through the fit parameters.

For the second issue, I discussed earlier how to track receptive field envelopes through time using a hierarchical Bayesian model and WinBUGS. The same basic idea can be used to fit static curves or surfaces. Simply cast the estimation problem as a Bayesian inference problem. For the 2D Gaussian for example:

zi ~ N(wi,sigma)
wi <-  a*exp(-((xi-x0).^2/2/sigmax^2 + (yi-y0).^2/2/sigmay^2)) + b
x0 ~ Uniform(min(xi),max(xi))
y0 ~ Uniform(min(yi),max(yi))
sigmax, sigmay, a, b ~ N(0,1000)

Because this model is analytically intractable, inference is done through MCMC sampling. I’ve used the dram (delayed rejection adaptive Metropolis) flavor of MCMC together with CODA functions available in the Spatial Econometrics toolbox to monitor convergence. Thus, the MCMC aspect is fully automatic, and it seems to be very reliable, although as they say in the WinBUGS documentation, MCMC is dangerous, use with caution.

2 responses to “Matlab code for 2D Gaussian surface fitting”

  1. Hi Patrick,
    This routine is a very useful program and seems work great for me. But I have some questions for my own purpose and maybe you can give me some suggestions.

    If I want to use this function to monitor a fixed point source over time (take 1000 images of the same bright spot over 1min). The intensity of the image will change over time but bright spot always shows up in the first frame. So I want to use the program to get the center position of the Gaussian surface (x0 and y0) from the first frame and share the center position for the rest 999 frames and see how the amplitude and sigma changes as a function of time. How can I do it ? Thanks for the help.

Leave a comment

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s