function [C,R] = incircle(x,y)
% incircle: compute the maximal in-circle of the polygonal convex hull of a set of points in the plane
%
% [C,R] = incircle(x,y)
%    Called with a pair of vectors of the same length, 
%    incircle computes the convex hull in 2-dimensions,
%    then finds the maximal radius in-circle that lies
%    entirely within the polygon of the convex hull.
%
% [C,R] = incircle(xy)
%    Called with a single nx2 array, incircle treats
%    this as a list of points in 2 dimensions. Each row
%    of xy is treated as a single point. The convex hull
%    is computed and the maximal in-circle then computed
%    for the convex hull.
%
% [C,R] = incircle(xy,edges)
%    Called with two arguments, the first of which is an
%    nx2 array of data points, and the second a list of
%    edges of the convex hull as would be generated by
%    convhull (or convhulln), the maximal in-circle is
%    generated for the hull as provided. No test is made
%    to assure the hull is truly convex.
%
% All input forms return a row vector of length 2 (C)
% which defines the center of the in-circle, and the
% radius (R) of the in-circle.
%
% Example:
% 
% [C,R] = incircle(rand(1000,1),rand(1000,1))
%
% C =
%      0.50191      0.49999
%
% R =
%      0.49788
%
% Example:
%
% x = randn(100,1) + 1;
% y = randn(100,1) - 2.5;
% H = convhull(x,y);
% [C,R] = incircle([x,y],H)
%
% C =
%      1.0212      -2.3458
%
% R =
%      2.1329
%
% t = linspace(0,2*pi,200);
% xc = R*cos(t) + C(1);
% yc = R*sin(t) + C(2);
% plot(x(H),y(H),'-r',x,y,'ko',xc,yc,'-b')
% axis equal
% xlabel X
% ylabel Y
% title 'Convex hull, with polygon in-circle'
% grid on
%
%
% See also: minboundcicle, convhull
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 6/14/09

if (nargin < 1) || (nargin > 2)
  error('INCIRCLE:improperarguments', ...
    'incircle requires exactly 1 or 2 arguments')
elseif (nargin == 2) && isvector(x) && isvector(y) && (numel(x) == numel(y))
  % a pair of vectors
  xy = [x(:),y(:)];
  % compute the hull. I prefer convhulln.
  edges = convhulln(xy);
elseif (nargin == 1) && (size(x,2) == 2)
  % a single list of points as rows of x
  xy = x;
  edges = convhulln(xy);
elseif (nargin == 2) && (size(x,2) == 2) && (size(y,2) == 2)
  % y must be a list of edges from convhulln
  xy = x;
  edges = y;
elseif (nargin == 2) && (size(x,2) == 2) && isvector(y) && (y(1) == y(end))
  % y must be a list of edges from convhull
  xy = x;
  I = y(:);
  % in case the first point was not wrapped in the polygon
  if I(1) ~= I(end)
    I = [I;I(1)];
  end
  edges = [I(1:(end-1)),I(2:end)];
else
  % none of the forms I allow seem to fit.
  % give up and throw an error.
  error('INCIRCLE:invaliddata', ...
    'x and y do not seem to fit into any of the allowed forms for incircle input')
end
ne = size(edges,1);

% the end points of each edge are...
A = xy(edges(:,1),:);
B = xy(edges(:,2),:);

% the normal vector to each edge
N = (B - A)*[0 1;-1 0];

% normalize to unit length
L = sqrt(sum(N.^2,2));
N = N./[L,L];

% a central point inside the hull itself
C0 = mean(A,1);

% ensure the normals point inwards. While
% I can use the signed area for a hull as
% generated by convhull (suggestion by Bruno)
% this test is also efficient, and it deals
% with the case of a hull provided from some
% unknown and less trusted source.
k = sum(N.*bsxfun(@minus,C0,A),2) < 0;
N(k,:) = -N(k,:);

% formulate the linear programming problem.
% given a point C inside the hull, the distance
% from C to the edge that contains A(i,:) is
% given by (dot(N(i,:),C - A(i,:))). If this
% number is negative for any edge of the hull,
% then C is outside of the hull.
%
% Now think of it as a set of slack variables,
% one for each edge of the hull. Given the
% unknown point C that defines the center of
% the in-circle,
%
%  dot(N(i,:),C-A(i,:)) - S(i) == 0
%
% Thus the vector S is of length ne, where ne is
% the number of edges in the convex hull. Every
% element of S must be positive.
%
% Create one more slack variable, a scalar R.
% R is the minimum value that S attains for a
% given point C. 
%
%   R >= 0
%   S(i) >= R
%
% The objective for our linear programming problem
% is simple. It is just -R. When we find the
% solution to the LP problem that satisfies the
% equality constraints above between C and S, the
% bound constraint on T, and the inequality
% constraints between S and R, the solution yields
% the maximal radius in-circle that fits inside
% the convex hull polygon.
%
% The unknowns for the LP are a list of ne+3
% variables. The first two variables represent
% the center coordinates of the circle. X(3) = R
% is the circle radius. The remainder of the
% variables are the slack variabls S(i).

% equality constraints defined by the dot products
% with the normal vectors.
Aeq = [N,zeros(ne,1),-eye(ne)];
beq = sum(N.*A,2);

% lower bounds only for the slack variables
LB = [-inf; -inf; 0; zeros(ne,1)];
% there are no upper bounds
UB = [];

% inequalities defined by the slack variable
% constraints
A = [zeros(ne,2),ones(ne,1),-eye(ne)];
b = zeros(ne,1);

% the objective is just -R
f = zeros(ne+3,1);
f(3) = -1;

% just turn off the output message
options = optimset('linprog');
options.Display = 'off';

% linprog does all of the hard work.
result = linprog(f,A,b,Aeq,beq,LB,UB,[],options);

% unpack the circle parameters
C = result(1:2)';
R = result(3);



