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

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

最简单的方法是目测！

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

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

# (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 valuesplt.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 polygonnvert = len(bounding_points)    # Points that are insidepoints_inside = []    # For every candidate position within the bounding boxfor idx, pos in enumerate(bounding_box_positions):testx, testy = (pos[0], pos[1])c = 0for i in range(0, nvert):j = i - 1 if i != 0 else nvert - 1if( ((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 polygonif 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 polygonnvert = len(bounding_points)    # Points that are insidepoints_inside = []    # For every candidate position within the bounding boxfor idx, pos in enumerate(bounding_box_positions):testx, testy = (pos[0], pos[1])inside = Falsefor i in range(0, nvert):j = i - 1 if i != 0 else nvert - 1# Point on horizontaol edgeif (verty[i] == testy ) and (verty[j] == testy) and between(testx, vertx[i],vertx[j]):inside = edgeAsInside                break# corner pointif (((verty[i] == testy) and (vertx[i] == testx)) or ((verty[j] == testy) and (vertx[j]==testx))):inside = edgeAsInside                breakif between(testy, verty[i], verty[j]):                if(testy==verty[i]<=verty[j]) or (testy==verty[j]<=verty[i]):                    continuec = (vertx[i]-testx)*(verty[j]-testy) - (vertx[j]-testx)*(verty[i]-testy)                # print(c)if c==0 or abs(c) < epsilon:inside = edgeAsInside                    breakif (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

往期精品(点击图片直达文字对应教程)

机器学习

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集