Digesting a polygon is a process in which you take an input polygon, and create another polygon that can contain the entirety of the input polygon, while having one side less than it.

The goal of this challenge is to generate the optimal polygon digest- a polygon digest that has the minimum area possible.

Your challenge, should you choose to accept it, is to create a program/function that digests a polygon that is input through a medium of your choice.

The program will accept an input 2D polygon in the following format:


Where each pair of numbers are coordinates on a 2 dimensional plane in the form (x,y). Each pair is a point on the polygon. The points are given in clockwise order. You can assume valid input. The program should accept decimal numbers too.

Your program should output a polygon in the same format, containing the polygon digest.

The requirements of the polygon digest are:

  1. Completely contains the input polygon
  2. Has the least possible area
  3. Has 1 side less than the input polygon

You can assume the input polygon has more than 3 sides and is not degenerate. You must handle both convex and concave polygons.

Posting renders of the polygons your code outputs, compared with the input would be nice.

Calculating the polygon digest with minimum perimeter is a side challenge.

Test cases: (no solutions yet...)

  • \$\begingroup\$ Are there known results on optimal polygon digests? Is there a fixed constant c for which all but c of the new polygon's sides overlap the old polygon's sides? \$\endgroup\$
    – xnor
    Nov 15, 2014 at 4:22
  • \$\begingroup\$ @xnor I do not know of any such algorithm that can calculate the optimal polygon digest, should i make this a code challenge then? I believe that the solution would have all sides the same as the original polygon, except for 2. But maybe i am wrong. \$\endgroup\$
    – vero
    Nov 15, 2014 at 4:25
  • \$\begingroup\$ Two is too low, consider a square turning into a triangle. Maybe 3? \$\endgroup\$
    – xnor
    Nov 15, 2014 at 4:26
  • \$\begingroup\$ @xnor I think so. I was thinking about the concave case and polygons with larger numbers of sides, my bad. \$\endgroup\$
    – vero
    Nov 15, 2014 at 4:29
  • \$\begingroup\$ An algorithm that i think could output the optimal result would be to iterate through every side of the polygon, 'un-truncating' it (making the two adjacent sides meet at a corner in a way that the adjacent sides retain their slope), and checking the new area of the polygon. Then select the one with the least area. (Again this would not work for a perfect square, haha...) \$\endgroup\$
    – vero
    Nov 15, 2014 at 4:31

1 Answer 1



