如何判断一个点是否落在一个平面内

最近有个朋友问起的一个问题，给定一个荷质比和保留时间构成的平面，想看下哪些化合物落在这个平面里面。

假如有下面一个平面，哪些点在平面里？哪些点在平面外？怎么判断呢？

最简单的方法是目测！

如果采用程序，其中一个方法是光线投射算法，如下图所示：从任意位置画一条到目标点的水平线，计算该水平线进出平面边界的次数，如果是偶数，则点在平面外；如果是奇数，则点在平面内。点正好落在定点或边界时会出现特殊的判断。

我们自己构造一个多边形，绘制几个点：

# (1.5,2)bounding_points = [(1,2),(2,1),(3,1.5),(3.5,2.5),(2.1,6.3),(1,3)]
bounding_box_positions = [(1.5,2), (1.5,4.5), (3,5)]

如下图

plotPolygenAndPoints(bounding_points, bounding_box_positions)

加上从左侧发过来的光线

plotPolygenAndPoints(bounding_points, bounding_box_positions, ray=True)

原始函数不能区分部分边界上的点

Stack Overflow上提供了一个函数polygon_ray_casting (代码见文末)实现ray casting算法，被转换成了各个程序语言的版本。

测试下这个函数，一个边界点(1.5,4.5)计算为在平面内，另一个边界点(2.5,1.25)却未落在边界内。

bounding_box_positions = [(1.5,2), (1,1.25),(1.5,4.5), (3,3),(4,3.1),(2.5,1.25),(3,5),(1.0,2),(3.5,2.5)]
bounding_box_positions = sorted(bounding_box_positions, key=lambda x:x[0])
polygon_ray_casting(bounding_points, bounding_box_positions)

[(3, 3), (3.5, 2.5)]

关于这个结果在评论区也有讨论，函数设计的目的不同，行为也会不同。在公共平台找到的代码和方法是给我们提供思路的，有些代码可以直接用，有些代码需要修改后用。使用前，要注意其使用条件，测试后再用于自己的数据。

修改后的函数可处理边界点

增加边界点的判断、定点的判断增加边界点的判断、定点的判断

测试凸多边形

# 默认落在边界也算polygon_ray_casting_modi(bounding_points, bounding_box_positions)

[(1.0, 2), (3, 3), (3.5, 2.5)]

plotPolygenAndPoints(bounding_points, bounding_box_positions)

# 落在边上不算在内polygon_ray_casting_modi(bounding_points, bounding_box_positions, edgeAsInside=False)

[(1.5, 2), (3, 3)]

测试凹多边形

# (1.5,2)bounding_points = [(1,2),(2,2.5),(3,1.5),(4,2),(2.1,4),(1,3)]
bounding_box_positions = [(1.5,2), (1,1.25), (1.5,2.5), (3,1.5), (3,3),(4,3.1), (1,2.1), (3.5,1.75)]
bounding_box_positions = sorted(bounding_box_positions, key=lambda x:x[0])
polygon_ray_casting_modi(bounding_points, bounding_box_positions)

[(1, 2.1), (1.5, 2.5), (3, 1.5), (3, 3), (3.5, 1.75)]

plotPolygenAndPoints(bounding_points, bounding_box_positions)

测试多种六边形

bounding_points = [(1,3),(2,3),(3,2),(2,1),(1,1),(0,2)]
bounding_box_positions = [(-1,3), (1.5,2.9), (1.5,3),(0,2),(1,3),(0.5,1.5),(2.5,1.5), (2.5,3),(0,1.5),(1,1.1)]
bounding_box_positions = sorted(bounding_box_positions, key=lambda x:x[0])
polygon_ray_casting_modi(bounding_points, bounding_box_positions)

[(0, 2), (0.5, 1.5), (1, 3), (1, 1.1), (1.5, 2.9), (1.5, 3), (2.5, 1.5)]

plotPolygenAndPoints(bounding_points, bounding_box_positions)

# (1.5,2)bounding_points = [(1,3),(2,3),(2,2),(3,2),(3,1),(1,1)]
bounding_box_positions = [(1.5,2), (1.5,3), (1,3), (3,1.25), (3,3)]
polygon_ray_casting_modi(bounding_points, bounding_box_positions)

[(1.5, 2), (1.5, 3), (1, 3), (3, 1.25)]

plotPolygenAndPoints(bounding_points, bounding_box_positions)

