Egauss2D_fit_act

%%%%%%%%
  To fit a 2D gaussian to the raster maps
  Use starting best guess values of:
   fit_params =
       [a1,sigmax,sigmay,x0,y0,theta]
  
  Expects input arrays from e.g. meshgrid XI, YI and ZI
   Code then makes into suitable vectors
  amp = amplitude
  sigmax = x-width
  sigmay = y-width
  x0 = offset in x
  y0 = offset in y0
  theta = angle of ellipse

 If no start points given in fit_params,  it will estimate from the data

 ACT 11/5/2012

%%%%%%%%%%%%%

0001 function [sf gof t] = gauss2D_fit_act(XI,YI,ZI,fit_params,chan,gfit,m)
0002 
0003 
0004 %%%%%%%%%
0005 %  To fit a 2D gaussian to the raster maps
0006 %  Use starting best guess values of:
0007 %   fit_params =
0008 %       [a1,sigmax,sigmay,x0,y0,theta]
0009 %
0010 %  Expects input arrays from e.g. meshgrid XI, YI and ZI
0011 %   Code then makes into suitable vectors
0012 %  amp = amplitude
0013 %  sigmax = x-width
0014 %  sigmay = y-width
0015 %  x0 = offset in x
0016 %  y0 = offset in y0
0017 %  theta = angle of ellipse
0018 %
0019 % If no start points given in fit_params,  it will estimate from the data
0020 %
0021 % ACT 11/5/2012
0022 %
0023 %%%%%%%%%%%%%%
0024 
0025 if(~exist('gfit'))
0026     gfit=[];
0027 end
0028 if(~exist('chan'))
0029     chan=1;
0030 end
0031 if(~exist('m')) % so we can count multiple raster fits
0032     m=1;
0033 end
0034 % Creat suitable vectors
0035 X = XI(:);
0036 Y = YI(:);
0037 Z = ZI(:);
0038 
0039 % Define the gaussian
0040 
0041 % %  General elliptical gaussian check code
0042 % % for theta = 0:pi/100:pi
0043 % %     a = (cos(theta)^2/2/sigma_x^2 + sin(theta)^2/2/sigma_y^2);
0044 % %     b = (-sin(2*theta)/4/sigma_x^2 + sin(2*theta)/4/sigma_y^2) ;
0045 % %     c = (sin(theta)^2/2/sigma_x^2 + cos(theta)^2/2/sigma_y^2);
0046 % %
0047 % %     [X, Y] = meshgrid(-5:.1:5, -5:.1:5);
0048 % %     Z = A*exp( - (a*(X-x0).^2 + 2*b*(X-x0).*(Y-y0) + c*(Y-y0).^2)) ;
0049 % %     surf(X,Y,Z);shading interp;view(-36,36);axis equal;drawnow
0050 % % end
0051 
0052 
0053     
0054    
0055 % Make sure we have a starting point for the I channels:
0056 
0057 %gauss2 = fittype('a1*exp( - ( (cos(theta)^2/2/sigma_x^2 + sin(theta)^2/2/sigma_y^2)*(x-x0).^2 + 2*(-sin(2*theta)/4/sigma_x^2 + sin(2*theta)/4/sigma_y^2)*(x-x0).*(y-y0) + (sin(theta)^2/2/sigma_x^2 + cos(theta)^2/2/sigma_y^2)*(y-y0).^2  ))',...
0058 %    'independent', {'x','y'},'dependent', 'z', 'coefficients',{'a1'},'problem',{'sigma_x','sigma_y','theta','x0','y0'});
0059 
0060 if(nargin<8)
0061     
0062     a1 = max(Z(:)); % height, determine from image. may want to subtract background
0063     sigma_x = 0.4; % guess width
0064     sigma_y = 0.4; % guess width
0065     x0 = mean(X); % guess position (center seems a good place to start)
0066     y0 = mean(Y);
0067     theta=0;
0068 end
0069 
0070 %Check for Nans in the data
0071 aa = isnan(Z);
0072 
0073 % Set some fit options
0074 options = fitoptions('Method','NonlinearLeastSquares');
0075 options.Exclude = aa;
0076 
0077 % Compute fit and plot figure
0078 chan
0079 if (chan==1 || chan==8)
0080     gauss2 = fittype('a1*exp( - ( (cosd(theta)^2/2/sigma_x^2 + sind(theta)^2/2/sigma_y^2)*(x-x0).^2 + 2*(-sind(2*theta)/4/sigma_x^2 + sind(2*theta)/4/sigma_y^2)*(x-x0).*(y-y0) + (sind(theta)^2/2/sigma_x^2 + cosd(theta)^2/2/sigma_y^2)*(y-y0).^2  ))',...
0081     'independent', {'x','y'},'dependent', 'z', 'coefficients',{'a1','sigma_x','sigma_y','theta','x0','y0'});
0082 
0083     options.StartPoint = [a1,sigma_x,sigma_y,x0,y0,theta]
0084     
0085     [sf gof t] = fit([X,Y],double(Z),gauss2,options);
0086 
0087 else
0088     disp('Fitting peak of polarization data only') % fix sigma_x, sigma_y, x0,y0,theta from I1 data.
0089 
0090     gauss2 = fittype('a1*exp( - ( (cosd(theta)^2/2/sigma_x^2 + sind(theta)^2/2/sigma_y^2)*(x-x0).^2 + 2*(-sind(2*theta)/4/sigma_x^2 + sind(2*theta)/4/sigma_y^2)*(x-x0).*(y-y0) + (sind(theta)^2/2/sigma_x^2 + cosd(theta)^2/2/sigma_y^2)*(y-y0).^2  ))',...
0091     'independent', {'x','y'},'dependent', 'z', 'coefficients',{'a1'},'problem',{'sigma_x','sigma_y','theta','x0','y0'});
0092 
0093     options.StartPoint = [a1];
0094     fixed_params={gfit{m,1}.sigma_x,gfit{m,1}.sigma_y,gfit{m,1}.theta,gfit{m,1}.x0,gfit{m,1}.y0} % use gaussian fit from the I1 channel to fix position and angle, then just fit amplitude
0095     
0096     [sf gof t] = fit([X,Y],double(Z),gauss2,options,'problem',fixed_params);
0097 
0098 end
0099 
0100 figure(20); clf; plot(sf,[X,Y],Z);
0101 
0102 %sf.x0 and sf.y0 is the center of gaussian.
0103 %sf.sigmax etc will get you the other parameters.