Re: the smallest circle that can enclosed a convex hull

"Bruno Luong" <b.luong@xxxxxxxxxxxxxxxxxxxx> wrote in message <h133ml$23h$1@xxxxxxxxxxxxxxxxxx>...
For the largest inner circle problem, I would suggest to use implicit equation of the circle:

f(x,y) = x^2 + y^2 + p*x + q*y + rho = 0 (C)

The convex set is
X = (x,y)'
A*X <= B
A is (m x 2) and B is (m x 1), A has minimal number of rows m.

Call the "extended" hull is H = {X : A*X = B}
Circle is inside if and only if
f(H) >= 0.

The square of the radius is rho (p^2+q^2)/4 - rho (could be mistake here) is to be maximize.
This is constrained quadratic programming. There is unique solution, and plenty of numerical methods to solve it.

Now back to your keyboard. Have fun.

Bruno

But, it would be such a pity to solve a QP when
an LP would do nicely.

[C,R] = incircle(rand(1000,2))
C =
0.50009 0.4977
R =
0.49719

(I'll post this on the FEX tomorrow after I finish
cleaning it up.)

John


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, 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 which
% defines the center of the in-circle, and the radius R
% of the in-circle.
%
%


if (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)



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...
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;

options = optimset('linprog');
options.Display = 'off';
[result,fval,exitflag] = linprog(f,A,b,Aeq,beq,LB,UB,[],options);

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

Relevant Pages

  • Re: find polygon
    ... perimeter of convex hull but would have some notches that intruded ... In most cases the perimiter of the convex hull isn't a triangle, ... the convex hull and doesn't have any straight edges, ...
    (comp.programming)
  • Another Difficult Combinatorial Geometry Problem
    ... in the plane is decomposable into a union of 3 or fewer convex sets. ... local non-convexity lie in the kernel of the set and in the planar, ... planar 3-convex sets. ... convex hull is a triangle or pentagon are relatively straightforward ...
    (sci.math.research)
  • Generalized convex perimeter of an n-d object
    ... One could define a "convex perimeter" of any bounded set of points S ... in the plane as p_2= the perimeter of the convex hull of S. ... Note that this generalized convex perimeter is a sort of dual to what ...
    (sci.math)
  • Re: convex hull
    ... > hull, then will the convex hull, defined in this manner, of S be convex ... What if S is the three-element set consiting of the vertices of a triangle? ...
    (sci.math)