%%%%%%%% To fit a 2D gaussian to the raster maps Use starting best guess values of 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 If no start points given, it will estimate from the data ACT 11/5/2012 %%%%%%%%%%%%%
0001 function [sf gof t] = gauss2D_fit_act(XI,YI,ZI,a1,sigmax, sigmay,x0,y0) 0002 0003 0004 %%%%%%%%% 0005 % To fit a 2D gaussian to the raster maps 0006 % Use starting best guess values of 0007 % 0008 % Expects input arrays from e.g. meshgrid XI, YI and ZI 0009 % Code then makes into suitable vectors 0010 % amp = amplitude 0011 % sigmax = x-width 0012 % sigmay = y-width 0013 % x0 = offset in x 0014 % y0 = offset in y0 0015 % 0016 % If no start points given, it will estimate from the data 0017 % 0018 % ACT 11/5/2012 0019 % 0020 %%%%%%%%%%%%%% 0021 0022 % Creat suitable vectors 0023 X = XI(:); 0024 Y = YI(:); 0025 Z = ZI(:); 0026 0027 % Define the gaussian 0028 0029 % % General elliptical gaussian check code 0030 % % for theta = 0:pi/100:pi 0031 % % a = (cos(theta)^2/2/sigma_x^2 + sin(theta)^2/2/sigma_y^2); 0032 % % b = (-sin(2*theta)/4/sigma_x^2 + sin(2*theta)/4/sigma_y^2) ; 0033 % % c = (sin(theta)^2/2/sigma_x^2 + cos(theta)^2/2/sigma_y^2); 0034 % % 0035 % % [X, Y] = meshgrid(-5:.1:5, -5:.1:5); 0036 % % Z = A*exp( - (a*(X-x0).^2 + 2*b*(X-x0).*(Y-y0) + c*(Y-y0).^2)) ; 0037 % % surf(X,Y,Z);shading interp;view(-36,36);axis equal;drawnow 0038 % % end 0039 0040 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 ))',... 0041 'independent', {'x','y'},'dependent', 'z' ); 0042 0043 0044 0045 0046 % Make sure we have a starting point 0047 0048 if(nargin<8) 0049 0050 a1 = min(Z(:)); % height, determine from image. may want to subtract background 0051 sigma_x = 0.4; % guess width 0052 sigma_y = 0.4; % guess width 0053 x0 = mean(X); % guess position (center seems a good place to start) 0054 y0 = mean(Y); 0055 theta=0; 0056 end 0057 0058 % Set some fit options 0059 options = fitoptions('Method','NonlinearLeastSquares'); 0060 options.StartPoint = [a1,sigma_x,sigma_y,theta,x0,y0]; 0061 0062 %Check for Nans in the data 0063 aa = isnan(Z); 0064 options.Exclude = aa; 0065 0066 % Compute fit and plot figure 0067 0068 [sf gof t] = fit([X,Y],double(Z),gauss2,options); 0069 0070 0071 %figure(20); clf; plot(sf,[X,Y],Z); 0072 0073 %sf.x0 and sf.y0 is the center of gaussian. 0074 %sf.sigmax etc will get you the other parameters.