initial commit

1 files changed, 101 insertions, 0 deletions

diff --git a/fitgaussian2D.m b/fitgaussian2D.m
new file mode 100644
index 0000000..a6f3d10
--- /dev/null
+++ b/fitgaussian2D.m
@@ -0,0 +1,101 @@
+%% a function to fit a thermal cloud 2-D
+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);
+% 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,:);
+
+
+%Here we do the 1D horizontal fit
+%This is to give better initial parameters to our 2D fit
+x1D = 1:sizex;
+%Initiall parameters
+ip1D = [cx,sx,pOD,background];
+fp1D = fminunc(@fitGaussian1D,ip1D,options,mx,x1D);
+
+cx = fp1D(1);
+sx = fp1D(2);
+PeakOD = fp1D(3);
+background = fp1D(4);
+
+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)';
+
+y1D = 1:sizey;
+ip1D = [cy,sy,PeakOD,background];
+%1D vertical fit
+fp1D = fminunc(@fitGaussian1D,ip1D,options,my,y1D);
+
+cy = 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
+LB = [1,1,1,1,0,-1.58,0];
+UB = [sizex,sizey,sizex*10,sizey*10,inf,1.58,inf];
+
+% Here we do the 2D fit
+options = optimset('lsqnonlin');
+options.Display = 'iter';
+options.TolFun = tol;
+initpar = [cx,cy,sx,sy,PeakOD,theta,background];
+%this is to avoid errors when doing fits on the image full of zeros
+if isnan(sum(sum(Gaussian2Dff(initpar,X,Y,m)))) < 1
+    % Here is the actual fit
+    fp = lsqnonlin( @(pt) Gaussian2Dff(pt,X,Y,m),initpar,LB,UB,options);
+else
+    fp = [0,0,0,0,0,0,0];
+end;
+
+
+
+% 2D Gauss function
+function [z] = Gaussian2Dff(p,X,Y,m);
+
+cx = p(1);
+cy = p(2);
+wx = p(3);
+wy = p(4);
+amp = p(5);
+theta = p(6);
+background=p(7);
+Xn = (X-cx)*cos(theta) - (Y-cy)*sin(theta);
+Yn = (X-cx)*sin(theta) + (Y-cy)*cos(theta);
+
+z = amp*(exp(-2*((Xn).^2./(wx^2)+(Yn).^2./(wy^2)))) + background - m;
\ No newline at end of file