Code covered by the BSD License

### Highlights from Auto Gaussian & Gabor fits

4.5
4.5 | 10 ratings Rate this file 13 Downloads (last 30 days) File Size: 12.4 KB File ID: #31485 Version: 1.6

# Auto Gaussian & Gabor fits

### Patrick Mineault (view profile)

18 May 2011 (Updated )

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

File Information
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 (R2009b)
17 Jun 2016 Christian Gonzalez-Capizzi

### Christian Gonzalez-Capizzi (view profile)

This keeps fitting a Gaussian to my Gabor like plot for some reason. I can't figure out why, can anyone help?

Comment only
08 Jun 2016 Christian Gonzalez-Capizzi

### Christian Gonzalez-Capizzi (view profile)

02 Mar 2015 Jim Morgenstern

### Jim Morgenstern (view profile)

i would like to use a mask to include only certain pixels in the calculations; does lsqcurvefit already allow for a mask? how would i implement such a mask here ? thanks

Comment only
02 Mar 2015 Jim Morgenstern

### Jim Morgenstern (view profile)

thanks for providing.
so i could control the warning outputs from my code.

11 Feb 2015 Patrick Mineault

### Patrick Mineault (view profile)

Tarun: the new URL for the DRAM code is http://helios.fmi.fi/~lainema/dram/dramcode.zip

Others: there's been some recent changes to the optimization toolbox, and it really depends on the version you're using. State the version of Matlab you're having trouble with and I might be able to help.

Comment only
08 Jan 2015 Tarun Narayan

### Tarun Narayan (view profile)

I was trying to use the mcmc error bar estimation method, but it appears that http://www.helsinki.fi/~mjlaine/dram/dramcode.zip is no longer accessible. Do you know what an appropriate fix is?
Thanks!

11 Oct 2014 md.

### md. (view profile)

i want to use sum of two Gaussian i.e [zi=a0*exp(-((xi-x0).^2/2/sigmax^2 +a1*exp(-((xi-x1).^2/2/sigmax1^2 + b] in the autoGaussianCurve(xi, yi) function.

19 Feb 2014 Ivan Mikhaylov

### Ivan Mikhaylov (view profile)

I also have the same error as described by Leo. How to fix it?

Comment only
18 Apr 2013 yuanmh

### yuanmh (view profile)

I have used the 'autoGaussianSurf' to estimate 2d Gaussian distribution,the result turns out to be not very good, the surface function is something like
zi = a*exp(-((xi-x0).^2/2/sigmax^2 + (yi-y0).^2/2/sigmay^2)) + b,but in my function the 'a' is equal to 1,and 'b' is zero,but the parameters 'a' and 'b' estimate from 'autoGaussianSurf' is 0.7 and -0.5.My question is if i have know the value of parameters 'a' and 'b' ,how can i change the program to make the 'autoGaussianSurf' fitting the surface better.

27 Jun 2012 Leo

### Leo (view profile)

I think this submission requires the Optimization toolbox?

I don't have it, and when I try to run autoGaussianCurve, I get the following error:

Error using optimset (line 203)
Unrecognized parameter name 'Jacobian'. Please see the
optimset reference page in the documentation for a list of
acceptable option parameters. Link to reference page.

Error in doFinalOptimization>getMLestimate (line 101)
opts = optimset('Display',display,'Jacobian','on');

Error in doFinalOptimization (line 8)
params =

Basically it seems that setting this option "Jacobian" is only valid if you have the Optimization toolbox... but I might be wrong. Any hints?

Comment only
26 Mar 2012 Patrick Mineault

### Patrick Mineault (view profile)

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.

Comment only
20 Mar 2012 Daniel

### Daniel (view profile)

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?

29 Feb 2012 Rakesh Chalasani

21 Feb 2012 Ben

### Ben (view profile)

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.

Comment only
20 Feb 2012 Patrick Mineault

### Patrick Mineault (view profile)

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.

Comment only
20 Feb 2012 Ben

### Ben (view profile)

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?

Comment only
17 Oct 2011 Qiuyan PENG

### Qiuyan PENG (view profile)

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.

17 Oct 2011 Patrick Mineault

### Patrick Mineault (view profile)

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.

Comment only
16 Oct 2011 Qiuyan PENG

### Qiuyan PENG (view profile)

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...

Comment only
16 Oct 2011 Qiuyan PENG

### Qiuyan PENG (view profile)

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.

Comment only
06 Sep 2011 Jan Churan

### Jan Churan (view profile)

02 Sep 2011 Patrick Mineault

### Patrick Mineault (view profile)

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.

Comment only
31 Aug 2011 Matt Munson

### Matt Munson (view profile)

Comment only
23 Aug 2011 Joseph Shtok

### Joseph Shtok (view profile)

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

Comment only
29 Jul 2011 Patrick Mineault

### Patrick Mineault (view profile)

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.

Comment only
29 Jul 2011 Federico Pinchetti

### Federico Pinchetti (view profile)

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?

Comment only
20 Jul 2011 Joseph Shtok

### Joseph Shtok (view profile)

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).

06 Jun 2011 Patrick Mineault

### Patrick Mineault (view profile)

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]

Comment only
06 Jun 2011 bmv

### bmv (view profile)

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
But even by eyes you can see that x0 and y0 = 1.5

Comment only
08 Jun 2011 1.2

Added mex version of Gibbs sampler and basic convergence diagnostics for MCMC

01 Sep 2011 1.3

Removed Gibbs sampling version of Gaussian surface fit (was unreliable). Added Gabor fitting. Changed function names. Uses better limits for sigmax, sigmay for Gaussian fit.

17 Sep 2011 1.4

Added Metropolis-Hastings method of estimating model params

21 Oct 2011 1.5

Added Gaussian curve fitting. Added MCMC and bootstrap based error bars for every function. Added support for tilted Gaussian

02 Nov 2011 1.6

Bug fixes in autoGaussianCurve, doFinalOptimization