Using matplotlib to solve Point in Polygone
<p>I am looking for an algorithm to check if a point is within a polygon or not.</p>

<p>I am currently using mplPath and contains_point() but it doesn't seem to work in some cases.</p>

<p>EDIT 16 Sept 2016:</p>

<p>Okay so I imporved my code by simply checking if the point where also on the edges. I still have some issues with the rectangle and bow-tie example though:</p>

<p>NEW CODE:</p>

<pre><code>#for PIP problem
import matplotlib.path as mplPath
import numpy as np
#for plot
import matplotlib.pyplot as plt


def plot(poly,points):
    bbPath = mplPath.Path(poly)
    #plot polygon
    plt.plot(*zip(*poly))

    #plot points
    xs,ys,cs = [],[],[]
    for point in points:
        xs.append(point[0])
        ys.append(point[1])
        color = inPoly(poly,point)
        cs.append(color)
        print point,":", color
    plt.scatter(xs,ys, c = cs , s = 20*4*2)

    #setting limits
    axes = plt.gca()
    axes.set_xlim([min(xs)-5,max(xs)+50])
    axes.set_ylim([min(ys)-5,max(ys)+10])

    plt.show()

def isBetween(a, b, c): #is c between a and b ?
    crossproduct = (c[1] - a[1]) * (b[0] - a[0]) - (c[0] - a[0]) * (b[1] - a[1])
    if abs(crossproduct) &gt; 0.01 : return False   # (or != 0 if using integers)

    dotproduct = (c[0] - a[0]) * (b[0] - a[0]) + (c[1] - a[1])*(b[1] - a[1])
    if dotproduct &lt; 0 : return False

    squaredlengthba = (b[0] - a[0])*(b[0] - a[0]) + (b[1] - a[1])*(b[1] - a[1])
    if dotproduct &gt; squaredlengthba: return False

    return True

def get_edges(poly):
    # get edges
    edges = []
    for i in range(len(poly)-1):
        t = [poly[i],poly[i+1]]
        edges.append(t)
    return edges

def inPoly(poly,point):
    if bbPath.contains_point(point) == True:
        return 1
    else:
        for e in get_edges(poly):
            if isBetween(e[0],e[1],point):
                return 1
    return 0
# TESTS ========================================================================
#set up poly
polys = {
1 : [[10,10],[10,50],[50,50],[50,80],[100,80],[100,10],[10,10]], # test rectangulary shape
2 : [[20,10],[10,20],[30,20],[20,10]], # test triangle
3 : [[0,0],[0,10],[20,0],[20,10],[0,0]], # test bow-tie
4 : [[0,0],[0,10],[20,10],[20,0],[0,0]] # test rect
}

#points to check
points = {
1 : [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)], # rectangulary shape test pts
2 : [[20,10],[10,20],[30,20],[-5,0],[20,15]] , # triangle  test pts
3 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,5]],  # bow-tie shape test pts
4 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,5]]  # rect shape test pts
}

#print bbPath.contains_points(points) #0 if outside, 1 if inside
for data in zip(polys.itervalues(),points.itervalues()):
    plot(data[0],data[1])
</code></pre>

<p>Outputs from new code:</p>

<p><a href="http://i.stack.imgur.com/tS9I9.png" rel="nofollow"><img src="http://i.stack.imgur.com/tS9I9.png" alt="enter image description here"></a></p>

<p>OLD CODE:</p>

<pre><code>#for PIP problem
import matplotlib.path as mplPath
import numpy as np
#for plot
import matplotlib.pyplot as plt

#set up poly
array = np.array([[10,10],[10,50],[50,50],[50,80],[100,80],[100,10]])
bbPath = mplPath.Path(array)

#points to check
points = [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)]


print bbPath.contains_points(points) #0 if outside, 1 if inside

#plot polygon
plt.plot(*zip(*array))

