Rev 210 | Blame | Compare with Previous | Last modification | View Log | RSS feed
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
% knowledge about the problem.
%
% 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.2: Added more efficient distance calculation, less memory-usage, non-maximum suppression.
% 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 = round(b)+1;
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
if (bestFits(end,end) < score)
%non-maximum suppression (based on center distance for now)
squaredDists = (bestFits(:,1)-x0).^2 + (bestFits(:,2)-y0).^2;
[minSqDist, minSqDistIdx] = min(squaredDists);
if(minSqDist < 20^2 && score > bestFits(minSqDistIdx, 6))
bestFits(minSqDistIdx,:) = [x0 y0 sqrt(aSq) (idx-1) atand((y1-y2)/(x1-x2)) score];
elseif(minSqDist >= 10^2)
bestFits(end,:) = [x0 y0 sqrt(aSq) (idx-1) atand((y1-y2)/(x1-x2)) score];
end
if params.numBest>1
[~,si]=sort(bestFits(:,end),'descend');
bestFits = bestFits(si,:);
end
end
end
end
Generated by GNU Enscript 1.6.6.