enclosingCircle

ENCLOSINGCIRCLE Find the minimum circle enclosing a set of points.

   CIRCLE = enclosingCircle(POINTS);
   computes the circle CIRCLE=[xc yc r] which encloses all the points POINTS
   given as a N-by-2 array.


   Rewritten from a file from
           Yazan Ahed (yash78@gmail.com)

   which was rewritten from a Java applet by Shripad Thite:
   http://heyoka.cs.uiuc.edu/~thite/mincircle/

   See also:
   circles2d, points2d, boxes2d, circumCircle

0001 function circle = enclosingCircle(pts)
0002 %ENCLOSINGCIRCLE Find the minimum circle enclosing a set of points.
0003 %
0004 %   CIRCLE = enclosingCircle(POINTS);
0005 %   computes the circle CIRCLE=[xc yc r] which encloses all the points POINTS
0006 %   given as a N-by-2 array.
0007 %
0008 %
0009 %   Rewritten from a file from
0010 %           Yazan Ahed (yash78@gmail.com)
0011 %
0012 %   which was rewritten from a Java applet by Shripad Thite:
0013 %   http://heyoka.cs.uiuc.edu/~thite/mincircle/
0014 %
0015 %   See also:
0016 %   circles2d, points2d, boxes2d, circumCircle
0017 %
0018 
0019 % ------
0020 % Author: David Legland
0021 % e-mail: david.legland@nantes.inra.fr
0022 % created the 07/07/2005.
0023 % Copyright 2012 INRA - Cepia Software Platform.
0024 
0025 
0026 % works on convex hull: it is faster
0027 pts = pts(convhull(pts(:,1), pts(:,2)), :);
0028 
0029 % call the recursive function
0030 circle = recurseCircle(size(pts, 1), pts, 1, zeros(3, 2));
0031 
0032 
0033 
0034 function circ = recurseCircle(n, p, m, b)
0035 %    n: number of points given
0036 %    m: an argument used by the function. Always use 1 for m.
0037 %    bnry: an argument (3x2 array) used by the function to set the points that
0038 %          determines the circle boundry. You have to be careful when choosing this
0039 %          array's values. I think the values should be somewhere outside your points
0040 %          boundary. For my case, for example, I know the (x,y) I have will be something
0041 %          in between (-5,-5) and (5,5), so I use bnry as:
0042 %                       [-10 -10
0043 %                        -10 -10
0044 %                        -10 -10]
0045 
0046 
0047 if m == 4
0048     circ = createCircle(b(1,:), b(2,:), b(3,:));
0049     return;
0050 end
0051 
0052 circ = [Inf Inf 0];
0053 
0054 if m == 2
0055     circ = [b(1,1:2) 0];
0056 elseif m == 3
0057     c = (b(1,:) + b(2,:))/2;
0058     circ = [c distancePoints(b(1,:), c)];
0059 end
0060 
0061 
0062 for i = 1:n
0063     if distancePoints(p(i,:), circ(1:2)) > circ(3)
0064         if sum(b(:,1)==p(i,1) & b(:,2)==p(i,2)) == 0
0065             b(m,:) = p(i,:);
0066             circ = recurseCircle(i, p, m+1, b);
0067         end
0068     end
0069 end
0070