#plot points
xs,ys,cs = [],[],[]
for point in points:
    xs.append(point[0])
    ys.append(point[1])
    cs.append(bbPath.contains_point(point))
plt.scatter(xs,ys, c = cs)

#setting limits
axes = plt.gca()
axes.set_xlim([0,120])
axes.set_ylim([0,100])

plt.show()
</code></pre>

<p>I come up with the following <a href="http://i.stack.imgur.com/MsJbg.png" rel="nofollow"><img src="http://i.stack.imgur.com/MsJbg.png" alt="graph"></a>. As you can see the three points surrounded in red are indicated as being outside the polygon (in blue) when I would expect them to be inside.</p>

<p>I also tried changing the radius value of the path <code>bbPath.contains_points(points, radius = 1.)</code> but that didn't make any difference.</p>

<p>Any help would be welcome.</p>

<p>EDIT :</p>

<p><a href="http://i.stack.imgur.com/L7Qam.png" rel="nofollow"><img src="http://i.stack.imgur.com/L7Qam.png" alt="enter image description here"></a>screenshots from the algorithm proposed in the answer to this question seem to show that it fails for other cases.</p>

<p>If doing it the long way isn't bad, you could go off of the principle: <strong>a point is inside a polygon if number of intersections is odd and similarly outside if number of intersections is even</strong>. Here's some sloppy python i put together with the test case you gave above.</p>

<pre class="lang-python prettyprint-override"><code># polygon points: [[x,y],...]
points = [[10,10],[10,50],[50,50],[50,80],[100,80],[100,10]]


# get edges
edges = []
for i in range(len(points)-1):
    t = [points[i],points[i+1]]
    edges.append(t)

# get min and max x values to use as the length of ray
xmax = max(points,key=lambda item:item[1])[0]
xmin = min(points,key=lambda item:item[1])[1]
dist = xmax-xmin

# return True if p1,p2,p3 are on the same line
def colinear(p1,p2,p3):
      return (p1[0]*(p2[1] - p3[1]) + p2[0]*(p3[1] - p1[1]) + p3[0]*(p1[1] - p2[1])) == 0

# return True if p1 is on the line segment p2-p3
def inRange(p1,p2,p3):
    dx = abs(p3[0]-p2[0])
    dy = abs(p3[1]-p2[1])
    if abs(p3[0]-p1[0])+abs(p2[0]-p1[0])==dx and abs(p3[1]-p1[1])+abs(p2[1]-p1[1])==dy:
        return True
    return False


# line segment intersection between
# (x1,y1)-(x1+dist,y1) and 
# (x3,y3)-(x4,y4)
def intersect(x1,y1,x3,y3,x4,y4):
    x2 = x1+dist
    B1 = x1-x2
    C1 = B1*y1

    A2 = y4-y3
    B2 = x3-x4
    C2 = A2*x3+B2*y3
    det =-A2*B1
    if(det == 0):
        return False
    x = (B2*C1 - B1*C2)/det
    if inRange((x,y1),(x3,y3),(x4,y4)):
        return True
    return False

# return True if point (x,y) is inside
# polygon defined above
def isInside(x,y):
    i = 0
    for edge in edges:
        # if (x,y) is on the edge return True
        if colinear((x,y),edge[0],edge[1]) and inRange((x,y),edge[0],edge[1]):
            return True
        # if both x values of edge are to the left of (x,y)
        # if both y values of edge are are above or bellow (x,y)
        #    then skip
        if edge[0][0] &lt; x and edge[1][0] &lt; x:
            continue
        if edge[0][1] &lt; y and edge[1][1] &lt; y:
            continue
        # if ray intersects edge, increment i
        if intersect(x,y,edge[0][0],edge[0][1],edge[1][0],edge[1][1]):
            i+=1
    if i%2==1:
        return True
    else:
        return False


l = [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)]
for p in l:
    print(isInside(p[0],p[1]))
</code></pre>

python algorithm python-2.7 matplotlib point-in-polygon