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.
% [C,R] = incircle(rand(1000,1),rand(1000,1))
% C =
% 0.50191 0.49999
% R =
% 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 =
% t = linspace(0,2*pi,200);
% xc = R*cos(t) + C(1);
% yc = R*sin(t) + C(2);
% 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: firstname.lastname@example.org
% Release: 1.0
% Release date: 6/14/09
if (nargin < 1) || (nargin > 2)
'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)];
edges = [I(1:(end-1)),I(2:end)];
% none of the forms I allow seem to fit.
% give up and throw an error.
'x and y do not seem to fit into any of the allowed forms for incircle input')
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
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);