本篇博客给出一种如何判断三维点,也就是一个空间点是否在多边形(三维空间的一个多边形)内的方法。

一、在二维情况下,如何判断一个点在多边形内

我这里使用python实现了文中的第二种方法,也就是将目标点和所有的顶点连接起来,然后计算相邻两边的夹角之和,判断是否等于360°,若是,这个点就在多边形内,若不是,就不再多边形内。根据这个思想,很容易根据下面的Python代码改写为其他语言的代码:

'''
@description 判断点point是否在由顶点数组vertices所指定的多边形内部
思想:将点point和多边形所有的顶点链接起来,计算相邻两边的夹角之和,
若和等于360°,那说明该点就在多边形内。
参考链接:http://www.html-js.com/article/1538
@param  point 待判断的点。有两个分量的List。
@param  vertices 多边形顶点数组,其中的前后相邻的元素在多边形上也
是相邻的。3个以上的二分量List(一个二分量List为一个顶点)组成的List。
@return 若在多边形之内或者在多边形的边界上,返回True,否则返回False
'''
def is_in_2d_polygon(point, vertices):
    px = point[0]
    py = point[1]
    angle_sum = 0

    size = len(vertices)
    if size < 3:
        raise ValueError("len of vertices < 3")
    j = size - 1
    for i in range(0, size):
        sx = vertices[i][0]
        sy = vertices[i][1]
        tx = vertices[j][0]
        ty = vertices[j][1]

        # 通过判断点到通过两点的直线的距离是否为0来判断点是否在边上
        # y = kx + b, -y + kx + b = 0
        k = (sy - ty) / (sx - tx + 0.000000000001)  # 避免除0
        b = sy - k * sx
        dis = fabs(k * px - 1 * py + b) / sqrt(k * k + 1)
        if dis < 0.000001:  # 该点在直线上
            if sx <= px <= tx or tx <= px <= sx:  # 该点在边的两个定点之间,说明顶点在边上
                return True

        # 计算夹角
        angle = atan2(sy - py, sx - px) - atan2(ty - py, tx - px)
        # angle需要在-π到π之内
        if angle >= math.pi:
            angle = angle - math.pi * 2
        elif angle <= -math.pi:
            angle = angle + math.pi * 2

        # 累积
        angle_sum += angle
        j = i

    # 计算夹角和于2*pi之差,若小于一个非常小的数,就认为相等
    return fabs(angle_sum - math.pi * 2) < 0.00000000001

二、如何判断一个三维点点在一个空间中的多边形内

算法思想。

首先计算point和多边形所在平面的距离,若距离大于0,则不在平面上,更不可能在多边形内部;然后把三维平面降为二维平面(直接删除一个分量),把三维点降为二维点,然后根据二维平面内判断点是否在多边形内的方法判断点是否在三维的多边形内。可以把三维平面降为二维平面的理由是,当把一个平面投影到某一坐标面时,若投影前点在多边形内,投影之后还是在多边形内;当然,上述结论的成立条件是,多边形所在平面A不能垂直于投影到的坐标面B。如果法向量中某一分量为0,说明多边形所在平面A垂直于某一坐标面B,如果将多边形投影到坐标面B,投影结果为一条线,就不能判断点是否在多边形内了。

接下来讨论投影方法。

考虑特殊情况,当多边形所在平面A的法向量为(0,0,a),a为一非0常量,则平面A平行于坐标面xOy,因此,删除point和顶点中的z分量,就能得到平面A在坐标面上xOy的投影,此时,可调用二维的函数来判断点是否在多边形内。类似的,当平面A的法向量为(0,a,0)、(a,0,0),a为非零常量,可分量通过删除point和顶点中的y分量或者x分量得到投影。

当平面A的法向量为(a,b,0),a、b为非零常量,则平面A是垂直于坐标面xOy,因此,不能投影到xOy面上,可以投影到另外两个坐标面上。因此,可以通过删除point和顶点中的x或y分量,就能得到一个多边形的投影。

若平面A的法向量的分量均不为0,则任意删除一个分量就可以完成投影。

综上所述,完成投影的方法是,删除point和顶点集中,法向量不为0的任意一个分量。

接下来给出Python代码,大家可根据此Python代码改写为其他语言的代码:

