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"""

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:
            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):

def plotpolygon(self):
    """Plots a polygon"""

def rotate(obj, th):
    """Rotate an object obj by angle 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
        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
                return False
            return False
        # 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
            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
    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]
                        p12 = p.vertices[0]
                    if k < p.n - 1:
                        p22 = p.vertices[k+1]
                        p.vertices[k+1] = p12
                        p22 = p.vertices[0]
                        p.vertices[0] = p12
                    if j < p.n - 1:
                        p.vertices[j+1] = p22
                        p.vertices[0] = p22
                    breakingflag = True
            if breakingflag:
        p.edges = polygon(p.vertices).edges

def shufflepolygon(p):
    """Randomely shuffle polygon 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
        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.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:
        currentscore = findarea(p, acc)
        print(hasbeenbestfor, bestscore, currentscore)
        plt.title("current area:" + str(currentscore))
        if currentscore > bestscore:
            bestscore = currentscore
            pbest = copy.deepcopy(p)
            hasbeenbestfor = 0
            plt.title("best area:" + str(bestscore))
            hasbeenbestfor += 1
    p.vertices = pbest.vertices
    p.edges = pbest.edges

if __name__ == '__main__':
    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)
    # but let's start with something hard
    randomtestpolygon(p, 30, 1000)
  • \$\begingroup\$ All of the vertices need to be used, right? \$\endgroup\$ Commented Jul 5, 2016 at 17:39
  • \$\begingroup\$ @200_success - Correct. \$\endgroup\$
    – Aguy
    Commented Jul 5, 2016 at 17:51
  • \$\begingroup\$ 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. \$\endgroup\$
    – Aguy
    Commented Jul 5, 2016 at 19:00

2 Answers 2


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
area_list = []

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


Refs: Shoelace formula implementation Arranging coordinates in a clockwise manner


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.