Polygon2 / Polygon /

Full commit
# -*- coding: utf-8 -*-

from cPolygon import Polygon
from math import sqrt, fabs, floor
from operator import add
from collections import defaultdict

def fillHoles(poly):
    Returns the polygon p without any holes.

        - p: Polygon
        new Polygon
    n = Polygon()
    [n.addContour(poly[i]) for i in range(len(poly)) if poly.isSolid(i)]
    return n

def pointList(poly, withHoles=1):
    Returns a list of all points of p.

        - p: Polygon
        list of points
    if not withHoles:
        poly = fillHoles(poly)
    return reduce(add, [list(c) for c in poly])

__left = lambda p: (p[1][0]*p[2][1]+p[0][0]*p[1][1]+p[2][0]*p[0][1]-
                    p[1][0]*p[0][1]-p[2][0]*p[1][1]-p[0][0]*p[2][1] >= 0)
def convexHull(poly):
    Returns a polygon which is the convex hull of p.

        - p: Polygon
        new Polygon
    points = list(pointList(poly, 0))
    u = [points[0], points[1]]
    for p in points[2:]:
        while len(u) > 2 and __left(u[-3:]):
            del u[-2]
    l = [points[0], points[1]]
    for p in points[2:]:
        while len(l) > 2 and __left(l[-3:]):
            del l[-2]
    return Polygon(u+l[1:-1])

def tile(poly, x=[], y=[], bb=None):
    Returns a list of polygons which are tiles of p splitted at the border values 
    specified in x and y. If you already know the bounding box of p, you may give 
    it as argument bb (4-tuple) to speed up the calculation.

        - p: Polygon
        - x: list of floats
        - y: list of floats
        - optional bb: tuple of 4 floats
        list of new Polygons
    if not (x or y):
        return [poly] # nothin' to do
    bb = bb or poly.boundingBox()
    x = [bb[0]] + [i for i in x if bb[0] < i < bb[1]] + [bb[1]]
    y = [bb[2]] + [j for j in y if bb[2] < j < bb[3]] + [bb[3]]
    cutpoly = []
    for i in range(len(x)-1):
        for j in range(len(y)-1):
    tmp = [c & poly for c in cutpoly]
    return [p for p in tmp if p]

def tileEqual(poly, nx=1, ny=1, bb=None):
    works like tile(), but splits into nx and ny parts.

        - p: Polygon
        - nx: integer
        - ny: integer
        - optional bb: tuple of 4 floats
        list of new Polygons
    bb = bb or poly.boundingBox()
    s0, s1 = bb[0], bb[2]
    a0, a1 = (bb[1]-bb[0])/nx, (bb[3]-bb[2])/ny 
    return tile(poly, [s0+a0*i for i in range(1, nx)], [s1+a1*i for i in range(1, ny)], bb)

def warpToOrigin(poly):
    Shifts lower left corner of the bounding box to origin.

        - p: Polygon
    x0, x1, y0, y1 = poly.boundingBox()
    poly.shift(-x0, -y0)

def centerAroundOrigin(poly):
    Shifts the center of the bounding box to origin.

        - p: Polygon
    x0, x1, y0, y1 = poly.boundingBox()
    poly.shift(-0.5*(x0+x1), -0.5*(y0+y1))

