Subversion Repositories seema-scanner

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.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 (no non-maxima suppresion)
        if (bestFits(end,end) < score)
            bestFits(end,:) = [x0 y0 sqrt(aSq) (idx-1) atand((y1-y2)/(x1-x2)) score];
            if params.numBest>1
                [~,si]=sort(bestFits(:,end),'descend');
                bestFits = bestFits(si,:);
            end
        end
    end
end