3d Gauss fitting

[Aishwarya]

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.

[Richard Willey]

Try the following

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

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

Post by Aishwarya
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.

[Aishwarya]

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+');
A = fread(fid,'double');
size(A);
%figure,imagesc(A)
fid = fopen('D6S3.bin','r+');
B = fread(fid,'double');
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???

[Steven Lord]

Post by Aishwarya
Hi,

*snip*

Post by Aishwarya
t = lsqcurvefit(fun,c_start,xdata,D1,options);

*snip*

Post by Aishwarya
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.

[ImageAnalyst]

Post by Aishwarya
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.

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

[Aishwarya]

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

I did not get that question reg z

[ImageAnalyst]

Post by Aishwarya
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 = sum(x *z) / sum(z)
yMean = sum(y * z) / sum(z)
and the standard deviations are the usual formulas
xStd = sqrt( mean(x^2 - xMean))
etc.
Of course you have to scan every pixel in the array to calculate these
values.

[Nathan Orloff]

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