--- /dev/null
+% Hough Transform for circles.
+% img : logical image, zeros are the circles edges
+% radius_range is a vector : [min, max]
+function [acc_votes, acc_radius1, acc_radius2, acc_alpha] = HTEllipse2(x_dir, y_dir, edges, radius_range)
+
+ [rows, columns] = size(edges);
+
+ acc_votes = zeros(rows, columns);
+ acc_radius1 = zeros(rows, columns);
+ acc_radius2 = zeros(rows, columns);
+ acc_alpha = zeros(rows, columns);
+
+ alpha_range = linspace(0, pi / 2, 8);
+ alpha_range(end) = [];
+% alpha_range = pi/8;
+
+ indexes_edges = find(edges)';
+ s = size(edges);
+
+ for alpha = alpha_range
+ rot_mat = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)];
+ rot_mat_inv = [cos(-alpha) -sin(-alpha); sin(-alpha) cos(-alpha)];
+ for r1 = radius_range(1):radius_range(2)
+ for r2 = radius_range(1):radius_range(2)
+ acc = zeros(rows, columns);
+ for i = indexes_edges
+ [y, x] = ind2sub(s, i);
+ v = [-x_dir(i); -y_dir(i)];
+ u = rot_mat_inv * v;
+ if u(1) == 0
+ if u(2) > 0
+ theta = pi / 2.;
+ else
+ theta = 3. * pi / 2.;
+ end
+ elseif u(2) == 0
+ if u(1) > 0
+ theta = 0;
+ else
+ theta = pi;
+ end
+ elseif u(2) < 0
+ t = u(2) / u(1);
+ theta = 2*pi + 2 * atan(...
+ (-r1 - r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
+ /(r2*t));
+ else
+ t = u(2) / u(1);
+ theta = 2 * atan(...
+ (-r1 + r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
+ /(r2*t));
+ end
+ C = round([x; y] - rot_mat * [r1*cos(theta); r2*sin(theta)]);
+ if C(1) < columns && C(1) > 0 && C(2) < rows && C(2) > 0
+ acc(C(2),C(1)) = acc(C(2),C(1)) + 1;
+ end
+ end
+
+ for x = 1:columns
+ for y = 1:rows
+ if acc(y, x) > acc_votes(y, x)
+ acc_votes(y, x) = acc(y, x);
+ acc_radius1(y, x) = r1;
+ acc_radius2(y, x) = r2;
+ acc_alpha(y, x) = alpha;
+ end
+ end
+ end
+ end
+ end
+ end
+
+ acc_votes = imgaussfilt(acc_votes, 1.0);
+end
+