Update KdTree tests.
[master-thesis.git] / src / Tests_hough / HTEllipse2.m
1 % Hough Transform for circles.
2 % img : logical image, zeros are the circles edges
3 % radius_range is a vector : [min, max]
4 function [acc_votes, acc_radius1, acc_radius2, acc_alpha] = HTEllipse2(x_dir, y_dir, edges, radius_range)
5
6     [rows, columns] = size(edges);
7     
8     acc_votes = zeros(rows, columns);
9     acc_radius1 = zeros(rows, columns);
10     acc_radius2 = zeros(rows, columns);
11     acc_alpha = zeros(rows, columns);
12     
13     alpha_range = linspace(0, pi / 2, 8);
14     alpha_range(end) = [];
15 %     alpha_range = pi/8;
16         
17     indexes_edges = find(edges)';
18     s = size(edges);
19     
20     for alpha = alpha_range
21         rot_mat = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)];
22         rot_mat_inv = [cos(-alpha) -sin(-alpha); sin(-alpha) cos(-alpha)];      
23         for r1 = radius_range(1):radius_range(2)
24             for r2 = radius_range(1):radius_range(2)
25                 acc = zeros(rows, columns);
26                 for i = indexes_edges
27                     [y, x] = ind2sub(s, i);
28                     v = [-x_dir(i); -y_dir(i)];
29                     u = rot_mat_inv * v;
30                     if u(1) == 0
31                         if u(2) > 0
32                             theta = pi / 2.;
33                         else
34                             theta = 3. * pi / 2.;
35                         end
36                     elseif u(2) == 0
37                         if u(1) > 0
38                             theta = 0;
39                         else
40                             theta = pi;
41                         end
42                     elseif u(2) < 0
43                         t = u(2) / u(1);
44                         theta = 2*pi + 2 * atan(...
45                             (-r1 - r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
46                             /(r2*t));
47                     else
48                         t = u(2) / u(1);
49                         theta = 2 * atan(...
50                             (-r1 + r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
51                             /(r2*t));
52                     end
53                     C = round([x; y] - rot_mat * [r1*cos(theta); r2*sin(theta)]);
54                     if C(1) < columns && C(1) > 0 && C(2) < rows && C(2) > 0
55                         acc(C(2),C(1)) = acc(C(2),C(1)) + 1;
56                     end
57                 end
58
59                 for x = 1:columns
60                     for y = 1:rows
61                         if acc(y, x) > acc_votes(y, x)
62                             acc_votes(y, x) = acc(y, x);
63                             acc_radius1(y, x) = r1;
64                             acc_radius2(y, x) = r2;
65                             acc_alpha(y, x) = alpha;
66                         end
67                     end
68                 end      
69             end
70         end
71     end
72     
73     acc_votes = imgaussfilt(acc_votes, 1.0);
74 end
75