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 leastsquares (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(((xix0).^2/2/sigmax^2 + (yiy0).^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 = (xix0)*cos(theta) + (yii0)*sin(theta);
yip =(xix0)*sin(theta) + (yii0)*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.
1.6.0.0  Bug fixes in autoGaussianCurve, doFinalOptimization 

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

1.4.0.0  Added MetropolisHastings method of estimating model params 

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

1.2.0.0  Added mex version of Gibbs sampler and basic convergence diagnostics for MCMC 
Create scripts with code, output, and formatted text in a single executable document.
Yanjun Han (view profile)
I have tried this example:
[Xi,Yi] = meshgrid(5:0.1:5,5:.1:5);
Zi = exp(((Xi0).^2+(Yi0).^2));
figure(1)
mesh(Xi,Yi,Zi);
axis equal
results = autoGaussianSurf(Xi,Yi,Zi);
% save results;
%load results;
a=results.a;
b=results.b;
x0=results.x0;
y0=results.y0;
sigmax=results.sigmax;
sigmay=results.sigmay;
[xi,yi] = meshgrid(5:0.01:5,5:0.01:5);
%x0=0; y0=0;
zi = a*exp(((xix0).^2/2/sigmax^2 + (yiy0).^2/2/sigmay^2)) + b;
hold on
plothandle= mesh(xi,yi,zi);
axis equal;
Now ,I found the fit result is not agree with my data, can everyone help to explain this phenomenon?
weikun (view profile)
Hi~~very glad to have these spendid codes. But how to use this for a 1D Gabor fit?
Christian GonzalezCapizzi (view profile)
This keeps fitting a Gaussian to my Gabor like plot for some reason. I can't figure out why, can anyone help?
Christian GonzalezCapizzi (view profile)
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
Jim Morgenstern (view profile)
thanks for providing.
I did add a line:
p.addOptional('Display', 'none');
so i could control the warning outputs from my code.
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.
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!
md. (view profile)
i want to use sum of two Gaussian i.e [zi=a0*exp(((xix0).^2/2/sigmax^2 +a1*exp(((xix1).^2/2/sigmax1^2 + b] in the autoGaussianCurve(xi, yi) function.
Here i need to change. Please help me.
Ivan Mikhaylov (view profile)
I also have the same error as described by Leo. How to fix it?
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(((xix0).^2/2/sigmax^2 + (yiy0).^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.
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 =
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?
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.
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?
Rakesh Chalasani (view profile)
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.
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.
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?
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.
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.
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...
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 = (xi4)*cos(pi/6) + yi*sin(pi/6);
yip =(xi4)*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/104)*cos(pi/6) + yi/10*sin(pi/6);
yip =(xi/104)*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.
Jan Churan (view profile)
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.
Matt Munson (view profile)
How would one go about adapting this to fit N gaussians?
Joseph Shtok (view profile)
Actually, this is the way you state in the definition.
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.
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?
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).
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.0130e05
x0: 1.5000
y0: 1.5000
sigmax: 0.5838
sigmay: 0.5838
sse: 1.1458e14
G: [3x3 double]
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
x0,y0=1; Your script.
But even by eyes you can see that x0 and y0 = 1.5