Python：查找点是否位于多边形边界上

Question

我有一个点i，并希望创建一个函数来确定该点是否位于多边形的边界上。

使用：

def point_inside_polygon(x, y, poly):
    """Deciding if a point is inside (True, False otherwise) a polygon,
    where poly is a list of pairs (x,y) containing the coordinates
    of the polygon's vertices. The algorithm is called the 'Ray Casting Method'"""
    n = len(poly)
    inside = False
    p1x, p1y = poly[0]
    for i in range(n):
        p2x, p2y = poly[i % n]
        if y > min(p1y, p2y):
            if y <= max(p1y, p2y):
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xinters = (y-p1y) * (p2x-p1x) / (p2y-p1y) + p1x
                    if p1x == p2x or x <= xinters:
                        inside = not inside
        p1x, p1y = p2x, p2y
    return inside

我可以确定点是否在多边形内部。

poly = [(0,0), (2,0), (2,2), (0,2)]
point_inside_polygon(1,1, poly)
True
point_inside_polygon(0,0, poly)
false
point_inside_polygon(2,0, poly)
False
point_inside_polygon(2,2, poly)
True
point_inside_polygon(0,2, poly)
True

如何编写一个函数来判断一个点是否位于多边形的边界上？

- Gianni Spear

---

嗨 - 我修复了代码并恢复了之前被删除的帖子，希望现在它更加清晰易懂。 - Steve Allison

Answer 1

大家都把事情复杂化了。以下是关于多边形的简短见解，假定您有一个距离函数和一个小的EPSILON。

def pointOnPolygon(point, poly):
    for i in range(len(poly)):
        a, b = poly[i - 1], poly[i]
        if abs(dist(a, point) + dist(b, point) - dist(a, b)) < EPSILON:
            return true
    return false

- Rico

Answer 2

将问题分解为三个步骤可能有所帮助：

1. 编写一个函数，可以确定点是否在线段上。
2. 计算构成多边形边界的所有线段。
3. 如果点位于任何线段上，则该点位于边界上。

这里是一些Python代码，假设您已经编写或找到了适合的isPointOnLineSegmentBetweenPoints候选者：

def pointOnPolygon(point, polygonVertices):
    n = len(polygonVertices)
    for i in range(n):
        p1 = polygonVertices[i]
        p2 = polygonVertices[-n+i+1]
        if isPointOnLineSegmentBetweenPoints(point, p1, p2):
            return true
    return false

- Craig Gidney

---

亲爱的Strilanc，感谢你。我在谷歌上找到了一些理论例子，但没有找到优化的Python代码。你有相关链接吗？谢谢。 - Gianni Spear

@Gianni 我已经更新了我的答案中的链接，提供了一个Python示例。当我搜索“python point on line segment”时，它是排名第一的结果。 - Craig Gidney

抱歉受星期五效应影响，但我没有找到任何东西。 - Gianni Spear

谢谢。我正在处理不同的工作。找到点和多边形线i之间的距离，如果为0.0，则表示该点在线i上。 - Gianni Spear

我不理解“isPointOnLineSegmentBetweenPoints”方法的含义？ - Gianni Spear

Answer 3

对于每一对相邻的顶点A和B：

1. 构造一个从A到B的向量，称为p
2. 现在构造一个从A到测试点X的向量，称为q
3. 一对向量的点积公式是p.q = |p||q|cosC，其中C是向量之间的夹角。
4. 因此，如果p.q/|p||q| == 1，则点AX和AB共线。在计算机上工作时，您需要1-p.q/|p||q| < 某个小值，具体取决于您希望精度有多高。
5. 还需要检查|q| < |p|（即X距离A比B更近）

如果4和5都成立，则您的点位于边界上。

编辑

我认为另一种方法是取测试点X，并构建一条垂直于A和B之间线的线。找出这条线和A->B线相交的位置。计算X到这个交点的距离，如果足够小，则将该点视为在线上。

编辑 -- 有趣的小练习！

之前发布了一些错误的代码，因为我记错了一些数学知识。在回家的火车上用Pythonista玩了一下，得到了这个基本可行的代码。由于在iPad上编辑帖子很麻烦，所以省略了数学证明。

没有进行太多测试，也没有测试除零等情况，请用户自行注意。

    # we determine the point of intersection X between
    # the line between A and B and a line through T
    # that is perpendicular to the line AB (can't draw perpendicular
    # in ascii, you'll have to imagine that angle between AB and XT is 90
    # degrees.
    #
    #       B
    #      /
    #.    X  
    #    / \
    #   /   T
    #  A
    # once we know X we can work out the closest the line AB
    # comes to T, if that distance is 0 (or small enough)
    # we can consider T to be on the line
    import math


    # work out where the line through test point t
    # that is perpendicular to ab crosses ab
    #
    # inputs must be 2-tuples or 2-element lists of floats (x,y)
    # returns (x,y) of point of intersection
    def intersection_of_perpendicular(a,b,t):

    if a[0] == b[0]:
            return (a[0],t[1])

    if a[1] == b[1]:
            return (t[0],a[1])

    m = (a[1] - b[1])/(a[0] - b[0]) #slope of ab

    x_inter = (t[1] - a[1] + m*a[0] + (1/m)*t[0])*m/(m**2 + 1)
    y_inter = m*(x_inter - a[0]) + a[1]
    y_inter2 = -(1/m)*(x_inter - t[0]) + t[1]

    #print '...computed ',m,(x_inter, y_inter), y_inter2
    return (x_inter, y_inter)

    # basic Pythagorean formula for distance between two points
    def distance(a,b):
        return math.sqrt( (a[0]-b[0])**2 + (a[1]-b[1])**2 )

    # check if a point is within the box defined by a,b at
    # diagonally opposite corners
    def point_in_box(a,b,t):
        xmin = min(a[0],b[0])
        xmax = max(a[0],b[0])
        ymin = min(a[1],b[1])
        ymax = max(a[1],b[1])

        x_in_bounds = True
        if xmax != xmin:
            x_in_bounds = xmin <= t[0] <= xmax
        y_in_bounds = True
        if ymax != ymin:
            y_in_bounds = ymin <= t[1] <= ymax
        return x_in_bounds and y_in_bounds

    # determine if point t is within 'tolerance' distance
    # of the line between a and b
    # returns Boolean
    def is_on_line_between(a,b,t,tolerance=0.01):
        intersect = intersection_of_perpendicular(a,b,t)
        dist = distance(intersect, t)
        in_bounds = point_in_box(a,b,t)
        return in_bounds and (dist < tolerance)


a = (0,0)
b = (2,2)
t = (0,2)

p = intersection_of_perpendicular(a,b,t)
bounded = point_in_box(a,b,t)
print 'd ',distance(p,t), ' p ',p, bounded

a = (0,2)
b = (2,2)
t = (1,3)

p = intersection_of_perpendicular(a,b,t)
bounded = point_in_box(a,b,t)
print 'd ',distance(p,t),' p ',p, bounded

a = (0.0,2.0)
b = (2.0,7.0)
t = (1.7,6.5)

p = intersection_of_perpendicular(a,b,t)
bounded = point_in_box(a,b,t)
on = is_on_line_between(a,b,t,0.2)
print 'd ',distance(p,t),' p ',p, bounded,on 

- Steve Allison

---

谢谢。我可以请求将其转换为Python代码，因为它不是很清楚:D！ - Gianni Spear

现在有点忙，如果今晚有时间，我可能会试一下。 - Steve Allison

好的，我尝试了第二个算法 - 不能保证，但看起来还不错。这对于周五下午的结束来说非常有趣 ;) - Steve Allison

抱歉Python格式混乱。下次在正式的电脑上会修复它。 - Steve Allison

Answer 4

我还没有测试过这个，但是一般的想法是：

