Python实现3D无限平面到Bounding Cube的裁剪方法

需求回顾 #

我们需要把由方程 Ax + By + Cz = D 定义的无限3D平面，裁剪到指定的边界立方体（由8个角点描述），最终输出裁剪后平面的四个角点（通常是一个凸四边形）。

核心数学思路 #

裁剪后的平面与立方体的交线是一个凸多边形（大概率是四边形，少数情况是三角形/线段，这里我们聚焦四边形场景），实现步骤如下：

• 计算平面与立方体边的交点：立方体共有12条边，每条边是两个相邻角点的连线。遍历所有边，判断边是否与平面相交，若相交则计算交点坐标。
• 对交点排序：由于交线是凸多边形，我们可以在平面内建立局部坐标系，将交点投影到该坐标系后按极角排序，确保点的顺序是顺时针/逆时针的闭合多边形顺序。
• 输出有效角点：筛选排序后的交点，提取出组成四边形的四个顶点。

Python 实现代码 #

import numpy as np

def line_plane_intersection(line_start, line_end, plane_params):
    """计算线段与平面的交点，若线段与平面无交点则返回None"""
    A, B, C, D = plane_params
    line_vec = np.array(line_end) - np.array(line_start)
    plane_normal = np.array([A, B, C])
    
    # 判断线段与平面是否平行，平行则无交点
    dot_product = np.dot(line_vec, plane_normal)
    if abs(dot_product) < 1e-8:
        return None
    
    # 计算线段上交点的参数t，t∈[0,1]时交点在线段上
    t = -(np.dot(np.array(line_start), plane_normal) + D) / dot_product
    if 0 <= t <= 1:
        intersection = np.array(line_start) + t * line_vec
        return intersection.tolist()
    else:
        return None

def clip_to_bounding_box(plane_params, bbox_points):
    """将无限平面裁剪到边界立方体，返回裁剪后平面的四个角点"""
    # 生成立方体的12条边（仅保留相邻角点的连线）
    edges = []
    for i in range(len(bbox_points)):
        for j in range(i+1, len(bbox_points)):
            # 判断两个点是否为相邻顶点：仅一个坐标分量存在差异
            diffs = [abs(p1 - p2) for p1, p2 in zip(bbox_points[i], bbox_points[j])]
            non_zero_diffs = [d for d in diffs if d > 1e-8]
            if len(non_zero_diffs) == 1:
                edges.append((bbox_points[i], bbox_points[j]))
    
    # 收集所有有效交点并去重
    intersections = []
    for start, end in edges:
        pt = line_plane_intersection(start, end, plane_params)
        if pt is not None:
            # 避免浮点精度导致的重复点
            if not any(np.allclose(pt, existing_pt, atol=1e-8) for existing_pt in intersections):
                intersections.append(pt)
    
    # 处理特殊情况：交点数量不符合四边形要求
    if len(intersections) != 4:
        raise ValueError("平面与立方体的交线不是四边形，请检查输入参数")
    
    # 在平面内建立局部坐标系，对交点按极角排序
    origin = np.array(intersections[0])
    normal = np.array(plane_params[:3])
    
    # 生成平面内的两个正交基向量
    if abs(normal[0]) < abs(normal[1]) and abs(normal[0]) < abs(normal[2]):
        vec1 = np.array([1, 0, 0])
    elif abs(normal[1]) < abs(normal[2]):
        vec1 = np.array([0, 1, 0])
    else:
        vec1 = np.array([0, 0, 1])
    # 将vec1投影到平面内
    vec1 = vec1 - np.dot(vec1, normal) * normal / np.dot(normal, normal)
    vec1 = vec1 / np.linalg.norm(vec1)
    # 第二个基向量由叉乘得到
    vec2 = np.cross(normal, vec1)
    vec2 = vec2 / np.linalg.norm(vec2)
    
    # 转换为局部坐标并排序
    local_coords = []
    for pt in intersections:
        delta = np.array(pt) - origin
        u = np.dot(delta, vec1)
        v = np.dot(delta, vec2)
        local_coords.append((u, v, pt))
    
    # 按极角排序，保证点的顺序为闭合多边形的顺序
    local_coords.sort(key=lambda x: np.arctan2(x[1], x[0]))
    
    # 提取排序后的最终点
    sorted_points = [coord[2] for coord in local_coords]
    
    return sorted_points

使用注意事项 #

1. 输入的 plane_params 对应平面方程 Ax + By + Cz + D = 0，如果你的原始方程是 Ax + By + Cz = D，需要把参数转为 [A, B, C, -D] 传入。
2. bbox_points 是包含8个3D坐标的列表，角点顺序不影响最终结果。
3. 代码中加入了浮点精度阈值（1e-8），避免因计算误差导致的重复点或误判。

内容的提问来源于stack exchange，提问作者molonepa