[SciPy-user] Triangulation of L-shaped domains

[John Hunter]

>>>>> "Robert" == Robert Kern <robert.kern at gmail.com> writes:

    Robert> Raycast with simulation of simplicity to handle degeneracy
    Robert> is probably your best bet.

    Robert>
    Robert> http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

    Robert> John's "add up the angles" approach is not really a good
    Robert> one. I frequently find it referred to in the literature as
    Robert> "the worst thing you could possibly do."  :-)

OK, I could only stand Robert's taunt for so long.  I added the pnpoly
routine to matplotlib svn extension code: there is a function to test
whether a a single x, y point is inside a polygon
(matplotlib.nxutils.pnpoly) and a function which takes a list of
points and returns a mask for the points inside
(matplotlib.nxutils.points_inside_poly).  When profiling against my
original "worst thing you could possibly do" implementation, this new
version is 15-35 times faster, in addition to being a better algorithm
in terms of handling the degenerate cases.  It is also considerably
faster than the pure numpy pnpoly implementation, since it avoids all
the temporaries and extra passes through the data of doing things at
the numeric level.

For 50 vertices and 1000 candidate inclusion points:

  nxutils.pnpoly vs pure numpy pnpoly: 50x speedup
  nxutils.points_inside_poly vs pure numpy pnpoly: 250x speedup

#!/usr/bin/env python

import time
import matplotlib.nxutils as nxutils
import matplotlib.numerix as nx

import numpy as N

def pnpoly(x, y, verts):
    """Check whether point is in the polygon defined by verts.
    
    verts - 2xN array
    point - (2,) array
    
    See http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html    
    """
    
    verts = verts.astype(float)
    xpi = verts[:,0]
    ypi = verts[:,1]    
    # shift
    xpj = xpi[N.arange(xpi.size)-1]
    ypj = ypi[N.arange(ypi.size)-1]
    
    possible_crossings = ((ypi <= y) & (y < ypj)) | ((ypj <= y) & (y < ypi))

    xpi = xpi[possible_crossings]
    ypi = ypi[possible_crossings]
    xpj = xpj[possible_crossings]
    ypj = ypj[possible_crossings]
    
    crossings = x < (xpj-xpi)*(y - ypi) / (ypj - ypi) + xpi
    
    return sum(crossings) % 2

numtrials, numverts, numpoints = 10, 50, 1000
verts = nx.mlab.rand(numverts, 2)
points = nx.mlab.rand(numpoints, 2)

mask1 = nx.zeros((numpoints, ))
mask2 = nx.zeros((numpoints, ))

t0 = time.time()
for i in range(numtrials):
    for j in range(numpoints):
        x, y = points[j]
        b = pnpoly(x, y, verts)
        mask1[j] = b
told =  time.time() - t0

t0 = time.time()
for i in range(numtrials):
    for j in range(numpoints):
        x, y = points[j]
        b = nxutils.pnpoly(x, y, verts)
        mask2[j] = b
tnew = time.time() - t0

t0 = time.time()
for i in range(numtrials):
    mask3 = nxutils.points_inside_poly(points, verts)
tmany = time.time() - t0
print told, tnew, told/tnew
print told, tmany, told/tmany

for v0, v1, v2 in zip(mask1, mask2, mask3):
    assert(v0==v1)
    assert(v1==v2)



JDH