1 files changed, 23 insertions, 44 deletions

diff --git a/fitgaussian2D.m b/fitgaussian2D.m
index c3a7051..82fc8fa 100644
--- a/fitgaussian2D.m
+++ b/fitgaussian2D.m
@@ -1,4 +1,3 @@
-%% a function to fit a thermal cloud 2-D
 function fp = fitgaussian2D(m,cx,cy,tol);
 %% m = image
 %% tol = fitting tolerance
@@ -10,60 +9,40 @@ background=0;
 
 %We determine the size of an image
 [sizey sizex] = size(m);
-% guessing the widths
-[tx,ty,sx,sy] = centerofmass(m);
-sx = 2.*sx;
-sy = 2.*sy;
-%peak value
-pOD = max(max(m));
 
-%This is to avoid any possible access to some sells outside of the matrix
-tmx = round(cy);
-if isnan(tmx) > 0
-    tmx = sizey/2;
-end;
-if round(cy) > sizey
-    tmx = sizey;
-end;
-if round(cy) < 1
-    tmx = 1;
-end;
-mx = m(tmx,:);
+% 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));
 
-%Here we do the 1D horizontal fit
-%This is to give better initial parameters to our 2D fit
+% improve size and position along x
+mx = m(cy,:);
 x1D = 1:sizex;
-%Initiall parameters
-ip1D = [cx,sx,pOD,background];
+ip1D = [cx,sx,peakAmpl,background];
 fp1D = fminunc(@fitGaussian1D,ip1D,options,mx,x1D);
 
-cx = fp1D(1);
-sx = fp1D(2);
-PeakOD = fp1D(3);
+% update guessed values
 background = fp1D(4);
+cx = round(fp1D(1));
+peakAmpl = fp1D(3);
 
-tmy = round(cx);
-if isnan(tmy) > 0
-    tmy = sizex/2;
-end;
-if round(cx) > sizex
-    tmy = sizex;
-end;
-if round(cx) < 1
-    tmy = 1;
-end;
-my = m(:,tmy)';
-
+% improve size and position along y
+my = m(:,cx)';
 y1D = 1:sizey;
-ip1D = [cy,sy,PeakOD,background];
-%1D vertical fit
+ip1D = [cy,sy,peakAmpl,background];
 fp1D = fminunc(@fitGaussian1D,ip1D,options,my,y1D);
 
-cy = fp1D(1);
+% update guessed values
+cy = round(fp1D(1));
 sy = fp1D(2);
-PeakOD = fp1D(3);
-background = fp1D(4);
+
+
 [X,Y] = meshgrid(1:sizex,1:sizey);
 
 %Lower and upper limits for the 2D fit
@@ -74,7 +53,7 @@ UB = [sizex,sizey,sizex*10,sizey*10,inf,pi/2,inf];
 options = optimset('lsqnonlin');
 options.Display = 'iter';
 options.TolFun = tol;
-initpar = [cx,cy,sx,sy,PeakOD,theta,background];
+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