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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
| function ellipse_t = fit_ellipse( x,y,axis_handle )
% initialize
orientation_tolerance = 1e-3;
% empty warning stack
warning( '' );
% prepare vectors, must be column vectors
x = x(:);
y = y(:);
% remove bias of the ellipse - to make matrix inversion more accurate. (will be added later on).
mean_x = mean(x);
mean_y = mean(y);
x = x-mean_x;
y = y-mean_y;
% the estimation for the conic equation of the ellipse
X = [x.^2, x.*y, y.^2, x, y ];
a = sum(X)/(X'*X);
% check for warnings
if ~isempty( lastwarn )
disp( 'stopped because of a warning regarding matrix inversion' );
ellipse_t = [];
return
end
% extract parameters from the conic equation
[a,b,c,d,e] = deal( a(1),a(2),a(3),a(4),a(5) );
% remove the orientation from the ellipse
if ( min(abs(b/a),abs(b/c)) > orientation_tolerance )
orientation_rad = 1/2 * atan( b/(c-a) );
cos_phi = cos( orientation_rad );
sin_phi = sin( orientation_rad );
[a,b,c,d,e] = deal(...
a*cos_phi^2 - b*cos_phi*sin_phi + c*sin_phi^2,...
0,...
a*sin_phi^2 + b*cos_phi*sin_phi + c*cos_phi^2,...
d*cos_phi - e*sin_phi,...
d*sin_phi + e*cos_phi );
[mean_x,mean_y] = deal( ...
cos_phi*mean_x - sin_phi*mean_y,...
sin_phi*mean_x + cos_phi*mean_y );
else
orientation_rad = 0;
cos_phi = cos( orientation_rad );
sin_phi = sin( orientation_rad );
end
% check if conic equation represents an ellipse
test = a*c;
switch (1)
case (test>0), status = '';
case (test==0), status = 'Parabola found'; warning( 'fit_ellipse: Did not locate an ellipse' );
case (test<0), status = 'Hyperbola found'; warning( 'fit_ellipse: Did not locate an ellipse' );
end
% if we found an ellipse return it's data
if (test>0)
% make sure coefficients are positive as required
if (a<0), [a,c,d,e] = deal( -a,-c,-d,-e ); end
% final ellipse parameters
X0 = mean_x - d/2/a;
Y0 = mean_y - e/2/c;
F = 1 + (d^2)/(4*a) + (e^2)/(4*c);
[a,b] = deal( sqrt( F/a ),sqrt( F/c ) );
long_axis = 2*max(a,b);
short_axis = 2*min(a,b);
% rotate the axes backwards to find the center point of the original TILTED ellipse
R = [ cos_phi sin_phi; -sin_phi cos_phi ];
P_in = R * [X0;Y0];
X0_in = P_in(1);
Y0_in = P_in(2);
else
% report an empty structure
ellipse_t = struct( ...
'a',[],...
'b',[],...
'phi',[],...
'X0',[],...
'Y0',[],...
'X0_in',[],...
'Y0_in',[],...
'long_axis',[],...
'short_axis',[],...
'status',status );
end
% check if we need to plot an ellipse with it's axes.
if (nargin>2) & ~isempty( axis_handle ) & (test>0)
% rotation matrix to rotate the axes with respect to an angle phi
R = [ cos_phi sin_phi; -sin_phi cos_phi ];
% the axes
ver_line = [ [X0 X0]; Y0+b*[-1 1] ];
horz_line = [ X0+a*[-1 1]; [Y0 Y0] ];
new_ver_line = R*ver_line;
new_horz_line = R*horz_line;
% the ellipse
theta_r = linspace(0,2*pi);
ellipse_x_r = X0 + a*cos( theta_r );
ellipse_y_r = Y0 + b*sin( theta_r );
rotated_ellipse = R * [ellipse_x_r;ellipse_y_r];
% pack ellipse into a structure
ellipse_t = struct( ...
'a',a,...
'b',b,...
'phi',orientation_rad,...
'X0',X0,...
'Y0',Y0,...
'X0_in',X0_in,...
'Y0_in',Y0_in,...
'long_axis',long_axis,...
'short_axis',short_axis,...
'status','',...
'rotated_ellipse',rotated_ellipse );
end |
reconstruction en 3D
et je travaille toujours sur le traitement des données brutes pour la reconstruction des ellipses.
Répondre avec citation
Envoyé par soldius87
Partager