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
|
function fp = fitgaussian2D(m,cx,cy,tol);
%% m = image
%% tol = fitting tolerance
%Here we prepare a set of options to do our fits
options = optimset('Display','off','TolFun',tol,'LargeScale','off');
theta = 0;
background=0;
%We determine the size of an image
[sizey sizex] = size(m);
% check guessed center position
cx=check_bound(cx,1,sizex);
cy=check_bound(cy,1,sizey);
cy=round(cy);
cx=round(cx);
% guessing the widths of the distribution
[tx,ty,sx,sy] = centerofmass(m);
%peak value
peakAmpl = max(max(m));
% improve size and position along x
mx = m(cy,:);
x1D = 1:sizex;
ip1D = [cx,sx,peakAmpl,background];
fp1D = fminunc(@fitGaussian1D,ip1D,options,mx,x1D);
% update guessed values
background = fp1D(4);
cx = round(fp1D(1));
peakAmpl = fp1D(3);
% improve size and position along y
my = m(:,cx)';
y1D = 1:sizey;
ip1D = [cy,sy,peakAmpl,background];
fp1D = fminunc(@fitGaussian1D,ip1D,options,my,y1D);
% update guessed values
cy = round(fp1D(1));
sy = fp1D(2);
[X,Y] = meshgrid(1:sizex,1:sizey);
%Lower and upper limits for the 2D fit
LB = [1,1,1,1,-inf,-pi/2,-inf];
UB = [sizex,sizey,sizex*10,sizey*10,inf,pi/2,inf];
% Here we do the 2D fit
options = optimset('lsqnonlin');
options.Display = 'iter';
options.TolFun = tol;
initpar = [cx,cy,sx,sy,peakAmpl,theta,background];
%this is to avoid errors when doing fits on the image full of zeros
if isnan(sum(sum(Gaussian2D(initpar,X,Y)-m))) < 1
% Here is the actual fit
fp = fit_image(@Gaussian2D, initpar, X, Y, m, LB, UB, options);
else
fp = [0,0,0,0,0,0,0];
end;
|
