搜索

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

生信宝典

共 7264字，需浏览 15分钟

 ·

2021-03-25 22:14

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

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


最简单的方法是目测！

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

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

# (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


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

机器学习

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



浏览 34
点赞
评论
收藏
分享

手机扫一扫分享

分享
举报
评论
图片
表情
推荐
如何快速判断一个行业是否值得加入？
卫星与网络
0
如何快速判断一个行业是否值得加入？
木木老贼
0
如何判断一个区块链项目是否优质？
小张讲区块链
0
算法1.如何判断一个链表是否是回文结构
叶子创业记
0
如何通过JS判断一个数组
前言 如何通过 js 判断一个数组，在面试中我们常常会遇到这个问题，本章主要基于几种判断方式及原理带
如何快速判断一个用户是否访问过我们的 APP？
程序员书单
0
如何快速判断一个用户是否访问过我们的 APP？
每天晒白牙
0
如何快速判断一个用户是否访问过我们的 APP？
码农沉思录
0
必读｜如何判断一个功能的复杂度？
小Q聊产品
0
如何在 Node.js 中判断一个文件是否存在？
程序员成长指北
0
点赞
评论
收藏
分享

手机扫一扫分享

分享
举报