function fp = fitgaussian2D(m,cx,cy,tol);
%% m = image
%% tol = fitting tolerance

%Here we prepare a set of options to do our fits
options = optimset('Display','off','TolFun',tol,'LargeScale','off');
theta = 0;
background=0;

%We determine the size of an image
[sizey sizex] = size(m);

% check guessed center position
cx=check_bound(cx,1,sizex);
cy=check_bound(cy,1,sizey);
cy=round(cy);
cx=round(cx);

% guessing the widths of the distribution
[tx,ty,sx,sy] = centerofmass(m);
%peak value
peakAmpl = max(max(m));

% improve size and position along x
mx = m(cy,:);
x1D = 1:sizex;
ip1D = [cx,sx,peakAmpl,background];
fp1D = fminunc(@fitGaussian1D,ip1D,options,mx,x1D);

% update guessed values
background = fp1D(4);
cx = round(fp1D(1));
peakAmpl = fp1D(3);

% improve size and position along y
my = m(:,cx)';
y1D = 1:sizey;
ip1D = [cy,sy,peakAmpl,background];
fp1D = fminunc(@fitGaussian1D,ip1D,options,my,y1D);

% update guessed values
cy = round(fp1D(1));
sy = fp1D(2);


[X,Y] = meshgrid(1:sizex,1:sizey);

%Lower and upper limits for the 2D fit
LB = [1,1,1,1,-inf,-pi/2,-inf];
UB = [sizex,sizey,sizex*10,sizey*10,inf,pi/2,inf];

% Here we do the 2D fit
options = optimset('lsqnonlin');
options.Display = 'iter';
options.TolFun = tol;
initpar = [cx,cy,sx,sy,peakAmpl,theta,background];
%this is to avoid errors when doing fits on the image full of zeros
if isnan(sum(sum(Gaussian2D(initpar,X,Y)-m))) < 1
    % Here is the actual fit
    fp = fit_image(@Gaussian2D, initpar, X, Y, m, LB, UB, options);
else
    fp = [0,0,0,0,0,0,0];
end;