# Digesting Polygons

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:

1,1,1,4,4,4,4,1


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

23.14,-20.795,16.7775,21.005,-15.5725,15.305,-19.6725,-15.515
21.9,-20.795,16.775,21.005,-15.5725,15.305,11.5,0.7,-19.6725,-15.515
23.14,-20.795,16.7775,21.005,-27.4,15.305,-19.6725,-15.515,-43.3,-30.3,-6.6,-36.6,-1.2,-25,-5.8,1.1
28.57,-20.795,17.8,21.005,-15.5725,15.305,-26.8,-15.515,-19.6725,-47.3,-1.2,-2.8
1,1,1,4,4,4,4,1

• 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? – xnor Nov 15 '14 at 4:22
• @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. – rodolphito Nov 15 '14 at 4:25
• Two is too low, consider a square turning into a triangle. Maybe 3? – xnor Nov 15 '14 at 4:26
• @xnor I think so. I was thinking about the concave case and polygons with larger numbers of sides, my bad. – rodolphito Nov 15 '14 at 4:29
• 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...) – rodolphito Nov 15 '14 at 4:31

# Python

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):
points.append(raw_points[i:i+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
angles.append(ang%(2*math.pi))

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]
else:
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:
try:
pe1 = line_intersection((pe0,points[i]),(points[(i+1)%len(points)],points[(i+2)%len(points)]))
except:
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]
else:
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:
result_list.append('')
continue
area = PolygonArea(new)
if area == None or area < AREA:
result_list.append('')
continue

if area < best_area[0]:
best_area = (area,new)
result_list.append(area)

last = ''
direction = 0
for e,result in enumerate(result_list):
if result == '' or last == '':
last = result
continue

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]
break
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]
break
else:
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)


## Results:

23.14,-20.795,16.7775,21.005,-15.5725,15.305,-19.6725,-15.515

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


21.9,-20.795,16.775,21.005,-15.5725,15.305,11.5,0.7,-19.6725,-15.515

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


23.14,-20.795,16.7775,21.005,-27.4,15.305,-19.6725,-15.515,-43.3,-30.3,-6.6,-36.6,-1.2,-25,-5.8,1.1

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


28.57,-20.795,17.8,21.005,-15.5725,15.305,-26.8,-15.515,-19.6725,-47.3,-1.2,-2.8

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


1,1,4,1,4,4,1,4*

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

• 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 – rodolphito Nov 17 '14 at 23:54
• @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. – KSab Nov 18 '14 at 0:07