Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
Copyright 2010 Josiah Carlson
Released into the public domain
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
Requires the Polygon library:
The underlying GPC library requires a license for commercial users:
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:
# 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.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)]))
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():
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment