------------------------------------------------------------------------------
fit_gauss2d function FitFun
Description: Fit a 2-D Gaussian to data of the form f(x,y).
f(x,y) = A*exp(-(a*(x-x0)^2+2*b*(x-x0)*(y-y0)+c*(y-y0)^2))
Input : - Vector describing the X axis of f(x,y).
- Vector describing the Y axis of f(x,y).
- Matrix or vector describing f(x,y).
- Errors in f(x,y). Default is 1 (equal weight).
If empty matrix then use default.
- Fit background B + A*exp(...) {'y' | 'n'}, default is 'y'.
- Linear least square method:
'chol' - Cholesky decomposition (default).
'orth' - Orthogonal decomposition.
- Fit elliptical gaussian {'y' | 'n'}, default is 'y'.
Output : - Structure containing the best fit solution.
The following fields are available:
.Par - Best fit parameters
[A; 1./(2.*(SigmaX.^2 + SigmaY.^2))]
.ParErr - Errors in best fit parameters.
.Chi2 - \chi^2 of fit.
.Dof - Number of degrees of freedom.
- Residuals matrix.
Tested : Matlab 7.10
By : Eran O. Ofek August 2010
URL : http://wise-obs.tau.ac.il/~eran/matlab.html
Example: % Construct a Gaussian with noise and fit:
A=0.01.*randn(100,100);
X = [1:1:100]; Y=[1:1:100]; [MatX,MatY]=meshgrid(X,Y);
A = A+5.*exp(-((MatX-45).^2 + (MatY-56).^2)./5)
[ParS,Resid]=fit_gauss2d(X,Y,Z,ErrZ,FitBack,Method);
Reliable: 2
------------------------------------------------------------------------------0001 function [ParS,Resid]=fit_gauss2d(X,Y,Z,ErrZ,FitBack,Method,Elliptical); 0002 %------------------------------------------------------------------------------ 0003 % fit_gauss2d function FitFun 0004 % Description: Fit a 2-D Gaussian to data of the form f(x,y). 0005 % f(x,y) = A*exp(-(a*(x-x0)^2+2*b*(x-x0)*(y-y0)+c*(y-y0)^2)) 0006 % Input : - Vector describing the X axis of f(x,y). 0007 % - Vector describing the Y axis of f(x,y). 0008 % - Matrix or vector describing f(x,y). 0009 % - Errors in f(x,y). Default is 1 (equal weight). 0010 % If empty matrix then use default. 0011 % - Fit background B + A*exp(...) {'y' | 'n'}, default is 'y'. 0012 % - Linear least square method: 0013 % 'chol' - Cholesky decomposition (default). 0014 % 'orth' - Orthogonal decomposition. 0015 % - Fit elliptical gaussian {'y' | 'n'}, default is 'y'. 0016 % Output : - Structure containing the best fit solution. 0017 % The following fields are available: 0018 % .Par - Best fit parameters 0019 % [A; 1./(2.*(SigmaX.^2 + SigmaY.^2))] 0020 % .ParErr - Errors in best fit parameters. 0021 % .Chi2 - \chi^2 of fit. 0022 % .Dof - Number of degrees of freedom. 0023 % - Residuals matrix. 0024 % Tested : Matlab 7.10 0025 % By : Eran O. Ofek August 2010 0026 % URL : http://wise-obs.tau.ac.il/~eran/matlab.html 0027 % Example: % Construct a Gaussian with noise and fit: 0028 % A=0.01.*randn(100,100); 0029 % X = [1:1:100]; Y=[1:1:100]; [MatX,MatY]=meshgrid(X,Y); 0030 % A = A+5.*exp(-((MatX-45).^2 + (MatY-56).^2)./5) 0031 % [ParS,Resid]=fit_gauss2d(X,Y,Z,ErrZ,FitBack,Method); 0032 % Reliable: 2 0033 %------------------------------------------------------------------------------ 0034 0035 Def.ErrZ = []; 0036 Def.FitBack = 'y'; 0037 Def.Method = 'chol'; 0038 Def.Elliptical = 'y'; 0039 if (nargin==3), 0040 ErrZ = Def.ErrZ; 0041 FitBack = Def.FitBack; 0042 Method = Def.Method; 0043 Elliptical = Def.Elliptical; 0044 elseif (nargin==4), 0045 FitBack = Def.FitBack; 0046 Method = Def.Method; 0047 Elliptical = Def.Elliptical; 0048 elseif (nargin==5), 0049 Method = Def.Method; 0050 Elliptical = Def.Elliptical; 0051 elseif (nargin==6), 0052 Elliptical = Def.Elliptical; 0053 elseif (nargin==7), 0054 % do nothing 0055 else 0056 error('Illegal number of input arguments'); 0057 end 0058 0059 0060 if (isempty(ErrZ)), 0061 ErrZ = 1; 0062 end 0063 SizeZ = size(Z); 0064 Nel = numel(Z); 0065 if (numel(ErrZ)==1), 0066 ErrZ = ErrZ.*ones(SizeZ); 0067 end 0068 0069 if (sum(SizeZ==1)>0), 0070 % Z is a vector 0071 VecX = X; 0072 VecY = Y; 0073 VecZ = Z; 0074 VecErrZ = ErrZ; 0075 else 0076 % Z is a matrix 0077 [MatX, MatY] = meshgrid(X,Y); 0078 VecX = reshape(MatX,[Nel, 1]); 0079 VecY = reshape(MatY,[Nel, 1]); 0080 VecZ = reshape(Z,[Nel, 1]); 0081 VecErrZ = reshape(ErrZ,[Nel, 1]); 0082 end 0083 0084 0085 % Z = A*exp(-(X.^2+Y.^2)./(2.*(SigmaX.^2 + SigmaY.^2))); 0086 % log(Z) = log(A) -(X.^2+Y.^2)./(2.*(SigmaX.^2 + SigmaY.^2)) 0087 0088 % transform to log(f(x,y)): 0089 LogVecErrZ = VecErrZ./VecZ; 0090 0091 switch lower(FitBack) 0092 case 'n' 0093 % no background - linear least squares 0094 LogVecZ = log(VecZ); 0095 0096 switch lower(Elliptical) 0097 case 'n' 0098 % Fit circular Gaussian 0099 % a*(X-X0)^2 + a*(Y-Y0)^2 0100 % = 0101 % X^2 (a) 0102 % Y^2 (a) 0103 % X (-2*a*X0) 0104 % Y (-2*a*Y0) 0105 % (a*X0^2 + a*Y0^2) 0106 %--- 0107 % A+(a*X0^2 + a*Y0^2), a ,-2*a*X0,-2*a*Y0 0108 H = [ones(Nel,1), -(VecX.^2+VecY.^2), VecX, VecY]; 0109 0110 case 'y' 0111 % Fit elliptical Gaussian 0112 % X^2 (a) 0113 % Y^2 (c) 0114 % X (-2*a*X0 - 2*b*Y0) 0115 % Y (-2*b*X0 - 2*c*Y0) 0116 % X*Y (2*b) 0117 % (a*X0^2 + c*Y0^2 + 2*b*X0*Y0) 0118 %--- 0119 % 0120 H = [ones(Nel,1), -VecX.^2, VecY.^2, VecX, VecY, VecX.*VecY]; 0121 0122 otherwise 0123 error('Unknown Elliptical option'); 0124 end 0125 0126 [Par,ParErr] = lscov(H,VecY,1./LogVecErrZ.^2,Method) 0127 0128 VecResid = exp(LogVecZ - H*Par); 0129 ParS.Par = Par; 0130 % ParS.Par(1) = exp(ParS.Par(1)); 0131 ParS.ParErr = ParErr; 0132 % ParS.ParErr(1) = ParS.ParErr(1).*exp(ParS.Par(1)); 0133 ParS.Chi2 = sum(VecResid.^2./(VecErrZ.^2)); 0134 Npar = length(ParS.Par); 0135 ParS.Dof = Nel - Npar; 0136 Resid = reshape(VecResid,SizeZ); 0137 0138 case 'y' 0139 % background - non-linear least squares 0140 % fminsearch... 0141 error('FitBack not implemented'); 0142 otherwise 0143 error('Unknown FitBack option'); 0144 end 0145 0146