# Find the largest area polygon built from a given list of vertices

This is code that gets a list of polygon vertices, non-ordered, and finds the order in which they should be arranged in order to create the polygon with the largest area.

There are two key elements in the approach:

1. Testing the area of the polygon is done by Monte Carlo calculation. Random points are dropped on the plane and we count the fraction that is inside the polygon. Note that the polygon is not necessarily concave.

2. Since going over all permutation of the vertices is too extensive, we employ a statistical trick. We have a greedy algorithm to unwind a polygon by identified pair of edges that cut each other and "untangle" them, then repeating until the polygon is completely untangled. Sometime untangling fails, but most of the time it's OK. Then we repeat shuffle-untangle cycle a low number of times, testing for the area, until we don't get a better candidate.

The code also plays a little with data structures, and have some plotting functions, etc.

I'm not sure I am too happy with the solution since the Monte Carlo findarea is slow and inaccurate. Also, the unwindpolygon untangling is not always successful and needs a bailout to break when no solution is found, which means that some tests are wasted. For large number of vertices it gets really slow.

import numpy as np
import matplotlib.pyplot as plt
import random
import copy

class point:
"""A point in 2d"""

def __init__(self, px, py):
self.x = float(px)
self.y = float(py)
self.type = "point"

def __repr__(self):
return "point(" + str(self.x) + ", " + str(self.y) + ")"

def rotate(self, theta):
"""Rotate the point around origin by angle theta"""
rotmat = np.array([[np.cos(theta), np.sin(theta)],
[-np.sin(theta), np.cos(theta)]])
self.x, self.y = np.matmul(rotmat, [self.x, self.y])

def plotpoints(points):
"""plots points"""
for k in range(len(points)):
plt.plot(points[k].x, points[k].y, '.')
plt.text(points[k].x, points[k].y, k)

class line:
"""A line in 2d"""

def __init__(self, p1, p2):
assert type(p1) is point and type(p2) is point,\
"Input must be points"
self.p1 = p1
self.p2 = p2
self.type = "line"

def __repr__(self):
return "line(" + str(self.p1) + ", " + str(self.p2) + ")"

def length(self):
return np.sqrt((self.p1.y - self.p2.y)**2 +
(self.p1.x - self.p2.x)**2)

def rotate(self, theta):
"""Rotate the line around origin by angle theta"""
self.p1.rotate(theta)
self.p2.rotate(theta)

def plotlines(lines):
"""Plots line"""
for k in range(len(lines)):
plt.plot([lines[k].p1.x, lines[k].p2.x],
[lines[k].p1.y, lines[k].p2.y], '-')

class polygon:
"""A polygon in 2d"""

def __init__(self, points):
self.n = len(points)
for p in points:
pass
assert type(p) is point,\
"Edges must be lines"
self.vertices = points
self.edges = []
for k in range(self.n - 1):
self.edges.append(line(points[k], points[k+1]))
self.edges.append(line(points[self.n - 1], points[0]))
self.type = polygon

def __repr__(self):
ms = [p for p in self.vertices]
return "polygon(" + str(ms) + ")"

def rotate(self, th):
"""rotate a polygon"""
for k in range(self.n):
self.vertices[k].rotate(th)
self.edges[k].rotate(th)

def plotpolygon(self):
"""Plots a polygon"""
plotpoints(self.vertices)
plotlines(self.edges)

def rotate(obj, th):
"""Rotate an object obj by angle th"""
obj.rotate(th)

def isininterval(x0, x1, x2):
"""Numerically test if x0 is in the open interval x1 to x2"""
eps = 1e-5
if x1+eps < x0 < x2-eps or x2+eps < x0 < x1-eps:
return True
else:
return False