'''
@description 判断空间点point是否在顶点集为vertices,法向量为normal的多边形内

@param point List(3) 三分量List,表示一个三维点。
@param normal List(3) 三分量List,表示一个三维向量。
@param vertices List(n>3)[List(3)] 由n个三维点组成的List。
@return 若在多边形之内或者在多边形的边界上,返回True,否则返回False。
'''
def is_in_3d_polygon(point, normal, vertices):
    # 复制数据
    local_v = []
    for v in vertices:
        local_v.append(v.copy())
    local_p = point.copy()

    # 求点到平面的距离,点到平面的距离公式
    na = normal[0]
    nb = normal[1]
    nc = normal[2]
    d = -(na * local_v[0][0] + nb * local_v[0][1] + nc * local_v[0][2])
    distance = fabs(na * local_p[0] + nb * local_p[1] + nc * local_p[2] + d) \
               / (sqrt(na * na + nb * nb + nc * nc))

    if distance > 0:  # 点不在平面上,肯定不在多边形内
        # print('distance > 0')
        return False

    index = 2  # 默认删除z分量
    if normal[0] != 0:  # 删除x分量
        index = 0
    elif normal[1] != 0:  # 删除y分量
        index = 1
    elif normal[2] != 0:  # 删除z分量
        index = 2
    else:
        raise ValueError('All elem in normal is zero')

    # 删除point和顶点集中指定的分量
    del local_p[index]
    for i in range(0, len(local_v) - 1):
        del local_v[i][index]
    # 调用二维的判断点是否在多边形内的方法。
    return is_in_2d_polygon(local_p, local_v)

三、最后,给出上述两个代码的测试用例

def test_2d():
    point1 = [0, 0]     # 多边形内
    point2 = [-2, -2]   # 多边形顶点
    point3 = [0, -2]    # 多边形边上的点
    point4 = [0, -2.1]  # 多边形外的点

    # 测试三角形
    vertices = [[-2, -2], [2, -2], [0, 4]]
    assert is_in_2d_polygon(point1, vertices)
    assert is_in_2d_polygon(point2, vertices)
    assert is_in_2d_polygon(point3, vertices)
    assert not is_in_2d_polygon(point4, vertices)

    # 测试四边形
    vertices = [[-2, -2], [2, -2], [2, 2], [-2, 2]]
    assert is_in_2d_polygon(point1, vertices)
    assert is_in_2d_polygon(point2, vertices)
    assert is_in_2d_polygon(point3, vertices)
    assert not is_in_2d_polygon(point4, vertices)

    # 测试凹五边形
    vertices = [[-2, -2], [2, -2], [2, 2], [0, 0], [-2, 2]]
    point5 = [0, 1]
    assert is_in_2d_polygon(point1, vertices)
    assert is_in_2d_polygon(point2, vertices)
    assert is_in_2d_polygon(point3, vertices)
    assert not is_in_2d_polygon(point4, vertices)
    assert not is_in_2d_polygon(point5, vertices)
    print('test_2d:all tests passed')


def test_3d():
    p1 = [0, 0, 1]      # 多边形内
    p2 = [-2, -2, 1]    # 多边形顶点
    p3 = [0, -2, 1]     # 多边形边上的点
    p4 = [0, -2.1, 1]   # 多边形外的点

    # 测试一个平行于xOy面的三角形
    v = [[-2, -2, 1], [2, -2, 1], [0, 2, 1]]  # 顶点数组
    n = [0, 0, 1]  # 法向量
    assert is_in_3d_polygon(p1, n, v)
    assert is_in_3d_polygon(p2, n, v)
    assert is_in_3d_polygon(p3, n, v)
    assert not is_in_3d_polygon(p4, n, v)

    # 测试一个倾斜的四边形
    v = [[-2, -2, -1], [2, -2, -1], [2, 2, 1], [-2, 2, 1]]  # 顶点数组
    n = get_normal(v[0], v[1], v[2])  # 法向量
    p1 = [0, 0, 0]      # 多边形内的点
    p2 = [-2, -2, -1]   # 多边形顶点
    p3 = [0, -2, -1]    # 多边形边上的点
    p4 = [0, -2.2, -1]  # 多边形平面外的点
    p5 = [-3, -2, -1]   # 多边形平面上的点,但是不在多边形内
    assert is_in_3d_polygon(p1, n, v)
    assert is_in_3d_polygon(p2, n, v)
    assert is_in_3d_polygon(p3, n, v)
    assert not is_in_3d_polygon(p4, n, v)
    assert not is_in_3d_polygon(p5, n, v)
    print('test_3d: all tests passed')


'''
通过平面上的三个点计算法向量。
'''
def get_normal(point_a, point_b, point_c):
    a = [0, 0, 0]
    b = [0, 0, 0]

    # 计算向量a
    a[0] = point_a[0] - point_b[0]
    a[1] = point_a[1] - point_b[1]
    a[2] = point_a[2] - point_b[2]

    # 计算向量b
    b[0] = point_a[0] - point_c[0]
    b[1] = point_a[1] - point_c[1]
    b[2] = point_a[2] - point_c[2]

    # 通过向量乘法得到法向量
    normal = [0, 0, 0]
    normal[0] = a[1] * b[2] - a[2] * b[1]
    normal[1] = a[2] * b[0] - a[0] * b[2]
    normal[2] = a[0] * b[1] - a[1] * b[0]

    # 如果法向量的三个分量都为0,说明三个点都在直线上
    if normal[0] == 0 and normal[1] == 0 and normal[2] == 0:
        return None
    else:
        return normal



if __name__ == '__main__':
    test_2d()
    test_3d()