如何在Python中提取点云内外边界多边形并筛选有效点

如果你的3D点云是平面结构（如图中带孔洞的平面点云），凸包算法只能拿到外边界，要提取所有内外边界多边形，可以按以下步骤实现：

1. 点云投影到2D平面 #

如果是3D点云，先投影到它所在的平面（比如XY平面），得到2D坐标。用Open3D和NumPy就能快速处理：

import numpy as np
import open3d as o3d

# 加载点云文件，替换成你的点云路径
pcd = o3d.io.read_point_cloud("your_point_cloud.pcd")
points = np.asarray(pcd.points)

# 假设点云位于XY平面，直接提取XY坐标作为2D点
points_2d = points[:, :2]

2. 提取2D内外边界多边形 #

用Delaunay三角剖分找出所有边界边，再把边连接成多边形，最后通过面积区分外边界和内边界（孔洞）：

from shapely.geometry import Polygon
from scipy.spatial import Delaunay

# 对2D点做Delaunay三角剖分
tri = Delaunay(points_2d)

# 筛选边界边：只属于一个三角形的边就是边界边
edges = set()
for simplex in tri.simplices:
    for i in range(3):
        # 边按点索引排序，避免重复统计(0,1)和(1,0)
        edge = tuple(sorted((simplex[i], simplex[(i+1)%3])))
        if edge in edges:
            edges.remove(edge)
        else:
            edges.add(edge)

# 将边界边连接成闭合多边形
def connect_boundary_edges(edges):
    polygons = []
    while edges:
        edge = edges.pop()
        poly_indices = [edge[0], edge[1]]
        while True:
            next_idx = None
            # 寻找下一个相连的点
            for e in list(edges):
                if e[0] == poly_indices[-1]:
                    next_idx = e[1]
                    edges.remove(e)
                    break
                elif e[1] == poly_indices[-1]:
                    next_idx = e[0]
                    edges.remove(e)
                    break
            if next_idx is None or next_idx == poly_indices[0]:
                break
            poly_indices.append(next_idx)
        # 转换为坐标点并闭合多边形
        poly_coords = points_2d[poly_indices]
        if not np.array_equal(poly_coords[0], poly_coords[-1]):
            poly_coords = np.vstack([poly_coords, poly_coords[0]])
        polygons.append(poly_coords)
    return polygons

# 获取所有2D多边形
polygons_2d = connect_boundary_edges(edges)

# 区分外边界和内边界：面积最大的是外边界，其余是内孔洞
polygons_2d.sort(key=lambda x: Polygon(x).area, reverse=True)
outer_boundary_2d = polygons_2d[0]
inner_boundaries_2d = polygons_2d[1:]

3. 映射回3D并可视化 #

把2D边界映射回3D空间，用Open3D绘制点云和边界：

# 外边界映射回3D
outer_boundary_3d = np.hstack([
    outer_boundary_2d,
    points[np.where((points[:, :2] == outer_boundary_2d).all(axis=1))][:, 2:3]
])

# 内边界映射回3D
inner_boundaries_3d = []
for inner_poly_2d in inner_boundaries_2d:
    inner_poly_3d = np.hstack([
        inner_poly_2d,
        points[np.where((points[:, :2] == inner_poly_2d).all(axis=1))][:, 2:3]
    ])
    inner_boundaries_3d.append(inner_poly_3d)

# 生成边界线集
outer_line_set = o3d.geometry.LineSet()
outer_line_set.points = o3d.utility.Vector3dVector(outer_boundary_3d)
outer_line_set.lines = o3d.utility.Vector2iVector([[i, i+1] for i in range(len(outer_boundary_3d)-1)])
outer_line_set.paint_uniform_color([1, 0, 0])  # 外边界标红

inner_line_sets = []
for inner_poly in inner_boundaries_3d:
    line_set = o3d.geometry.LineSet()
    line_set.points = o3d.utility.Vector3dVector(inner_poly)
    line_set.lines = o3d.utility.Vector2iVector([[i, i+1] for i in range(len(inner_poly)-1)])
    line_set.paint_uniform_color([0, 1, 0])  # 内边界标绿

# 可视化
o3d.visualization.draw_geometries([pcd, outer_line_set] + inner_line_sets)

4. 过滤点云（移除外边界外/内边界内的点） #

用shapely判断点的位置，保留外边界内且所有内边界外的点：

# 创建多边形对象
outer_polygon = Polygon(outer_boundary_2d)
inner_polygons = [Polygon(poly) for poly in inner_boundaries_2d]

# 过滤点
filtered_points = []
for point_2d, point_3d in zip(points_2d, points):
    point = (point_2d[0], point_2d[1])
    # 先判断是否在外边界内
    if outer_polygon.contains(point):
        # 再判断是否不在任何内边界内
        inside_hole = False
        for inner_poly in inner_polygons:
            if inner_poly.contains(point):
                inside_hole = True
                break
        if not inside_hole:
            filtered_points.append(point_3d)

# 生成过滤后的点云并可视化
filtered_pcd = o3d.geometry.PointCloud()
filtered_pcd.points = o3d.utility.Vector3dVector(np.array(filtered_points))
o3d.visualization.draw_geometries([filtered_pcd])

非平面3D点云的处理方案 #

如果是曲面点云，需要先将点云网格化（用Open3D的泊松重建），再提取网格的边界边：

# 泊松重建生成网格
mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd, depth=9)
# 提取边界线集
boundary_lines = mesh.compute_boundary_loop()
# 可视化
o3d.visualization.draw_geometries([pcd, boundary_lines])

依赖安装 #

执行以下命令安装所需库：

pip install numpy open3d shapely scipy

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