|
Re: Fitting to a Gaussian using nlinfit
Posted:
May 23, 2013 1:10 PM
|
|
This will probably do it.
function [fitresult, zfit, fiterr, zerr, resnorm] = fmgaussfit(xx,yy,zz) %fmgaussfit(xx,yy,zz) Custom Guassian Fit algorithm using lsqcurvefit % A 3D gaussian fitting routine with error propagation and uncertainties.
%% Condition the data [xData, yData, zData] = prepareSurfaceData( xx, yy, zz ); xyData = {xData,yData};
%% Set up the startpoint [amp, ind] = max(zData); % amp is the amplitude. xo = xData(ind); % guess that it is at the maximum yo = yData(ind); % guess that it is at the maximum ang = 45; % angle in degrees. sy = 1; sx = 1; zo = median(zData(:))-std(zData(:)); xmax = max(xData)+2; ymax = max(yData)+2; xmin = min(xData)-2; ymin = min(yData)-2;
%% Set up fittype and options. Lower = [0, 0, 0, 0, xmin, ymin, 0]; Upper = [Inf, 90, Inf, Inf, xmax, ymax, Inf]; % angles greater than 90 are redundant StartPoint = [amp, ang, sx, sy, xo, yo, zo];%[amp, sx, sxy, sy, xo, yo, zo];
tols = 1e-16; options = optimset('Algorithm','levenberg-marquardt',... 'Display','off',... 'MaxFunEvals',5e2,... 'MaxIter',5e2,... 'TolX',tols,... 'TolFun',tols,... 'TolCon',tols ,... 'UseParallel','always');
%% perform the fitting [fitresult,resnorm,residual] = ... lsqcurvefit(@gaussian2D,StartPoint,xyData,zData,Lower,Upper,options); [fiterr, zfit, zerr] = gaussian2Duncert(fitresult,residual,jacobian,xyData);
end
function z = gaussian2D(par,xy) z = par(7) + ... par(1)*exp(-(((xy{1}-par(5)).*cosd(par(2))+(xy{2}-par(6)).*sind(par(2)))./par(3)).^2-... ((-(xy{1}-par(5)).*sind(par(2))+(xy{2}-par(6)).*cosd(par(2)))./par(4)).^2); end
function [dpar,zf,dzf] = gaussian2Duncert(par,resid,xy) J = guassian2DJacobian(par,xy); parci = nlparci(par,resid,'Jacobian',J); dpar = (diff(parci,[],2)./2)'; [zf,dzf] = nlpredci(@gaussian2D,xy,par,resid,'Jacobian',J); end
function J = guassian2DJacobian(par,xy) x = xy{1}; y = xy{2}; J(:,1) = exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2); J(:,2) = -par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(3).^2 - (2.*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2); J(:,3) = (2.*par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2)./par(3)^3; J(:,4) = (2.*par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2)./par(4)^3; J(:,5) = par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*cosd(par(2)).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))))./par(3).^2 - (2.*sind(par(2)).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2); J(:,6) = par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*cosd(par(2)).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2 + (2.*sind(par(2)).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))))./par(3).^2); J(:,7) = ones(size(x)); end
|
|