Subversion Repositories seema-scanner

% ========================================================================================= 
% Author: Andre A. Boechat
% File: hough_ellipse.m
% Date: February 27, 2014, 10:12:22 AM
% Description: Implementation of the classical Hough Transform algorithm to detect 
% ellipses.
% 
% Reference: 
% [1] Y. Xie and Q. Ji, “A new efficient ellipse detection method,” Pattern
% Recognition, 2002. Proceedings. 16th …, vol. 2, pp. 957–960, 2002.
%
% [2] http://en.wikipedia.org/wiki/Hough_transform
% [3] http://en.wikipedia.org/wiki/Randomized_Hough_transform
% ========================================================================================= 
%
% Usage:
% [PARAMETERS] = hough_ellipse(IMG, MIN2A = 10, MIN_VOTES = 10)
%
% Inputs:
% IMG is the inputs image. Images composed of only edges are better.
% MIN2A is the minimum length of the major axis (default 10).
% MIN_VOTES is the minimum number of votes on a "b" value (half-length of the minor 
%   axis) to validate the existence of an ellipse.
%
% Outputs:
% PARAMETERS is the parameters of the best fitted ellipse on the image and is 
%   composed of [x0 y0 a b alpha] (ellipse's center, major and minor 
%   half-length axis and orientation, respectively).
% 
% 
% Possible improvement: to return a matrix with all best fitted ellipses on the 
% image. It adds another loop to the algorithm.

function [parameters] = hough_ellipse(img, min2a, min_votes)

    [width height] = size(img);
    %% Finding all nonzero pixels of the image, possible ellipse's pixels.
    [ys xs] = find(img);
    pixels = [xs ys];
    %% Accumulator for the minor axis' half-length. The indexes correspond to the
    %% possible b values. 
    %% TODO: the data structure can be improved (tree with the possible values?).
    acc = zeros(1, max(width, height)/2);
    %% ij1, ij2 are indexes of (x1, y1) and (x2, y2), following the reference [1].
    for ij1 = 1:(length(xs)-1)
        for ij2 = length(xs):-1:(ij1+1)
            x1 = pixels(ij1, 1);
            y1 = pixels(ij1, 2);
            x2 = pixels(ij2, 1);
            y2 = pixels(ij2, 2);
            d12 = norm([x1 y1] - [x2 y2]);
            acc = acc * 0;
            if  x1 - x2 && d12 > min2a
                %% Center
                x0 = (x1 + x2)/2;
                y0 = (y1 + y2)/2;
                %% Half-length of the major axis.
                a = d12/2;
                %% Orientation.
                alpha = atan((y2 - y1)/(x2 - x1));
                %% Distances between the two points and the center.
                d01 = norm([x1 y1] - [x0 y0]);
                d02 = norm([x2 y2] - [x0 y0]);
                for ij3 = 1:length(xs)
                    %% The third point must be a different point, obviously.
                    if (ij3 == ij1) && (ij3 == ij2)
                        continue;
                    end
                    x3 = pixels(ij3, 1);
                    y3 = pixels(ij3, 2);
                    d03 = norm([x3 y3] - [x0 y0]);
                    %% Distance f
                    if  d03 >= a
                        continue;
                    end
                    f = norm([x3 y3] - [x2 y2]);
                    %% Estimating the half-length of the minor axis.
                    cos2_tau = ((a^2 + d03^2 - f^2) / (2 * a * d03))^2;
                    sin2_tau = 1 - cos2_tau;
                    b = round(sqrt((a^2 * d03^2 * sin2_tau) /...
                                (a^2 - d03^2 * cos2_tau)));
                    %% Changing the score of the accumulator, if b is a valid value.
                    %% NOTE: the accumulator's length gives us the biggest expected
                    %% value for b, which means, in this current implementation,
                    %% we wouldn't detect ellipses whose half of minor axis is
                    %% greater than the image's size (look at the acc's declaration).
                    if b > 0 && b <= length(acc)
                        acc(b) = acc(b)+1;
                    end
                    end
                %% Taking the highest score.
                [sv si] = max(acc);
                if sv > min_votes
                    %% Ellipse detected!
                    %% The index si gives us the best b value.
                    parameters = [x0 y0 a si alpha];
                    return;
                end
                    end
                    end
                    end
    printf("No ellipses detected!\n");
end