def isintersect(line1, line2):
"""Test if two line segments l1 and l2 intersect.
Lines are considered open so intersect at line edge
counts as false."""
l1 = copy.deepcopy(line1)
l2 = copy.deepcopy(line2)
if l1.length() == 0 or l2.length() == 0:
return False
while (l1.p1.x == l1.p2.x) or (l2.p1.x == l2.p2.x):
# Rotate the whole world by 45deg to avoid dealing
# with infinite a's which comes as a result of a
# line being vertical
rotate(l1, 45*np.pi/180)
rotate(l2, 45*np.pi/180)
# Solve for the line formula for both segments
a1 = (l1.p2.y - l1.p1.y) / (l1.p2.x - l1.p1.x)
a2 = (l2.p2.y - l2.p1.y) / (l2.p2.x - l2.p1.x)
b1 = l1.p1.y - a1 * l1.p1.x
b2 = l2.p1.y - a2 * l2.p1.x
if a1 == a2:
if b1 == b2:
# This is a special case where two line segments
# are on the same line
if (isininterval(l2.p1.x, l1.p1.x, l1.p2.x) or
isininterval(l2.p2.x, l1.p1.x, l1.p2.x) or
isininterval(l1.p1.x, l2.p1.x, l2.p2.x) or
isininterval(l1.p2.x, l2.p1.x, l2.p2.x)):
return True
else:
return False
else:
return False
else:
# The intersection's x
x = - (b2 - b1) / (a2 - a1)
# If the intersection x is within the interval of each segment
# then there is an intersection
if (isininterval(x, l1.p1.x, l1.p2.x) and
isininterval(x, l2.p1.x, l2.p2.x)):
return True
else:
return False

def anyintersect(p):
"""Test if there is any intersect left at the polygon p
or is it completly untangled"""
assert type(p) is polygon, "Input must be polygon."
return any([isintersect(p.edges[j], p.edges[k])
for j in range(p.n)
for k in range(p.n)])

def unwindpolygon(p, bailout = 50):
"""Greadily untangle a polygon by un crossing intersecting
edges"""
t = 0
while anyintersect(p) and t < bailout:
t += 1
breakingflag = False
for j in range(p.n):
for k in range(p.n):
if isintersect(p.edges[j], p.edges[k]):
if j < p.n - 1:
p12 = p.vertices[j+1]
else:
p12 = p.vertices[0]
if k < p.n - 1:
p22 = p.vertices[k+1]
p.vertices[k+1] = p12
else:
p22 = p.vertices[0]
p.vertices[0] = p12
if j < p.n - 1:
p.vertices[j+1] = p22
else:
p.vertices[0] = p22
breakingflag = True
break
if breakingflag:
break
p.edges = polygon(p.vertices).edges

def shufflepolygon(p):
"""Randomely shuffle polygon vertices"""
random.shuffle(p.vertices)
p.edges = polygon(p.vertices).edges

# not used. Good for testing be rearranging manually a polygon
def rearrangepolygon(p, per):
"""p is polygon. per is a permutation of indices)"""
p.vertices = [p.vertices[k] for k in per]
p.edges = polygon(p.vertices).edges

def findbbox(pol):
"""Find the bounding box of a polygon"""
maxx = minx = pol.vertices[0].x
maxy = miny = pol.vertices[0].y
for k in range(pol.n):
minx = np.min([minx, pol.vertices[k].x])
maxx = np.max([maxx, pol.vertices[k].x])
miny = np.min([miny, pol.vertices[k].y])
maxy = np.max([maxy, pol.vertices[k].y])
return [minx, maxx, miny, maxy]

def isinside(pt, pol):
"""Test whether a point pt is inside polygon pol. The test
is done by taking a line from pt to the edges of the bounding box
and test if the number of intersections with the polygon edges is
odd (pt is inside) or even (pt is outside)"""
_, mx, _, _ = findbbox(pol)
l = line(pt, point(mx + 1, pt.y + 1))
# line is not completly horizontal to avoid specuial cases in "isintersect"
res = sum([isintersect(l, pol.edges[k]) for k in range(pol.n)])
if not res % 2:
return False
else:
return True

def findarea(p, n=1000):
"""Finds the area by monte-carlo. Drop points and count the fraction
that is inside the polygon"""
minx, maxx, miny, maxy = findbbox(p)
area = (maxx - minx) * (maxy - miny)
fra = 0
for _ in range(n):
ptx = minx + (maxx - minx)*np.random.rand()
pty = miny + (maxy - miny)*np.random.rand()
pt = point(ptx, pty)
fra += isinside(pt, p)
return area * fra / n

