alternate_tile.py

'''
Copyright 2010 Josiah Carlson
Released into the public domain

alternate_tile.py

This will generate the integer-grid tiles of a provided polygon, similar in
fashion to Polygon.Util.tile(), though this uses a BSP-like algorithm to
partition the space for faster overall execution on large polygons with many
tiles.

Discussed:
http://dr-josiah.blogspot.com/2010/08/binary-space-partitions-and-you.html

Requires the Polygon library:
http://polygon.origo.ethz.ch/

The underlying GPC library requires a license for commercial users:
http://www.cs.man.ac.uk/~toby/alan/software/
'''

from collections import defaultdict
import math

import Polygon

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 tile(polygon):
    '''
    When provided with a Polygon(), this function will yield tiles of the
    original polygon on the integer grid.
    '''
    _int = int
    _floor = math.floor

    work = [polygon]
    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:
            continue

        # 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
            continue

        # find the split
        if len(xs) < 2:
            spx, countx = xs.items()[0]
            countx *= 3
        else:
            spx, countx = _find_split(xs)
        if len(ys) < 2:
            spy, county = ys.items()[0]
            county *= 3
        else:
            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.Polygon([(minx, miny), (minx, maxy+1),
                                 (spx, maxy+1), (spx, miny)]))
            work.append(polygon &
                Polygon.Polygon([(spx, miny), (spx, maxy+1),
                                 (maxx+1, maxy+1), (maxx+1, miny)]))
        else:
            work.append(polygon &
                Polygon.Polygon([(minx, miny), (minx, spy),
                                 (maxx+1, spy), (maxx+1, miny)]))
            work.append(polygon &
                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():
            work.append(work.pop(-2))