def pointOnBorder(x, y, poly):
    n = len(poly)
    for(i in range(n)):
        p1x, p1y = poly[i]
        p2x, p2y = poly[(i + 1) % n]
        v1x = p2x - p1x
        v1y = p2y - p1y #vector for the edge between p1 and p2
        v2x = x - p1x
        v2y = y - p1y #vector from p1 to the point in question
        if(v1x * v2y - v1y * v2x == 0): #if vectors are parallel 
            if(v2x / v1x > 0): #if vectors are pointing in the same direction
                if(v1x * v1x + v1y * v1y >= v2x * v2x + v2y * v2y): #if v2 is shorter than v1
                    return true
    return false

- Alex

---

但这是针对所有多边形还是仅适用于矩形？ - Gianni Spear

它应该适用于所有多边形。此算法的作用是迭代多边形的边缘，并检查点是否在该边缘上。 - Alex

Answer 5

对于一个凸多边形，你可以通过按顺时针顺序排序顶点并为每个顶点存储顶点与多边形内部点c之间的角度，在O(log n)时间内解决问题。然后对于查询点x，您可以获取从c到x的角度，并进行二分查找以找到唯一的相邻顶点对(v1,v2)，使得x的角度介于v1和v2的角度之间。然后x要么在边缘(v1,v2)上，要么不在边界上。

如果您有一个更复杂的多边形，则可以尝试通过添加一些内部边缘（例如首先三角剖分，然后删除边缘以获得更大的凸多边形）将多边形分解为凸多边形的联合体；如果您获得的凸多边形数量很少（比如k），则可以测试每个凸多边形以查看点是否在边缘上，总运行时间为O(k lg n)，其中n是多边形中顶点的总数。

或者，如果您不担心使用额外的空间，并且您真的想快速确定您是否在边缘上，那么您可以通过沿着每条边添加额外的“顶点”来将每条边分成等距段；很难说需要多少个（听起来像是一个有趣的数学问题），但显然，如果您在每条边上添加足够多的额外顶点，则可以通过在您的顶点集合中找到最近的邻居（原始顶点和您添加的顶点）来确定点必须位于哪条边上，然后只需测试最近邻居顶点所在的一两条边即可。如果您使用二维kd树（您可以将树作为预处理步骤构建，然后树支持快速的k最近邻查询），则可以非常快速地在2-d中找到k最近邻居，而kd树仅使用线性空间。

- user2566092

Answer 6

在这里找到了解决方案：http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html

这是代码：

# Improved point in polygon test which includes edge
# and vertex points

def point_in_poly(x,y,poly):

   # check if point is a vertex
   if (x,y) in poly: return "IN"

   # check if point is on a boundary
   for i in range(len(poly)):
      p1 = None
      p2 = None
      if i==0:
         p1 = poly[0]
         p2 = poly[1]
      else:
         p1 = poly[i-1]
         p2 = poly[i]
      if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
         return "IN"

   n = len(poly)
   inside = False

   p1x,p1y = poly[0]
   for i in range(n+1):
      p2x,p2y = poly[i % n]
      if y > min(p1y,p2y):
         if y <= max(p1y,p2y):
            if x <= max(p1x,p2x):
               if p1y != p2y:
                  xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
               if p1x == p2x or x <= xints:
                  inside = not inside
      p1x,p1y = p2x,p2y

   if inside: return "IN"
   else: return "OUT"

# Test a vertex for inclusion
polygon = [(-33.416032,-70.593016), (-33.415370,-70.589604),
(-33.417340,-70.589046), (-33.417949,-70.592351),
(-33.416032,-70.593016)]
lat= -33.416032
lon= -70.593016

print point_in_poly(lat, lon, polygon)

# test a boundary point for inclusion
poly2 = [(1,1), (5,1), (5,5), (1,5), (1,1)]
x = 3
y = 1
print point_in_poly(x, y, poly2)

- Gustavo Chavez