COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

3d Gauss fitting

• Email
• Follow

```Hi ppl
I have a 100 x 100 data which is in the form of a gaussian distribution.The data has too much noise and I want to fit a 3d gaussian to the data. I have the surface fitting tool box but it says that x,y,z dimensions should be same. But my data is like x=1x100 y=1x100 and z is 100x100. Any help would be highly appreciated.
```
 0

See related articles to this posting

```Try the following

[X,Y] = meshgrid(x,y)

X = X(:)
Y = Y(:)
Z = z(:)

"Aishwarya " <icemails@yahoo.co.in> wrote in message
news:hdv3dg\$8in\$1@fred.mathworks.com...
> Hi ppl
> I have a 100 x 100 data which is in the form of a gaussian
> distribution.The data has too much noise and I want to fit a 3d gaussian
> to the data. I have the surface fitting tool box but it says that x,y,z
> dimensions should be same. But my data is like x=1x100 y=1x100 and z is
> 100x100. Any help would be highly appreciated.

```
 0

```On Nov 17, 4:09=A0pm, "Aishwarya " <icema...@yahoo.co.in> wrote:
> Hi ppl
> I have a 100 x 100 data which is in the form of a gaussian distribution.T=
he data has too much noise and I want to fit a 3d gaussian to the data. I h=
ave the surface fitting tool box but it says that x,y,z dimensions should b=
e same. But my data is like x=3D1x100 y=3D1x100 and z is 100x100. Any help =
would be highly appreciated.

--------------------------------
How is this 3D data?  You have two independent variables, x and y, and
a value - this makes it a 2D function.  Repeat: a 2D function.

What is your z variable?  Some kind of image?  It's a 2D array of
100x100 but what is it?
```
 0

```Ya.It is a 2d function as in z is a function of x and y.

I did not get that question reg z
```
 0

```On Nov 18, 1:48=A0am, "Aishwarya " <icema...@yahoo.co.in> wrote:
> Ya.It is a 2d function as in z is a function of x and y.
>
> I did not get that question reg z

--------------------------------------------------
Maybe there's a built in function but you can do it the regular,
normal manual way where
xMean =3D sum(x *z) / sum(z)
yMean =3D sum(y * z) / sum(z)
and the standard deviations are the usual formulas
xStd =3D sqrt( mean(x^2 - xMean))
etc.
Of course you have to scan every pixel in the array to calculate these
values.
```
 0

```Hi,
I tried the following code in MATLAB 2009b :
clc;
clear all;
[X,Y] = meshgrid(1:100,1:100);
% x=1:1:100;
% y=1:1:100;
x=X(:);
y=Y(:);
xdata = {x,y};
fid = fopen('S3D6.bin','r+');
size(A);
%figure,imagesc(A)
fid = fopen('D6S3.bin','r+');
size(B);
%figure,imagesc(B)
C=A.*B;
size(C);
D= reshape(C,100,100,100);
D1=D(:,:,10);
D2=D1(:);
fun = @(c,xdata) c(1)*exp(-((c(2)*(xdata{1}-30)^2)+(2*c(3)*(xdata{1}-30)*(xdata{2}-40))+(c(4)*(xdata{2}-40)^2)));
c_start=[30 50 30 50];
options=optimset('TolFun',1e-12,'TolX',1e-12,'MaxFunEvals',40000,'MaxIter',50000);
t = lsqcurvefit(fun,c_start,xdata,D1,options);
% [INLP,ILP] = fminspleas(funlist,NLPstart,xdata,D1)
mesh(D1);
op= t(1)*exp(-((t(2)*(xdata{1}-30)^2)+(2*t(3)*(xdata{1}-30)*(xdata{2}-40))+(t(4)*(xdata{2}-40)^2)));
figure,mesh(op);

and im getting the error as: lsqcurvefit accepts input of type double when i include options.
When I dont include options it takes the initial params as the fitting params since the tolerance is 1e-6 and my data is e-20 range.

Any idea???
```
 0

```"Aishwarya " <icemails@yahoo.co.in> wrote in message
news:he2osv\$a73\$1@fred.mathworks.com...
> Hi,
> I tried the following code in MATLAB 2009b :

*snip*

>     t = lsqcurvefit(fun,c_start,xdata,D1,options);

*snip*

> and im getting the error as: lsqcurvefit accepts input of type double when
> i include options.
> When I dont include options it takes the initial params as the fitting
> params since the tolerance is 1e-6 and my data is e-20 range.
>
> Any idea???

Yes -- you're calling LSQCURVEFIT incorrectly.  If you don't want to specify
some of the inputs, you can't just leave them out.  Use empty matrices for
the bound inputs instead of omitting them.  Reread HELP LSQCURVEFIT for a
description of the correct calling syntax.

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ

```
 0

```Hi,

I had the same problem. I saw a bunch of solutions out there. I wrote my own solution with fmincon first. And then I realized that lsqcurvefit could be tricked into working.

Here is my solution.

Best of luck,

Nate

function [fitresult, zfit, fiterr, zerr, resnorm] = fmgaussfit(xx,yy,zz)
%fmgaussfit(xx,yy,zz) Custom Guassian Fit algorithm using lsqcurvefit
%   A 3D gaussian fitting routine with error propagation and uncertainties.

%% Condition the data
[xData, yData, zData] = prepareSurfaceData( xx, yy, zz );
xyData = {xData,yData};

%% 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;
xmin = min(xData)-2;
ymin = min(yData)-2;

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

tols = 1e-16;
options = optimset('Algorithm','levenberg-marquardt',...
'Display','off',...
'MaxFunEvals',5e2,...
'MaxIter',5e2,...
'TolX',tols,...
'TolFun',tols,...
'TolCon',tols ,...
'UseParallel','always');

%% perform the fitting
[fitresult,resnorm,residual] = ...
lsqcurvefit(@gaussian2D,StartPoint,xyData,zData,Lower,Upper,options);
[fiterr, zfit, zerr] = gaussian2Duncert(fitresult,residual,jacobian,xyData);

end

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

function [dpar,zf,dzf] = gaussian2Duncert(par,resid,xy)
J = guassian2DJacobian(par,xy);
parci = nlparci(par,resid,'Jacobian',J);
dpar = (diff(parci,[],2)./2)';
[zf,dzf] = nlpredci(@gaussian2D,xy,par,resid,'Jacobian',J);
end

function J = guassian2DJacobian(par,xy)
x = xy{1}; y = xy{2};
J(:,1) = exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2);
J(:,2) = -par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(3).^2 - (2.*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2);
J(:,3) = (2.*par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2)./par(3)^3;
J(:,4) = (2.*par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2)./par(4)^3;
J(:,5) = par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*cosd(par(2)).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))))./par(3).^2 - (2.*sind(par(2)).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2);
J(:,6) = par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*cosd(par(2)).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2 + (2.*sind(par(2)).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))))./par(3).^2);
J(:,7) = ones(size(x));
end

function [sse, zf] = gaussian2Derr(par,xy,z)
zf = gaussian2D(par,xy);
zerr = zf-z;
sse = sum(zerr.^ 2);
end
```
 0
Reply nathan.orloff (55) 5/23/2013 4:59:13 PM

7 Replies
500 Views

Similar Articles

12/6/2013 6:40:03 AM
page loaded in 4042 ms. (0)