最近有个朋友问起的一个问题，给定一个荷质比和保留时间构成的平面，想看下哪些化合物落在这个平面里面。
假如有下面一个平面，哪些点在平面里？哪些点在平面外？怎么判断呢？
最简单的方法是目测！
如果采用程序，其中一个方法是光线投射算法，如下图所示：从任意位置画一条到目标点的水平线，计算该水平线进出平面边界的次数，如果是偶数，则点在平面外；如果是奇数，则点在平面内。点正好落在定点或边界时会出现特殊的判断。
我们自己构造一个多边形，绘制几个点：
# (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 最终用的函数
def between(a,x,y,equal=True, epsilon=0.000001):
if equal: return (x=a>=y) or (abs(a-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]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参考
- https://stackoverflow.com/questions/8721406/how-to-determine-if-a-point-is-inside-a-2d-convex-polygon/23223947#23223947
- https://github.com/mikolalysenko/robust-point-in-polygon
- 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?page=1&tab=votes#tab-top
- https://en.m.wikipedia.org/wiki/Point_in_polygon
- 完整代码 https://gitee.com/ct5869/shengxin-baodian
服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net
机房租用，北京机房租用，IDC机房托管， http://www.e1idc.net