__vImp = lambda p: ((sqrt((p[1][0]-p[0][0])**2 + (p[1][1]-p[0][1])**2) +
                     sqrt((p[2][0]-p[1][0])**2 + (p[2][1]-p[1][1])**2)) *
def reducePoints(cont, n):
    Remove points of the contour 'cont', return a new contour with 'n' points.
    *Very simple* approach to reduce the number of points of a contour. Each point 
    gets an associated 'value of importance' which is the product of the lengths 
    and absolute angle of the left and right vertex. The points are sorted by this 
    value and the n most important points are returned. This method may give 
    *very* bad results for some contours like symmetric figures. It may even 
    produce self-intersecting contours which are not valid to process with 
    this module.

        - contour: list of points
        new list of points
    if n >= len(cont):
        return cont
    cont = list(cont)
    cont.insert(0, cont[-1])
    a = [(__vImp(cont[i-1:i+2]), i) for i in range(1, len(cont)-1)]
    ind = [x[1] for x in a[len(cont)-n-2:]]
    return [cont[i] for i in ind]

__linVal = lambda p: (p[1][0]-p[0][0])*(p[2][1]-p[0][1])-(p[1][1]-p[0][1])*(p[2][0]-p[0][0])
def prunePoints(poly):
    Returns a new Polygon which has exactly the same shape as p, but unneeded 
    points are removed. The new Polygon has no double points or points that are 
    exactly on a straight line.

        - p: Polygon
        new Polygon
    np = Polygon()
    for x in range(len(poly)): # loop over contours
        c = list(poly[x])
        c.insert(0, c[-1])
        # remove double points
        i = 1
        while (i < (len(c))):
            if c[i] == c[i-1]:
                del c[i]
                i += 1
        # remove points that are on a straight line
        n = []
        for i in range(1, len(c)-1):
            if __linVal(c[i-1:i+2]) != 0.0:
        if len(n) > 2:
            np.addContour(n, poly.isHole(x))
    return np

def cloneGrid(poly, con, xl, yl, xstep, ystep):
    Create a single new polygon with contours that are made from contour con from 
    polygon poly arranged in a xl-yl-grid with spacing xstep and ystep.

        - poly: Polygon
        - con: integer
        - xl: integer
        - yl: integer
        - xstep: float
        - ystep: float
        new Polygon
    p = Polygon(poly[con])
    for xi in range(xl):
        for yi in range(yl):
            p.cloneContour(0, xi*xstep, yi*ystep)
    return p

# following functions are contributed by Josiah Carlson <>
# the tileBSP() function is much faster for large polygons and a large number 
# of tiles than the original tile()
# background information can be found at:
# the original code is located here:

def _find_split(count_dict):
    When provided a dictionary of counts for the number of points inside each
    of the unit grid rows/columns, this function will return the best column
    choice so as to come closest to cutting the points in half.  It will also
    return the score, lower being better.

    Returns: (cutoff, score)
    # find the prefix sums
    tmp = {}
    for i in xrange(min(count_dict), max(count_dict)+1):
        tmp[i] = tmp.get(i-1, 0) + count_dict.get(i, 0)
    by_index = sorted(tmp.items())
    # the target number of points
    midpoint = by_index[-1][1] // 2
    # calculate how far off from the target number each choice would be
    totals = []
    for i in xrange(1, len(by_index)):
        totals.append(abs(by_index[i-1][1] - midpoint))
    # choose the best target number
    mi = min(totals)
    index = totals.index(mi)
    return by_index[index+1][0], totals[index]

def _single_poly(polygon, dim, maxv):
    for poly in polygon:
        if max(pt[dim] for pt in poly) > maxv:
            return False
    return True

def tileBSP(p):
    This generator function returns tiles of a polygon. It will be much 
    more efficient for larger polygons and a large number of tiles than the 
    original tile() function. For a discussion see:

        - p: Polygon
        tiles of the Polygon p on the integer grid
    _int = int
    _floor = floor

    work = [p]
    while work:
        # we'll use an explicit stack to ensure that degenerate polygons don't
        # blow the system recursion limit
        polygon = work.pop()

        # find out how many points are in each row/column of the grid
        xs = defaultdict(_int)
        ys = defaultdict(_int)
        for poly in polygon:
            for x,y in poly:
                xs[_int(_floor(x))] += 1
                ys[_int(_floor(y))] += 1

        # handle empty polygons gracefully
        if not xs:

        # handle top and right-edge border points
        mvx = max(max(x for x,y in poly) for poly in polygon)
        vx = _int(_floor(mvx))
        if len(xs) > 1 and mvx == vx:
            xs[vx-1] += xs.pop(vx, 0)
        mvy = max(max(y for x,y in poly) for poly in polygon)
        vy = _int(_floor(mvy))
        if len(ys) > 1 and mvy == vy:
            ys[vy-1] += ys.pop(vy, 0)

        # we've got a single grid, yield it
        if len(xs) == len(ys) == 1:
            yield polygon

        # find the split
        if len(xs) < 2:
            spx, countx = xs.items()[0]
            countx *= 3
            spx, countx = _find_split(xs)
        if len(ys) < 2:
            spy, county = ys.items()[0]
            county *= 3
            spy, county = _find_split(ys)

        # get the grid bounds for the split
        minx = min(xs)
        maxx = max(xs)
        miny = min(ys)
        maxy = max(ys)

        # actually split the polygon and put the results back on the work
        # stack
        if (countx < county and not _single_poly(polygon, 0, minx + 1.0)) or _single_poly(polygon, 1, miny + 1.0):
            work.append(polygon &
                Polygon([(minx, miny), (minx, maxy+1),
                         (spx, maxy+1), (spx, miny)]))
            work.append(polygon &
                Polygon([(spx, miny), (spx, maxy+1),
                         (maxx+1, maxy+1), (maxx+1, miny)]))
            work.append(polygon &
                Polygon([(minx, miny), (minx, spy),
                         (maxx+1, spy), (maxx+1, miny)]))
            work.append(polygon &
                Polygon([(minx, spy), (minx, maxy+1),
                         (maxx+1, maxy+1), (maxx+1, spy)]))

        # Always recurse on the smallest set, which is a trick to ensure that
        # the stack size is O(log n) .
        if work[-2].nPoints() < work[-1].nPoints():

def gpfInfo(fileName):
    Get information on a gpc/gpf file.

        - fileName: name of the file to read
        - contours: number of contours
        - holes: number of holes (if contained)
        - points: total number of points
        - withHoles: file contains hole-flags
    f = open(fileName, 'r')
    contours = int(f.readline())
    holes  = 0
    points = 0
    withHoles = True
    for c in range(contours):
        pp = int(f.readline())
        x = 0
        if c == 0:
            # check for hole-flags
                holes += int(f.readline())
                withHoles = False
                x = 1
            if withHoles:
                holes += int(f.readline())
        [ f.readline() for p in range(pp-x) ]
        points += pp
    return contours, holes, points, withHoles