Speeding up point-in-polygon for Python

In researching collision detection methods for non-Rect-like objects I came across point-in-polygon -- specifically, the even-odd rule.

The one presently on Wikipedia claims it is written in Python, but it is a little bit pedantic:

# x, y -- x and y coordinates of point
# a list of tuples [(x, y), (x, y), ...]
def isPointInPath(x, y, poly):
num = len(poly)
i = 0
j = num - 1
c = False
for i in range(num):
if  ((poly[i][1] > y) != (poly[j][1] > y)) and \
(x < (poly[j][0] - poly[i][0]) * (y - poly[i][1]) / (poly[j][1] - poly[i][1]) + poly[i][0]):
c = not c
j = i
return c


This didn't particularly remind me of Python, so I spent my lunch break working on this. I don't consider myself a master chef optimizer, but I came up with something slightly more attractive:

def fast_pt_in(pt, poly):
testX, testY = pt
last = -1
for now in xrange(0, len(poly)):
nowX, nowY = poly[now]
lastX, lastY = poly[last]
if (nowY > testY) != (lastY > testY) and (testX < (lastX - nowX) * (testY - nowY)/(lastY - nowY) + nowX):
last += 1


A couple of cProfile tests reveal mine to be faster - both when returning True and False - when done 2**16 times (I figure… 60 FPS, 20 minutes of gameplay, assuming it's brute forcing the collisions without other things to prevent overchecking of course).

In the example below the 'old version' is pt_in and the new version is fast_pt_in. cProfile is executing the following:

cProfile.run("shp = [(x, y) for x in (10, 30) for y in (10, 30)]; [func(point, shape) for x in range(0, 2**16)]")


where func is either function, and point is either one inside the shape (15,15) or outside(5, 5).

Stick\$ python pt_in.py

---TRUE hit in Rect---
=== fast_pt_in ===
131075 function calls in 0.211 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.021    0.021    0.211    0.211 <string>:1(<module>)
65536    0.184    0.000    0.188    0.000 pt_in.py:64(fast_pt_in)
65536    0.004    0.000    0.004    0.000 {len}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
1    0.001    0.001    0.001    0.001 {range}

=== pt_in ===
196611 function calls in 0.312 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.021    0.021    0.312    0.312 <string>:1(<module>)
65536    0.266    0.000    0.291    0.000 pt_in.py:1(pt_in)
65536    0.005    0.000    0.005    0.000 {len}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
65537    0.021    0.000    0.021    0.000 {range}

---FALSE hit in Rect---
=== pt_in ===
196611 function calls in 0.173 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.021    0.021    0.173    0.173 <string>:1(<module>)
65536    0.126    0.000    0.152    0.000 pt_in.py:1(pt_in)
65536    0.005    0.000    0.005    0.000 {len}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
65537    0.021    0.000    0.021    0.000 {range}

=== fast_pt_in ===
131075 function calls in 0.145 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.021    0.021    0.145    0.145 <string>:1(<module>)
65536    0.119    0.000    0.123    0.000 pt_in.py:64(fast_pt_in)
65536    0.004    0.000    0.004    0.000 {len}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
1    0.000    0.000    0.000    0.000 {range}


And that's just with Rect-like shapes; as the number of points in the polygon increases, mine eeks out further and further ahead.

This is all well and good, but I want to make sure I'm getting a lot out of this one. Are there any further changes I can make to this (apart from refactoring it into C, lulz) to eek a little more out of it? I believe I am assigning variables as efficiently as possible to get repeated lookup out of the loop unlike the original, but apart from that I am uncertain if there are other reasonable changes to make.

EDIT: The use of xrange prompts me to clarify - I'm using Python 2.7.

Bonus question -- I'm right, right; that Wiki example is super-unpythonic. Should I replace it? :/

• Mar 7 '14 at 22:21
• Which Python version are you using? The wikipedia version uses the range function which will create a list with all elements in Python version 2.7.x and below. In Python 3+ range returns an iterator function only creating one element during iterator. Your code uses xrange instead, a function only available in Python 2. Try running both in python 3. Mar 7 '14 at 22:51
• Ooh I guess I need to specify, my bad :/ :/ Python 2.7 Mar 7 '14 at 23:03
• Ok try running the wikipedia code with xrange instead of range. Mar 8 '14 at 0:03

Before testing for performance, test for correctness. The way you have defined your shape in the test results in wrong answers:

shp = [(x, y) for x in (10, 30) for y in (10, 30)]


This results in the following shape:

[(10, 10), (10, 30), (30, 10), (30, 30)]


If you plot this, you will notice it looks like this:

(10, 30)  (30, 30)
B----D
|    |
A----C
(10, 10)  (30, 10)


The polygon A-B-C-D is not a convex polygon. What you want instead is A-B-D-C:

shape = [(10, 10), (10, 30), (30, 30), (30, 10)]


Now we have a shape that is a convex polygon, and we can use your two points for verification:

inside = (15, 15)
outside = (5, 5)


To make it easier to compare both functions, I've changed the signature of isPointInPath to be the same as fast_pt_in and just unpack pt to x and y in the first line (the rest of the code is the same):

def isPointInPath(pt, poly):
x, y = pt
# ... rest of the code ...


And because we can possibly have multiple candidates, let's define a global list of functions that we test:

candidates = [isPointInPath, fast_pt_in]


Now, we can use nosetests to write a unit test. What we want to check is that the function returns True when given inside and shape as parameters and False when given outside and shape as parameters (read up on test generators in the nose docs to find out how this test works):

def test_correctness():
def test_call(func, point, shape, result):
assert func(point, shape) == result

for func in candidates:
yield test_call, func, inside, shape, True
yield test_call, func, outside, shape, False


Assuming we have both functions, and those variables (shape, inside, outside, candidates) in a file point_in_poly.py, then we can run it with nosetests (use pip install nose to install nose if you don't have it yet):

% python -m nose point_in_poly.py -v
point_in_poly.test_correctness(<function isPointInPath at 0x10a2aaaa0>, (15, 15), [(10, 10), (10, 30), (30, 30), (30, 10)], True) ... ok
point_in_poly.test_correctness(<function isPointInPath at 0x10a2aaaa0>, (5, 5), [(10, 10), (10, 30), (30, 30), (30, 10)], False) ... ok
point_in_poly.test_correctness(<function fast_pt_in at 0x10ac272a8>, (15, 15), [(10, 10), (10, 30), (30, 30), (30, 10)], True) ... ok
point_in_poly.test_correctness(<function fast_pt_in at 0x10ac272a8>, (5, 5), [(10, 10), (10, 30), (30, 30), (30, 10)], False) ... ok

----------------------------------------------------------------------
Ran 4 tests in 0.001s

OK


Ok, so now we can be sure the functions actually return a correct result for these test cases (if you set shape to your original definition, you see it failing).

Let's get to the performance test. For that, you don't need to use profile, Python provides a nice tool called timeit. At the end of the file add the following:

if __name__ == '__main__':
import timeit
for func in candidates:
duration = timeit.timeit('%s(%r, %r); %s(%r, %r)' % (func.__name__, inside, shape,
func.__name__, outside, shape),
setup='from %s import %s' % (__name__, func.__name__))
print '%-30s: %.5fs' % (func.__name__, duration)


What this does is for each candidate function it calls the function with a point inside and one outside and run this 1000000 times and returns the number of seconds passed, then running it results in:

% python point_in_poly.py
isPointInPath                 : 3.21318s
fast_pt_in                    : 2.50431s


So yeah, fast_pt_in is a bit faster. Now, for the question about increasing performance: You can avoid having to use range together and just remember the last point, then you don't need to index into poly at all:

def faster_pt_in(pt, poly):
testX, testY = pt
last = poly[-1]
for now in poly:
nowX, nowY = now
lastX, lastY = last
if (nowY > testY) != (lastY > testY) and (testX < (lastX - nowX) * (testY - nowY)/(lastY - nowY) + nowX):
last = now


Add this to candidates like so:

candidates = [isPointInPath, fast_pt_in, faster_pt_in]


Use nosetests to see if it still gives the right results (it does) and then run the performance check:

% python point_in_poly.py
isPointInPath                 : 3.22605s
fast_pt_in                    : 2.42880s
faster_pt_in                  : 1.72670s


Now of course there might be other improvements, but with this framework you should now be able to easily add new variants, see if they still do the right thing, and then compare their performance.