diff options
| -rw-r--r-- | Gauss2DRotFit.m | 61 |
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); |