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
|
%% a function to fit a thermal cloud 2-D
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);
% guessing the widths
[tx,ty,sx,sy] = centerofmass(m);
sx = 2.*sx;
sy = 2.*sy;
%peak value
pOD = max(max(m));
%This is to avoid any possible access to some sells outside of the matrix
tmx = round(cy);
if isnan(tmx) > 0
tmx = sizey/2;
end;
if round(cy) > sizey
tmx = sizey;
end;
if round(cy) < 1
tmx = 1;
end;
mx = m(tmx,:);
%Here we do the 1D horizontal fit
%This is to give better initial parameters to our 2D fit
x1D = 1:sizex;
%Initiall parameters
ip1D = [cx,sx,pOD,background];
fp1D = fminunc(@fitGaussian1D,ip1D,options,mx,x1D);
cx = fp1D(1);
sx = fp1D(2);
PeakOD = fp1D(3);
background = fp1D(4);
tmy = round(cx);
if isnan(tmy) > 0
tmy = sizex/2;
end;
if round(cx) > sizex
tmy = sizex;
end;
if round(cx) < 1
tmy = 1;
end;
my = m(:,tmy)';
y1D = 1:sizey;
ip1D = [cy,sy,PeakOD,background];
%1D vertical fit
fp1D = fminunc(@fitGaussian1D,ip1D,options,my,y1D);
cy = fp1D(1);
sy = fp1D(2);
PeakOD = fp1D(3);
background = fp1D(4);
[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,PeakOD,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;
|
