Verwenden von Matplotlib zum Lösen von Punkten in Polygone

Question

Ich suche nach einem Algorithmus, um zu prüfen, ob ein Punkt innerhalb eines Polygons liegt oder nicht.

Ich verwende derzeit mplPath und contains_point(), aber es scheint in einigen Fällen nicht zu funktionieren.

EDIT 16. September 2016:

Okay, so imporved ich meinen Code, indem Sie einfach, wenn der Punkt, wo die Überprüfung auch an den Rändern. Ich habe immer noch einige Probleme mit dem Rechteck und Fliege Beispiel wenn:

NEW 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) > 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 < 0 : return False 

    squaredlengthba = (b[0] - a[0])*(b[0] - a[0]) + (b[1] - a[1])*(b[1] - a[1]) 
    if dotproduct > 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]) 

Ausgänge von neuem Code:

OLD 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() 

Ich komme mit den folgenden . Wie Sie sehen können, sind die drei rot umrandeten Punkte als außerhalb des Polygons (in blau) angegeben, wenn ich erwarte, dass sie innen liegen.

Ich habe auch versucht, den Radius Wert des Pfades bbPath.contains_points(points, radius = 1.) zu ändern, aber das machte keinen Unterschied.

Jede Hilfe wäre willkommen.

EDIT:

Screenshots aus dem Algorithmus in der Antwort auf diese Frage vorgeschlagen scheinen zu zeigen, dass es auch für andere Fälle versagt.

Answer 1

Okay, also habe ich es endlich geschafft, es formschön zu machen.

#for PIP problem 
import matplotlib.path as mplPath 
import numpy as np 
#for plot 
import matplotlib.pyplot as plt 
import shapely.geometry as shapely 

class MyPoly(shapely.Polygon): 
    def __init__(self,points): 
     super(MyPoly,self).__init__(points) 
     self.points = points 
     self.points_shapely = [shapely.Point(p[0],p[1]) for p in points] 

def convert_to_shapely_points_and_poly(poly,points): 
    poly_shapely = MyPoly(poly) 
    points_shapely = [shapely.Point(p[0],p[1]) for p in points] 
    return poly_shapely,points_shapely 

def plot(poly_init,points_init): 
    #convert to shapely poly and points 
    poly,points = convert_to_shapely_points_and_poly(poly_init,points_init) 

    #plot polygon 
    plt.plot(*zip(*poly.points)) 

    #plot points 
    xs,ys,cs = [],[],[] 
    for point in points: 
     xs.append(point.x) 
     ys.append(point.y) 
     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.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y) 
    if abs(crossproduct) > 0.01 : return False # (or != 0 if using integers) 

    dotproduct = (c.x - a.x) * (b.x - a.x) + (c.y - a.y)*(b.y - a.y) 
    if dotproduct < 0 : return False 

    squaredlengthba = (b.x - a.x)*(b.x - a.x) + (b.y - a.y)*(b.y - a.y) 
    if dotproduct > squaredlengthba: return False 

    return True 

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


def inPoly(poly,point): 
    if poly.contains(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 clockwise 
5 : [[0,0],[20,0],[20,10],[0,10],[0,0]] # test rect counter-clockwise 
} 

#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,2],[30,8]], # rect shape test pts 
5 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,2],[30,8]] # 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]) 

Answer 2

Wenn es der lange Weg ist nicht schlecht, könnten Sie von dem Prinzip gehen: ein Punkt ist in einem Polygon, wenn die Anzahl der Kreuzungen ungerade ist und ähnlich außerhalb, wenn die Anzahl der Kreuzungen sogar ist. Hier ist ein schlampiger Python, den ich zusammen mit dem Testfall, den Sie oben angegeben haben, zusammengestellt habe.

# 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] < x and edge[1][0] < x: 
      continue 
     if edge[0][1] < y and edge[1][1] < 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]))