File Exchange

image thumbnail

Fit 2D Gaussian with Optimization Toolbox

version 1.2.0.0 (2.55 KB) by Nathan Orloff
Fit a 2D rotated gaussian. http://en.wikipedia.org/wiki/Gaussian_function
5.0
10 Ratings

21 Downloads

Updated 30 May 2013

View Version History

View License

FMGAUSSFIT performs a gaussian fit on 3D data (x,y,z).
[fitresult,..., rr] = fmgaussfit(xx,yy,zz) uses ZZ for the surface
height. XX and YY are vectors or matrices defining the x and y
components of a surface. If XX and YY are vectors, length(XX) = n and
length(YY) = m, where [m,n] = size(Z). In this case, the vertices of the
surface faces are (XX(j), YY(i), ZZ(i,j)) triples. To create XX and YY
matrices for arbitrary domains, use the meshgrid function. FMGAUSSFIT
uses the lsqcurvefit tool, and the OPTIMZATION TOOLBOX. The initial
guess for the gaussian is places at the maxima in the ZZ plane. The fit
is restricted to be in the span of XX and YY.
See:
http://en.wikipedia.org/wiki/Gaussian_function

Examples:
To fit a 2D gaussian:
[fitresult, zfit, fiterr, zerr, resnorm, rr] =
fmgaussfit(xx,yy,zz);
See also SURF, OPTIMSET, LSQCURVEFIT, NLPARCI, NLPREDCI.

Cite As

Nathan Orloff (2021). Fit 2D Gaussian with Optimization Toolbox (https://www.mathworks.com/matlabcentral/fileexchange/41938-fit-2d-gaussian-with-optimization-toolbox), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (16)

Thomas Nowack

Thanks, it's a really useful!
I was just wondering what will be returned as fitresult(1) % amplitude in case of a z-offset (e.g. image with darkcount > 0 for the area surrounding the Gaussian peak)? Will it return the difference between peak-maximum and 0 or the difference between peak-maximum and the darkcount/level (> 0) of surrounding intensity?
I tried to find it in the code, but I'm still a beginner in Matlab...
Any help is greatly appreciated!

Shira Taragin

Luis Oviedo

Kristoffer Kjærnes

Sergey

Could you please tell me, how should I modify the code to fit data with symmetric gaussian, i.e. i want sigmas to be equal ; I am not sure I matched code in gaussian2D with formulas correctly.

Pavlo Kliuiev

thank you!

Juan Garcia

Hi Nathan,

You're right it is a constant, but if you want to actually retrieve the sigmas of the gaussian you need that constant.

Try the simple case
x = -100:100; y = x; [xx, yy] = meshgrid(x,y);
sigma = 20;
zz = exp(-(xx.^2 + yy.^2) ./ (2*sigma^2));

You obtain fitresult(3:4) = 28.2847 when you would like it equal to 20. It can be corrected afterwards by dividing fitresult(3:4)./sqrt(2), but just to be aware of it.

Regards

Nathan Orloff

Hi Juan,

I think that Zo is just an offset. As to the normalization of the sigma, I think it is just a constant. I omitted it for readability.

Thanks GBorg!

GBorghesan

The fitresults is

mu=fitresult(5:6)
angle=fitresult(2)%in deg
sigma= fitresult(3:4)
hg=fitresult(1) % amplitude ;

and can be reused in function
gaussian2D=@(par,x,y)(par(7) + ...
par(1)*exp(-(((x-par(5)).*cosd(par(2))+(y-par(6)).*sind(par(2)))./par(3)).^2-...
((-(x-par(5)).*sind(par(2))+(y-par(6)).*cosd(par(2)))./par(4)).^2));

XYZ

Juan Garcia

Nice function. I was having a look at it and I have a question and a comment.

Question: What is the meaning of z0 or par(7) in the gaussian2D function.

Comment: I would say that both sigmas in the gaussian2D function lack the multiplication by 2, or to maintain the power of 2 to the whole fraction, by sqrt(2).

Then:
function z = gaussian2D(par,xy)
% compute 2D gaussian
z = par(7) + ...
par(1) * exp(-(((xy{1}-par(5)).*cosd(par(2))+(xy{2}-par(6)).*sind(par(2)))./(sqrt(2)*par(3))).^2-...
((-(xy{1}-par(5)).*sind(par(2))+(xy{2}-par(6)).*cosd(par(2)))./(sqrt(2)*par(4))).^2);
end

Eric T

I get the following warning:
Warning: The Levenberg-Marquardt algorithm does not handle bound constraints; using the
trust-region-reflective algorithm instead.
> In lsqncommon at 83
In lsqcurvefit at 252
In fmgaussfit at 56

Nathan Orloff

fitresult = fit parameters used to generate the fit.

zfit = the z values of the fit

fiterr = error in the fitparameters assuming a confidence interval.

zerr = error in the z values assuming the uncertainty in the fit paramteres

resnorm = residual

rr = reduced chi-squared.

Olhado

This bit of code is doing exactly what I want I think but I wondering if it might be possible for the Author to expand on the meaning of the values the function outputs or point me in the direction as to where I can find the answer. Thanks a Ton

Christopher

Nathan Orloff

I edited this a little. The fit is better with z bounds. I will post it soon... If you have already downloaded this just add the following.

%% Set up the startpoint
[amp, ind] = max(zData); % amp is the amplitude.
xo = xData(ind); % guess that it is at the maximum
yo = yData(ind); % guess that it is at the maximum
ang = 45; % angle in degrees.
sy = 1;
sx = 1;
zo = median(zData(:))-std(zData(:));
xmax = max(xData)+2;
ymax = max(yData)+2;
zmax = amp*2; % amp is the amplitude.
xmin = min(xData)-2;
ymin = min(yData)-2;
zmin = min(zData)/2; % amp is the amplitude.

%% Set up fittype and options.
Lower = [0, eps, 0, 0, xmin, ymin, zmin];
Upper = [Inf, 180, Inf, Inf, xmax, ymax, zmax]; % angles greater than 90 are redundant
StartPoint = [amp, ang, sx, sy, xo, yo, zo];%[amp, sx, sxy, sy, xo, yo, zo];

MATLAB Release Compatibility
Created with R2012b
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!