Digesting Polygons [closed]

Question

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

Answer 1

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