1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
|
function [fitresult, gof, fout] = gauss2DRotFit(im,varargin)
%gauss2DRotFit(im,fop,title)
% 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/background
% amplitude of the 2D gaussian
% centroid X of the 2D gaussian
% centroid Y of the 2D gaussian
% angle of rotation for the 2D gaussian
% width X of the 2D gaussian
% width Y of the 2D gaussian
%
% Output:
% 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
[X,Y]=meshgrid(1:sy,1:sx);
mystr='Untitled fit';
[xData, yData, zData] = prepareSurfaceData( X,Y, im );
% Set up fittype and options.
% 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')
opts.Lower=s.Lower;
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
opts.Display = 'Off';
% Fit model to data.
[fitresult, gof,fout] = fit( [xData, yData], zData, ft, opts );
opts
gof
fout
imgfit = fitresult(xData, yData);
imgfit = reshape(imgfit, sx, sy);
% Create a figure for the plots.
figure( 1);
imagesc(1:sy, 1:sx, imgfit);
title('Fitted image');
xlabel( 'X' );
ylabel( 'Y' );
zlabel( 'image fitted' );
colormap(jet);
% Plot residuals.
figure(2);
imagesc(1:sy, 1:sx, imgfit-im);
title('Residuals of fitted image');
xlabel( 'X' );
ylabel( 'Y' );
zlabel( 'residuals' );
colormap(jet);
|
