Auto Gaussian & Gabor fits

Fit 1D/2D Gaussian or Gabor to a curve/surface without start guesses for params

Description

Auto Gaussian & Gabor Surface fit
---
Functions to fit a 1D Gaussian to a curve and a 2D Gaussian or Gabor to a surface. The routines are automatic in the sense that they do not require the specification of starting guesses for the model parameters. This is done by evaluating the quality of fit for many different choices of parameters then refining the most promising set of params through least-squares (exhaustive search followed by refinement).

All functions support 2 methods for computing error bars on the parameters: bootstrapping and MCMC.

autoGaussianSurf(xi,yi,zi) fits a 2D Gaussian to a surface, defined as:

zi = a*exp(-((xi-x0).^2/2/sigmax^2 + (yi-y0).^2/2/sigmay^2)) + b

It can also fit a tilted 2d Gaussian or isotropic 2d Gaussian.

autoGaborSurf(xi,yi,zi) fits a Gabor, defined as:

zi = a*exp(-(xip,.^2+yip.^2)/2/sigma^2)*cos(2*pi*xip/lambda + phase) + b

Where:
xip = (xi-x0)*cos(theta) + (yi-i0)*sin(theta);
yip =-(xi-x0)*sin(theta) + (yi-i0)*cos(theta);

The Gabor fit calls autoGaussianSurf internally, using the fact that the absolute value of a Gabor in the Fourier domain is a Gaussian.

autoGaussianCurve(xi,zi) fits a 1D Gaussian to a curve.

• MATLAB release: MATLAB 7.9 (2009b)

[bmv]

Something wrong with the definition of x0 and y0
A simple example:
    0.0247 0.1074 0.0247
    0.1074 0.4661 0.1074
    0.0247 0.1074 0.0247
x0,y0=1; Your script.
 But even by eyes you can see that x0 and y0 = 1.5

[Patrick Mineault]

I think you are misinterpreting the results. The coordinate system that the method depends on what you put as xi and yi. So if your convention is that the first pixel at the top is (0.5, 0.5), this code works correctly:

r = [0.0247 0.1074 0.0247
     0.1074 0.4661 0.1074
     0.0247 0.1074 0.0247];
[xi,yi] = meshgrid(0.5:2.5,0.5:2.5);
results = autoGaussianSurfML(xi,yi,r)

results =

         a: 0.4662
         b: -8.0130e-05
        x0: 1.5000
        y0: 1.5000
    sigmax: 0.5838
    sigmay: 0.5838
       sse: 1.1458e-14
         G: [3x3 double]

[Joseph Shtok]

The code works great, up to one small correction: in the autoGaussianSurfML.m you use the command
>>[~,pos] = max(x);
to find the position of the maximal value in x. My 2008 Matlab did not eat this, so I replace the command with
>>[maxval,pos] = max(x);
(the maxlval is not used later).

[Federico Pinchetti]

Matlab 2007b:

??? Undefined function or method 'range' for input arguments of type 'double'.

Error in ==> autoGaussianSurfML at 35
    rg = (range(boundx)+range(boundy))/2;

how do I solve?

[Patrick Mineault]

range is in the stats toolbox. If you don't have it you can write it yourself:

range = @(x) max(x(:))-min(x(:));

Only 2009b and onwards are supported. It won't work in 2007b or 2008 without some changes.

[Joseph Shtok]

Actually, this is the way you state in the definition.

[Matt Munson]

How would one go about adapting this to fit N gaussians?

[Patrick Mineault]

A standard approach to fitting a curve or surface (including 2 gaussians) is to minimize the squared error between the function and the measured data. You can use lsqcurvefit for this purpose.
One complication is that the squared error might have multiple local minima, and you can easily get stuck with suboptimal parameters. You can avoid this by specifying start parameters close to the true minimum (which is annoying) or by performing the optimization several times with different start parameters and selecting the best fit out of the lot (which is expensive). For a single Gaussian it turns out that there is a third, smarter way of avoiding local minima; but this only works for a single Gaussian. So this is not the best script to start with to fit 2 gaussians to a surface.

[Jan Churan]

 

[Qiuyan PENG]

It seems that if the step size of xi and yi is not integer, there are some problems.

For example, if I run the following code, the fitting plot is quite bad.

[xi,yi] = meshgrid(-10:0.1:10,-20:0.1:20);
xip = (xi-4)*cos(pi/6) + yi*sin(pi/6);
yip =-(xi-4)*sin(pi/6) + yi*cos(pi/6);
zi = exp(-(xip.^2+yip.^2)/2/2^2).*cos(xip*2*pi/5+pi/3) + .2*randn(size(xi));
results = autoGaborSurf(xi,yi,zi);
subplot(1,2,1);imagesc(xi(:),yi(:),zi);
subplot(1,2,2);imagesc(xi(:),yi(:),results.G);

However, if I replace the first thee lines with the following, it works perfectly.

[xi,yi] = meshgrid(-100:100,-200:200);
xip = (xi/10-4)*cos(pi/6) + yi/10*sin(pi/6);
yip =-(xi/10-4)*sin(pi/6) + yi/10*cos(pi/6);

However, they are basically the same thing for zi, except for the xi and yi step size.

I guess this is due to the transformation into frequency domain when you do not consider the sampling effect?

Thanks.

[Qiuyan PENG]

Also, is it the case that the "autoGaussianSurf" only supports gaussian with diagonal covariance matrix? I know that "autoGaborSurf" supports orientation, but it would always assume some sinusoidal modulation...

[Patrick Mineault]

Quiyan: you're right, that's a bug. The issue is in line 67: it should go lambda = 1/sf*da and not lambda = 1/sf/da.

Yes, Gaussian curve fitting is currently restricted to diagonal covariance matrix, though perhaps I'll add tilt to the next version.

[Qiuyan PENG]

Thanks, Patrick. The bug is fixed as you suggested.

I've finished implementing the tilted version (also with some aspect ratio) for Gaussian by modifying you Gabor fitting functions. Thank you.

[Ben]

Hi Patrick,

Thanks for the code.

From my tests, the code always fit a Gaussian surface for any code given to it. Thus the question comes: how to determine whether the fit is acceptable or not?

[Patrick Mineault]

Ben, look at the r2 values that Matlab gives you. If it's 1, then it's a perfect fit, if it's 0, then not. The other cue is if you compute error bars using one of the many methods that's included in the code, and your error bars are huge, that's a signal of misfit.

[Ben]

Patric, is the error bars you mentioned "sse" or "sse0"? Besides, using bootstrap, the program gave out negative r2 while, for the save data, it gave positive r2 using mcmc. Why negative r2? Which "error bar" is more trustable? Thanks.

[Rakesh Chalasani]

 

[Daniel]

Hi Patric, very nice program, however I think it is not acting "robustly" for me. It keeps conking out after about 10 iterations for each parameter, saying a local minimum is found, and saying that the iteration has yielded a tolerance less than my TolX... however I can't seem to create my own options structure with tolerances that this will use. The upshot is my R2 is negative and my errors are huge! Any thoughts?

[Patrick Mineault]

Ben, bootstrapping might work incorrectly if the Gaussian bump takes only a few pixels. In general, MCMC should be more reliable. The only reason I would use bootstrapping over MCMC is that MCMC is mathematically complicated and it's a bit tough to explain in the methods section of a paper.

Daniel, it's difficult to say without having access to the data you have. It's possible that there's a bug in the code that only appears for very specific datasets. Try to boil it down to the simplest example that shows the bug and send me the code to replicate the problem.