Code covered by the BSD License

# Ellipse Detection Using 1D Hough Transform

### Martin Simonovsky (view profile)

• 1 file
• 4.16667

28 Nov 2011 (Updated )

An efficient ellipse detector based on Hough voting

ellipseDetection(img, params)
function bestFits = ellipseDetection(img, params)
% ellipseDetection: Ellipse detection
%
% Overview:
% --------
% Fits an ellipse by examining all possible major axes (all pairs of points) and
% getting the minor axis using Hough transform. The algorithm complexity depends on
% the number of valid non-zero points, therefore it is beneficial to provide as many
% restrictions in the "params" input arguments as possible if there is any prior
%
% The code is reasonably fast due to (optional) randomization and full code vectorization.
% However, as the algorithm needs to compute pairwise point distances, it can be quite memory
% intensive. If you get out of memory errors, either downsample the input image or somehow
% decrease the number of non-zero points in it.
% It can deal with big amount of noise but can have severe problem with occlusions (major axis
% end points need to be visible)
%
% Input arguments:
% --------
% img
%   - One-channel input image (greyscale or binary).
% params
%   - Parameters of the algorithm:
%       minMajorAxis: Minimal length of major axis accepted.
%       maxMajorAxis: Maximal length of major axis accepted.
%        rotation, rotationSpan: Specification of restriction on the angle of the major axis in degrees.
%                                If rotationSpan is in (0,90), only angles within [rotation-rotationSpan,
%                               rotation+rotationSpan] are accepted.
%       minAspectRatio: Minimal aspect ratio of an ellipse (in (0,1))
%       randomize: Subsampling of all possible point pairs. Instead of examining all N*N pairs, runs
%                  only on N*randomize pairs. If 0, randomization is turned off.
%       numBest: Top numBest to return
%       uniformWeights: Used to prefer some points over others. If false, accumulator points are weighted
%                       by their grey intensity in the image. If true, the input image is regarded as binary.
%       smoothStddev: In order to provide more stability of the solution, the accumulator is convolved with
%                     a gaussian kernel. This parameter specifies its standard deviation in pixels.
%
% Return value:
% --------
% Returns a matrix of best fits. Each row (there are params.numBest of them) contains six elements:
% [x0 y0 a b alpha score] being the center of the ellipse, its major and minor axis, its angle in degrees and score.
%
% Based on:
% --------
% - "A New Efficient Ellipse Detection Method" (Yonghong Xie Qiang , Qiang Ji / 2002)
% - random subsampling inspired by "Randomized Hough Transform for Ellipse Detection with Result Clustering"
%   (CA Basca, M Talos, R Brad / 2005)
%
% Update log:
% --------
% 1.1: More memory efficient code, better documentation, more parameters, more solutions possible, example code.
% 1.0: Initial version
%
%
% Author: Martin Simonovsky
% e-mail: <mys007@seznam.cz>
% Release: 1.1
% Release date: 25.7.2013
%
%
% --------
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% default values
if nargin==1; params=[]; end
% - parameters to contrain the search
if ~isfield(params,'minMajorAxis');     params.minMajorAxis = 10; end
if ~isfield(params,'maxMajorAxis');     params.maxMajorAxis = 200; end
if ~isfield(params,'rotation');            params.rotation = 0; end
if ~isfield(params,'rotationSpan');        params.rotationSpan = 0; end
if ~isfield(params,'minAspectRatio');    params.minAspectRatio = 0.1; end
if ~isfield(params,'randomize');        params.randomize = 2; end
% - others
if ~isfield(params,'numBest');            params.numBest = 3; end
if ~isfield(params,'uniformWeights');   params.uniformWeights = true; end
if ~isfield(params,'smoothStddev');        params.smoothStddev = 1; end

eps = 0.0001;
bestFits = zeros(params.numBest,6);
params.rotationSpan = min(params.rotationSpan, 90);
H = fspecial('gaussian', [params.smoothStddev*6 1], params.smoothStddev);

[Y,X]=find(img);
Y = single(Y); X = single(X);
N = length(Y);

fprintf('Possible major axes: %d * %d = %d\n', N, N, N*N);

% compute pairwise distances between points (memory intensive!) and filter
% TODO: do this block-wise, just appending the filtered results (I,J)
distsSq = bsxfun(@minus,X,X').^2 + bsxfun(@minus,Y,Y').^2;
[I,J] = find(distsSq>=params.minMajorAxis^2 & distsSq<=params.maxMajorAxis^2);
idx = I<J;
I = uint32(I(idx)); J = uint32(J(idx));

fprintf('..after distance constraint: %d\n', length(I));

% compute pairwise angles and filter
if params.rotationSpan>0
tangents = (Y(I)-Y(J)) ./ (X(I)-X(J));
tanLo = tand(params.rotation-params.rotationSpan);
tanHi = tand(params.rotation+params.rotationSpan);
if tanLo<tanHi
idx = tangents > tanLo & tangents < tanHi;
else
idx = tangents > tanLo | tangents < tanHi;
end
I = I(idx); J = J(idx);
fprintf('..after angular constraint: %d\n', length(I));
else
fprintf('..angular constraint not used\n');
end

npairs = length(I);

% compute random choice and filter
if params.randomize>0
perm = randperm(npairs);
pairSubset = perm(1:min(npairs,N*params.randomize));
clear perm;
fprintf('..after randomization: %d\n', length(pairSubset));
else
pairSubset = 1:npairs;
end

% check out all hypotheses
for p=pairSubset
x1=X(I(p)); y1=Y(I(p));
x2=X(J(p)); y2=Y(J(p));

%compute center & major axis
x0=(x1+x2)/2; y0=(y1+y2)/2;
aSq = distsSq(I(p),J(p))/4;
thirdPtDistsSq = (X-x0).^2 + (Y-y0).^2;
K = thirdPtDistsSq <= aSq; % (otherwise the formulae in paper do not work)

%get minor ax propositions for all other points
fSq = (X(K)-x2).^2 + (Y(K)-y2).^2;
cosTau = (aSq + thirdPtDistsSq(K) - fSq) ./ (2*sqrt(aSq*thirdPtDistsSq(K)));
cosTau = min(1,max(-1,cosTau)); %inexact float arithmetic?!
sinTauSq = 1 - cosTau.^2;
b = sqrt( (aSq * thirdPtDistsSq(K) .* sinTauSq) ./ (aSq - thirdPtDistsSq(K) .* cosTau.^2 + eps) );

%proper bins for b
idxs = ceil(b+eps);

if params.uniformWeights
weights = 1;
else
weights = img(sub2ind(size(img),Y(K),X(K)));
end
accumulator = accumarray(idxs, weights, [params.maxMajorAxis 1]);

%a bit of smoothing and finding the most busy bin
accumulator = conv(accumulator,H,'same');
accumulator(1:ceil(sqrt(aSq)*params.minAspectRatio)) = 0;
[score, idx] = max(accumulator);

%keeping only the params.numBest best hypothesis (no non-maxima suppresion)
if (bestFits(end,end) < score)
bestFits(end,:) = [x0 y0 sqrt(aSq) idx atand((y1-y2)/(x1-x2)) score];
if params.numBest>1
[~,si]=sort(bestFits(:,end),'descend');
bestFits = bestFits(si,:);
end
end
end
end