The code is somewhat messy, but it basically tries three methods on every point in the polygon and keeps whichever one is best. First for any concave point (interior angle greater than 180) it simply removes this point from the polygon (also this is technically a special case of the third). Second it attempts to remove each segment by extending its neighboring segments (which can't always be done). Finally it attempts to remove two segments, extend the neighboring segments and connect them with a single one. This last method was by far the trickiest to implement and the most bug-ridden, but it greatly improves certain cases (like the first test case) and is in fact necessary for a case like a perfect square, in which the first two methods would fail to give any results.

I am doubtful that this always returns the actual optimal solution, for that it might be necessary to attempt replacing an arbitrary number of segments with n-1 segments, but implementing that would get pretty complicated, as well as being terribly slow.

import math

## Copied from http://stackoverflow.com/questions/24467972/calculate-area-of-polygon-given-x-y-coordinates

def PolygonArea(corners):
    n = len(corners) # of 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][2]
    if area < 0:
        return None  # Don't give ccw polygons
    area = abs(area) / 2.0
    return area

## Copied from http://stackoverflow.com/questions/20677795/find-the-point-of-intersecting-lines

def line_intersection(line1, line2):
    xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
    ydiff = (line1[0][3] - line1[1][4], line2[0][5] - line2[1][6])

    def det(a, b):
        return a[0] * b[1] - a[1] * b[0]

    div = det(xdiff, ydiff)
    if div == 0:
        raise Exception('lines do not intersect')

    d = (det(*line1), det(*line2))
    x = det(d, xdiff) / div
    y = det(d, ydiff) / div
    return x, y

def dist((x0,y0),(x1,y1)):
    return math.hypot(x0-x1,y0-y1)
def angle(x1, y1, x2, y2):
    # Use dotproduct to find angle between vectors
    # This always returns an angle between 0, pi
    numer = (x1 * x2 + y1 * y2)
    denom = math.sqrt((x1 ** 2 + y1 ** 2) * (x2 ** 2 + y2 ** 2))
    return math.acos(numer / denom) 

def cross_sign(x1, y1, x2, y2):
    # True if cross is positive
    # False if negative or zero
    return x1 * y2 > x2 * y1

def ccw(A,B,C):
    return (C[1]-A[1])*(B[0]-A[0]) > (B[1]-A[1])*(C[0]-A[0])
def segment_intersect(A,B,C,D):
    return ccw(A,C,D) != ccw(B,C,D) and ccw(A,B,C) != ccw(A,B,D)
## Set stuff up

raw_points = input()
points = []
for i in range(0,len(raw_points),2):

AREA = PolygonArea(points)
angles = []
for i in range(len(points)):
    p1 = points[i-1]
    ref = points[i]
    p2 = points[(i+1)%len(points)]
    x1, y1 = p1[0] - ref[0], p1[1] - ref[1]
    x2, y2 = p2[0] - ref[0], p2[1] - ref[1]

    ang = angle(x1, y1, x2, y2)
    if cross_sign(x1, y1, x2, y2):
        ang = -ang
best_area = ('','')

## Try removing any concave points

for i in range(len(points)):
    if angles[i] > math.pi:
        new = points[:i]+points[i+1:]
        area = PolygonArea(new)
        if area < best_area[0]:
            best_area = (area,new)
## Try removing a single segment by extending neighboring segments

for i in range(len(points)):
    if angles[i] < math.pi and angles[(i+1)%len(points)] < math.pi and angles[i]+angles[(i+1)%len(points)] > math.pi:
        intersect = line_intersection((points[i-1],points[i]),(points[(i+1)%len(points)],points[(i+2)%len(points)]))
        if i+1 == len(points):
            new = points[1:i]+[intersect]
            new = points[:i]+[intersect]+points[i+2:]
        area = PolygonArea(new)
        if area is not None and area < best_area[0]:
            best_area = (area,new)

## Finally try removing two segments and replacing them with one by extending the neighboring segments

def extended(i,length):
    p0 = points[i-1]
    p1 = points[i]
    ang = math.atan2(p1[1]-p0[1],p1[0]-p0[0])
    return p0[0] + math.cos(ang)*length, p0[1] + math.sin(ang)*length
def get_extended_points(i,length):
    if angles[i-1] > math.pi or angles[(i+1)%len(points)] > math.pi:
        return None

    pe0 = extended(i-1,length)
    p1 = points[i]; p2 = points[(i+1)%len(points)]
    new_angle = (math.atan2(p2[1]-p1[1],p2[0]-p1[0]) - math.atan2(pe0[1]-p1[1],pe0[0]-p1[0])) % (2*math.pi)
    if new_angle + angles[(i+1)%len(points)] > math.pi:
            pe1 = line_intersection((pe0,points[i]),(points[(i+1)%len(points)],points[(i+2)%len(points)]))
            return None
        for j in range(len(points)):
            t0 = points[j-1]; t1 = points[j]
            if segment_intersect(t0,t1,pe0,points[i-1]) or segment_intersect(t0,t1,pe1,points[(i+1)%len(points)]):
                return None
        if dist(pe0,points[i-2]) < dist(points[i-1],points[i-2]) or dist(pe1,points[(i+2)%len(points)]) < dist(points[(i+1)%len(points)],points[(i+2)%len(points)]):
            return None
        if i+1 == len(points):
            new = [pe1] + points[1:i-1] + [pe0]
        elif i == 0:
            new = [pe1] + points[2:i-1] + [pe0]
            new = points[:i-1] + [pe0, pe1] + points[i+2:]
        return new
    return None
for i in range(len(points)):
    bounds = [0.0,100.0]
    ITER = 3
    while ITER:
        ITER -= 1
        result_list = []
        for e in range(100):
            extend = e*(bounds[1]-bounds[0])/100. + bounds[0]
            new = get_extended_points(i,extend)
            if new == None:
            area = PolygonArea(new)
            if area == None or area < AREA:
            if area < best_area[0]:
                best_area = (area,new)
        last = ''
        direction = 0
        for e,result in enumerate(result_list):
            if result == '' or last == '':
                last = result
            if result < last:
                direction = -1
            elif direction == -1:
                low = (e-2)*(bounds[1]-bounds[0])/100. + bounds[0]
                high = e*(bounds[1]-bounds[0])/100. + bounds[0]
                bounds = [low,high]
            elif direction == 0:
                low = (e-1)*(bounds[1]-bounds[0])/100. + bounds[0]
                high = e*(bounds[1]-bounds[0])/100. + bounds[0]
                bounds = [low,high]
            low = e*(bounds[1]-bounds[0])/100. + bounds[0]
            high = bounds[1]
            bounds = [low,high]

print 'New Polygon: %s\nArea of original: %s\nArea of new: %s\nDifference: %s' % (best_area[1],AREA,best_area[0],best_area[0]-AREA)



New Polygon: [(23.14, -20.795), (13.509300122415228, 42.47623848849403), (-45.19115435248799, -12.367823474893163)]
Area of original: 1364.81275
Area of new: 2121.11857765
Difference: 756.305827651

enter image description here


New Polygon: [(21.9, -20.795), (16.775, 21.005), (-15.5725, 15.305), (-19.6725, -15.515)]
Area of original: 894.99775
Area of new: 1342.125225
Difference: 447.127475

enter image description here


New Polygon: [(23.14, -20.795), (16.7775, 21.005), (-27.4, 15.305), (-19.6725, -15.515), (-43.3, -30.3), (1.0767001474878615, -37.91779866292026), (-5.8, 1.1)]
Area of original: 1886.22665625
Area of new: 1934.3095735
Difference: 48.0829172453

enter image description here


New Polygon: [(28.57, -20.795), (17.8, 21.005), (-33.042021159455786, 12.321217825787757), (-19.6725, -47.3), (-1.2, -2.8)]
Area of original: 1643.63435625
Area of new: 1896.08947014
Difference: 252.455113887

enter image description here


New Polygon: [(1.0000000000000004, -2.0), (6.999999999999999, 4.0), (1, 4)]
Area of original: 9.0
Area of new: 18.0
Difference: 9.0

*The example in the question was not the same clockwise-ness as the rest

enter image description here

  • \$\begingroup\$ Wow... I'm quite surprised. I'll give the question some more time before I mark you as the accepted answer. I'm still trying to prepare a random polygon generation tool... (got carried away making it perfectly random, by generating a random set of points and then rejecting the ones that are invalid, multiple times). +1 \$\endgroup\$
    – vero
    Nov 17, 2014 at 23:54
  • \$\begingroup\$ @Rodolvertice Yeah this was a fun one to try and do. Before I began writing anything I did try and solve a couple on paper but I honestly can't figure out any way you could prove any results except maybe for some very trivial specific cases. I'm very interested to see if anyone comes up with better solutions than mine. \$\endgroup\$
    – KSab
    Nov 18, 2014 at 0:07

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