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