A common problem in handling geometric data is determining the center of a given polygon. This is not quite so easy as it sounds as there is not a single definition of center that makes sense in all cases. For instance, sometimes computing the center of a polygon's bounding box may be sufficient. In some instances this may give a point on an edge (consider a right triangle). If the given polygon is non-convex, that point may not even be inside or on its boundary.

This post looks at computing Chebyshev centers for arbitrary convex polygons. We employ essentially the same model as in Boyd & Vandenberghe's *Convex Optimization* text, but using Gurobi instead of CVXOPT.

Consider a polygon defined by the intersection of a finite number of half-spaces, $Au \le b$. We assume we are given the set of vertices, $V$, in clockwise order around the polygon. $E$ is the set of edges connecting these vertices. Each edge in $E$ defines a boundary of the half-space $a_i'u \le b_i$.

$E = \{((1,1),(2,5)), ((2,5),(5,4)), ((5,4),(6,2)), ((6,2),(4,1)), ((4,1),(1,1))\}$

The Chebyshev center of this polygon is the center point $(x, y)$ of the maximum radius inscribed circle. That is, if we can find the largest circle that will fit inside our polygon without going outside its boundary, its center is the point we are looking for. Our decision variables are the center $(x, y)$ and the maximum inscribed radius, $r$.

In order to do this, we consider the edges independently. The long line segment below shows an arbitrary edge, $a_i'u \le b_i$. The short line connected to it is orthogonal in the direction $a$. $(x, y)$ satisfies the inequality.

The shortest distance from $(x, y)$ will be in the direction of $a$. We'll call this distance $r$. If we were to move the edge so it had the same slope but went through $(x, y)$, its distance from $a_i'u = b_i$ would be $r||a_i||_2$. Thus we can add a constraint of the form $a_i'u + r||a_i||_2 \le b_i$ for each edge and maximize the value of $r$ as our objective function.

max $r$

subject to $(y_i-y_j)x + (x_j-x_i)y + r\sqrt{(x_j-x_i)^2 + (y_j-y_i)^2} \le (y_i-y_j)x_i + (x_j-x_i)y_i\, \forall\, ((x_i,y_i), (x_j,y_j)) \in E$

As this is linear, we can solve it using any LP solver. The following code does so with Gurobi.

#!/usr/bin/env python from gurobipy import Model, GRB from math import sqrt vertices = [(1,1), (2,5), (5,4), (6,2), (4,1)] edges = zip(vertices, vertices[1:] + [vertices[0]]) m = Model() r = m.addVar() x = m.addVar(lb=-GRB.INFINITY) y = m.addVar(lb=-GRB.INFINITY) m.update() for (x1, y1), (x2, y2) in edges: dx = x2 - x1 dy = y2 - y1 m.addConstr((dx*y - dy*x) + (r * sqrt(dx**2 + dy**2)) <= dx*y1 - dy*x1) m.setObjective(r, GRB.MAXIMIZE) m.optimize() print 'r = %.04f' % r.x print '(x, y) = (%.04f, %.04f)' % (x.x, y.x)

The model output shows our center and its maximum inscribed radius.

r = 1.7466 (x, y) = (3.2370, 2.7466)

Question for the reader: in certain circumstances, such as rectangles, the Chebyshev center is ambiguous. How might one get around this ambiguity?

Another definition of Chebyshev center is the center of the smallest enclosing circle. It yields a non ambiguous answer for rectangles.

ReplyDeleteActually, computing the center of the containing circle yields to a non convex MIQCP when stated naively. One could implement Nimrod 's algorithm in that case: Megiddo, Nimrod (1983), "Linear-time algorithms for linear programming in R3 and related problems", SIAM Journal on Computing 12 (4): 759–776, doi:10.1137/0212052, MR 721011.

ReplyDeleteI've taken another route here: https://www.ibm.com/developerworks/community/blogs/jfp/entry/chebyshev_centers?lang=en