# (1.5,2)bounding_points = [(1,2),(3,2),(3,3),(4,3),(4,1),(1,1)]
bounding_box_positions = [(1.5,2), (1,1.25), (3,3),(4,3.1)]
polygon_ray_casting_modi(bounding_points, bounding_box_positions)

[(1.5, 2), (1, 1.25), (3, 3)]

plotPolygenAndPoints(bounding_points, bounding_box_positions)

for i in bounding_box_positions:
    print(i)
    print(wn_PnPoly(i, bounding_points))

(1.5, 2)
0
(1, 1.25)
-1
(3, 3)
0
(4, 3.1)
0

用到的绘图函数

import matplotlib.pyplot as pltdef plotPolygenAndPoints(bounding_points, bounding_box_positions, ray=False):
    coord = bounding_points[:]
    coord.append(coord[0]) #repeat the first point to create a 'closed loop'

    xs, ys = zip(*coord) #create lists of x and y values

    plt.figure()
    plt.plot(xs,ys) 
    for x,y in bounding_box_positions:
        plt.plot(x,y,'o')        if ray:
            plt.plot((0,x),(y,y))
            plt.arrow(x, y, 0.3, 0, shape='full', lw=1, length_includes_head=True, head_width=.1)

    plt.show()

用到的判断函数

def polygon_ray_casting(bounding_points, bounding_box_positions):
    # https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
    # https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon
    # Arrays containing the x- and y-coordinates of the polygon's vertices.
    vertx = [point[0] for point in bounding_points]
    verty = [point[1] for point in bounding_points]    # Number of vertices in the polygon
    nvert = len(bounding_points)    # Points that are inside
    points_inside = []    # For every candidate position within the bounding box
    for idx, pos in enumerate(bounding_box_positions):
        testx, testy = (pos[0], pos[1])
        c = 0
        for i in range(0, nvert):
            j = i - 1 if i != 0 else nvert - 1
            if( ((verty[i] > testy ) != (verty[j] > testy))   and
                    (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
                c += 1
        # If odd, that means that we are inside the polygon
        if c % 2 == 1: 
            points_inside.append(pos)    return points_inside

最终用的函数

def between(a,x,y,equal=True, epsilon=0.000001):
    if equal:        return (x<=a<=y) or (x>=a>=y) or (abs(a-x)<epsilon) or (abs(a-y)<epsilon)    else:        return (x<a<y) or (x>a>y)def polygon_ray_casting_modi(bounding_points, bounding_box_positions, edgeAsInside=True, epsilon=0.000001):
    # https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
    # https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon
    # Arrays containing the x- and y-coordinates of the polygon's vertices.
    vertx = [point[0] for point in bounding_points]
    verty = [point[1] for point in bounding_points]    # Number of vertices in the polygon
    nvert = len(bounding_points)    # Points that are inside
    points_inside = []    # For every candidate position within the bounding box
    for idx, pos in enumerate(bounding_box_positions):
        testx, testy = (pos[0], pos[1])

        inside = False
        for i in range(0, nvert):
            j = i - 1 if i != 0 else nvert - 1

            # Point on horizontaol edge
            if (verty[i] == testy ) and (verty[j] == testy) and between(testx, vertx[i],vertx[j]):
                inside = edgeAsInside                break

            # corner point
            if (((verty[i] == testy) and (vertx[i] == testx)) or ((verty[j] == testy) and (vertx[j]==testx))):
                inside = edgeAsInside                break

            if between(testy, verty[i], verty[j]):                if(testy==verty[i]<=verty[j]) or (testy==verty[j]<=verty[i]):                    continue
                c = (vertx[i]-testx)*(verty[j]-testy) - (vertx[j]-testx)*(verty[i]-testy)                # print(c)
                if c==0 or abs(c) < epsilon:
                    inside = edgeAsInside                    break
                if (verty[i]<verty[j]) == (c>0):
                    inside = not inside        # If odd, that means that we are inside the polygon
        #print(inside)
        if inside:
            points_inside.append(pos)    return points_inside

参考

1. https://stackoverflow.com/questions/8721406/how-to-determine-if-a-point-is-inside-a-2d-convex-polygon/23223947#23223947
2. https://github.com/mikolalysenko/robust-point-in-polygon
3. https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html
4. https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon?page=1&tab=votes#tab-top
5. https://en.m.wikipedia.org/wiki/Point_in_polygon
6. 完整代码 https://gitee.com/ct5869/shengxin-baodian