def randomtestpolygon(p, t=20, acc = 200):
"""Repeatedly shuffle and untangle, keeping record of the largest
area result, for t trials. Some plottings are added to
show what's going on."""
hasbeenbestfor = 0
bestscore = findarea(p, acc)
pbest = copy.deepcopy(p)
plt.figure()
plt.subplot(121)
plotpolygon(p)
plt.axis('square')
plt.title("best area:" + str(bestscore))
print("trial, largest area, current area")
# Look for the configuration that can hold against competition
# for t trials
while hasbeenbestfor < t:
shufflepolygon(p)
unwindpolygon(p)
currentscore = findarea(p, acc)
print(hasbeenbestfor, bestscore, currentscore)
plt.subplot(122)
plt.cla()
plotpolygon(p)
plt.axis('square')
plt.title("current area:" + str(currentscore))
plt.draw()
plt.pause(0.1)
if currentscore > bestscore:
bestscore = currentscore
pbest = copy.deepcopy(p)
hasbeenbestfor = 0
plt.subplot(121)
plt.cla()
plotpolygon(p)
plt.axis('square')
plt.title("best area:" + str(bestscore))
else:
hasbeenbestfor += 1
p.vertices = pbest.vertices
p.edges = pbest.edges

if __name__ == '__main__':
plt.close('all')
n = 7  # number of vertices
# This is just random vertices generator. The details are not important
# This way ensures that a nice untangled polygon is likely
r = 10+100*np.random.rand(n,)
th = [(k + 10*n/360*np.random.randn())*np.pi/180
for k in np.linspace(0, 360, n)]
vertices = [point(r_*np.cos(th_), r_*np.sin(th_))
for r_, th_ in zip(r, th)]
p = polygon(vertices)
shufflepolygon(p)
randomtestpolygon(p, 30, 1000)

• All of the vertices need to be used, right? – 200_success Jul 5 '16 at 17:39
• @200_success - Correct. – Aguy Jul 5 '16 at 17:51
• One thing I should definitely do (post posting thoughts) is calculate the bounding box only once and not for every pt. Have it as a property of the polygon object. – Aguy Jul 5 '16 at 19:00

Your algorithm seems to rely on randomness more than necessary.

One intuitively obvious conclusion is that the polygon with the largest area should be a simple polygon (i.e., no self-intersecting edges). You have implemented that, which is a good start.

It it easy to calculate the area of any simple polygon exactly using the shoelace formula. The algorithm is $O(\left|V\right|)$. No Monte-Carlo simulation is needed.

The convex hull of the points should provide a skeleton of the solution. The question is, what to do with the interior points? Intuitively, I think that the solution should involve computing the convex hull of the remaining interior points, then interspersing the traversal of the outer polygon with the traversal of the inner polygon in such a way that the area is maximized. Then, recurse as necessary if the inner polygon also has interior points.

I can think if 2 approaches to solve this:

1. Brute force approach : Calculate the greatest area polygon starting from a triangle by applying shoelace formula on all possible combination of points. A python implementation is given below.

2. Remove all the internal points which are not contributing to the convex hull of the points: this can be done by arranging all the points in a sorted manner and then deciding for each point, this decision can be taken on the basis of finding a triangle encompassing the point.

import math
import itertools

def clockwiseangle_and_distance(point):
vector = [point[0]-origin[0], point[1]-origin[1]]
lenvector = math.hypot(vector[0], vector[1])
if lenvector == 0:
return -math.pi, 0
normalized = [vector[0]/lenvector, vector[1]/lenvector]
dotprod  = normalized[0]*refvec[0] + normalized[1]*refvec[1]
diffprod = refvec[1]*normalized[0] - refvec[0]*normalized[1]
angle = math.atan2(diffprod, dotprod)
if angle < 0:
return 2*math.pi+angle, lenvector
return angle, lenvector

# shoelace formula implementation
def PolygonArea(corners):
n = len(corners)
area = 0.0
for i in range(n):
j = (i + 1) % n
area += corners[i][0] * corners[j][1]
area -= corners[j][0] * corners[i][1]
area = abs(area) / 2.0
return area

# idea is to find the area of biggest polygon in the given set of points

pts = [(10, 10), (0, 0), (0, 10), (5,5), (10,0)]
#pts = [(6, 39), (28, 25), (28,13),(31,3),(11,19),(31,17),(26,19),(18,13),(30,11),(25,20)]

origin = pts[0]
refvec = [0, 1]
init_set_size = 3
max_set_size = len(pts)+1
#print(max_set_size)
area_list = []

for i in range(init_set_size, max_set_size, 1):
#print(i)
comb = list(itertools.combinations(pts, i))
for n in comb:
corners = sorted(n, key=clockwiseangle_and_distance)
area_list.append(PolygonArea(corners))

area_list.sort(reverse=True)
print(area_list[0])