Gauss fit takes care of initial values if they are not provided

1 files changed, 39 insertions, 22 deletions

diff --git a/Gauss2DRotFit.m b/Gauss2DRotFit.m
index a669723..84af35e 100644
--- a/Gauss2DRotFit.m
+++ b/Gauss2DRotFit.m
@@ -1,12 +1,12 @@
-function [fitresult, gof] = gauss2DRotFit(im,varargin)
+function [fitresult, gof, fout] = gauss2DRotFit(im,varargin)
 %gauss2DRotFit(im,fop,title)
-%  Create a fit.
+%  Create a 2D Gaussian  fit (distribution might be squeezed and rotated).
 %  Input :
 %   
 %      im  : input image of 2D gaussian
 %      fop : fitoptions - a structure with 3 fields : Lower,StartPoint and
 %            Upper. Each of the fields is an array of 7 values : 
-%            offset
+%             offset/background
 %             amplitude of the 2D gaussian
 %             centroid X of the 2D gaussian
 %             centroid Y of the 2D gaussian
@@ -18,12 +18,19 @@ function [fitresult, gof] = gauss2DRotFit(im,varargin)
 %      fitresult : a fit object representing the fit.
 %      gof : structure with goodness-of fit info.
 %
+%  Created by:  Steffen B. Petersen (12 Dec 2012)
+%  Modified by: Eugeniy E. Mikhailov (Nov 2014)
+
+ 
 
 if nargin==0
     errordlg('gauss2DRotFit has to be called with at least 1 argument');
     return;
 end
 [sx,sy,sz]=size(im);
+if sz==1
+    im=double(im);
+end
 if sz==3
     im=rgb2gray(im);
 end
@@ -35,10 +42,13 @@ mystr='Untitled fit';
 [xData, yData, zData] = prepareSurfaceData( X,Y, im );
 
 % Set up fittype and options.
-ft = fittype( 'a + b*exp(-(((x-c1)*cosd(t1)+(y-c2)*sind(t1))/w1)^2-((-(x-c1)*sind(t1)+(y-c2)*cosd(t1))/w2)^2)', 'independent', {'x', 'y'}, 'dependent', 'z' );
+% Note we will fit the expression for intensity (not field) of the Gaussian beam,
+% thus extra 2 in exponents.
+ft = fittype( 'a + b*exp(-2*(((x-c1)*cosd(t1)+(y-c2)*sind(t1))/w1)^2 - 2*((-(x-c1)*sind(t1)+(y-c2)*cosd(t1))/w2)^2)', 'independent', {'x', 'y'}, 'dependent', 'z' );
 
 
 opts = fitoptions( ft );
+flagMakeStartPointGuess = true;
 if numel(varargin)>0
     s=varargin{1};
     if isfield(s,'Lower')
@@ -46,11 +56,24 @@ if numel(varargin)>0
     end
     if isfield(s,'StartPoint')
         opts.StartPoint=s.StartPoint;
+        flagMakeStartPointGuess = false;
     end
     if isfield(s,'Upper')
         opts.Upper=s.Upper;
     end
 end
+if ( flagMakeStartPointGuess )
+    % let's attempt to guess starting point
+    a0 = min(zData);
+    b0 = max(zData) - a0;
+    cx0=mean(xData.*zData)/mean(zData);
+    cy0=mean(yData.*zData)/mean(zData);
+    t=0;
+    % see above note about intensity of Gaussian beam
+    w10=sqrt(2* mean((xData-cx0).^2.*zData)/mean(zData) );
+    w20=sqrt(2* mean((yData-cy0).^2.*zData)/mean(zData) );
+    opts.StartPoint= [a0, b0, cx0, cy0, t, w10, w20]
+end
 if numel(varargin)==2
     mystr=varargin{2};
 end
@@ -65,32 +88,26 @@ opts.Display = 'Off';
 opts
 gof
 fout
-% Create a figure for the plots.
-figure( 1);
-set(gcf,'Name', '2D Gaussian fit 1' );
-set(gcf,'Position',[10,10,900,600]);
 
-% Plot fit with data.
+imgfit = fitresult(xData, yData);
+imgfit = reshape(imgfit, sx, sy);
+
 
-h = plot( fitresult, [xData, yData], zData);
-legend( h, mystr, 'im vs. X, Y',  'Location', 'NorthEast' );
-% Label axes
+% Create a figure for the plots.
+figure( 1);
+imagesc(1:sy, 1:sx, imgfit);
+title('Fitted image');
 xlabel( 'X' );
 ylabel( 'Y' );
-zlabel( 'im' );
-grid on
+zlabel( 'image fitted' );
 colormap(jet);
 
 % Plot residuals.
 figure(2);
-set(gcf,'Name', 'Gaussian fit Residuals');
-set(gcf,'Position',[100,100,900,600]);
-h = plot( fitresult, [xData, yData], zData, 'Style', 'Residual');
-legend( h, [mystr  ' - residuals'], 'Location', 'NorthEast' );
-% Label axes
+imagesc(1:sy, 1:sx, imgfit-im);
+title('Residuals of fitted image');
 xlabel( 'X' );
 ylabel( 'Y' );
-zlabel( 'im' );
-grid on
-figure(1);
+zlabel( 'residuals' );
+colormap(jet);