Code covered by the BSD License  

Highlights from
Auto Gaussian & Gabor fits

4.5

4.5 | 6 ratings Rate this file 57 Downloads (last 30 days) File Size: 12.4 KB File ID: #31485
image thumbnail

Auto Gaussian & Gabor fits

by

 

18 May 2011 (Updated )

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

| Watch this File

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)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (22)
19 Feb 2014 Ivan Mikhaylov

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

18 Apr 2013 yuanmh

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

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 =
getMLestimate(curvefun,'Iter',p0,X,z,lb,ub,mask);

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

26 Mar 2012 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.

20 Mar 2012 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?

29 Feb 2012 Rakesh Chalasani  
21 Feb 2012 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.

20 Feb 2012 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.

20 Feb 2012 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?

17 Oct 2011 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.

17 Oct 2011 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.

16 Oct 2011 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...

16 Oct 2011 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.

06 Sep 2011 Jan Churan  
02 Sep 2011 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.

31 Aug 2011 Matt Munson

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

23 Aug 2011 Joseph Shtok

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

29 Jul 2011 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.

29 Jul 2011 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?

20 Jul 2011 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).

06 Jun 2011 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]

06 Jun 2011 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

Updates
08 Jun 2011

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

01 Sep 2011

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

Added Metropolis-Hastings method of estimating model params

21 Oct 2011

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

02 Nov 2011

Bug fixes in autoGaussianCurve, doFinalOptimization

Contact us