如何将PolyData提取的离散边缘连接为独立陆地特征并输出为PolyData

核心思路 #

针对3D场景下的分离线段连接与独立区域提取，可借助PyVista的拓扑分析能力，通过点的邻接关系拼接线段，再结合面重建、连通性分割得到独立的陆地PolyData。

步骤1：线段拓扑连接 #

将分离线段通过公共点拼接为连续多段线：

1. 构建点到关联线段的映射字典，记录每个点所属的所有线段。
2. 遍历未处理线段，从起点出发，不断寻找共享端点的下一段线段，双向延伸直至无法连接，形成连续PolyLine。
3. 过滤重复线段，保证每条PolyLine唯一且连续。

示例代码：

import pyvista as pv
import numpy as np

# 替换为你的特征边缘PolyData
feature_edges = ...  

# 构建点-线段映射
point_to_cells = {}
for i, cell in enumerate(feature_edges.cells):
    if len(cell) != 3:  # 仅处理Line类型单元格
        continue
    p0, p1 = cell[1], cell[2]
    point_to_cells.setdefault(p0, []).append(i)
    point_to_cells.setdefault(p1, []).append(i)

# 连接线段为PolyLine
visited_cells = set()
polylines = []

for cell_idx in range(len(feature_edges.cells)):
    if cell_idx in visited_cells:
        continue
    cell = feature_edges.cells[cell_idx]
    current_points = [cell[1], cell[2]]
    visited_cells.add(cell_idx)
    
    # 向前延伸线段
    while True:
        last_point = current_points[-1]
        next_candidates = [cid for cid in point_to_cells.get(last_point, []) if cid not in visited_cells]
        if not next_candidates:
            break
        next_cid = next_candidates[0]
        next_cell = feature_edges.cells[next_cid]
        current_points.append(next_cell[2] if next_cell[1] == last_point else next_cell[1])
        visited_cells.add(next_cid)
    
    # 向后延伸线段
    while True:
        first_point = current_points[0]
        prev_candidates = [cid for cid in point_to_cells.get(first_point, []) if cid not in visited_cells]
        if not prev_candidates:
            break
        prev_cid = prev_candidates[0]
        prev_cell = feature_edges.cells[prev_cid]
        current_points.insert(0, prev_cell[1] if prev_cell[2] == first_point else prev_cell[2])
        visited_cells.add(prev_cid)
    
    if len(current_points) >= 2:
        polylines.append([len(current_points)] + current_points)

# 生成连接后的PolyData
connected_edges = pv.PolyData()
connected_edges.points = feature_edges.points
connected_edges.lines = np.array(polylines, dtype=np.int64)

步骤2：闭合线段重建面 #

针对闭合PolyLine，通过投影+三角化生成面（若点集近似共面）：

# 筛选闭合线段
closed_polylines = [pline for pline in polylines if pline[1] == pline[-1]]

# 创建闭合线段PolyData
closed_edges = pv.PolyData()
closed_edges.points = feature_edges.points
closed_edges.lines = np.array(closed_polylines, dtype=np.int64)

# 投影到平面（示例为XY平面，可根据数据调整）
proj_points = closed_edges.points.copy()
proj_points[:, 2] = 0
proj_poly = pv.PolyData(proj_points, closed_edges.lines)

# 三角化生成面并映射回3D坐标
triangulated = proj_poly.triangulate()
triangulated.points = closed_edges.points[triangulated.point_ids]

步骤3：分割独立陆地PolyData #

通过连通性分析分割不同区域：

# 提取所有连通区域
regions = triangulated.connectivity(largest=False)
region_ids = regions.cell_data['RegionId']

# 生成独立陆地PolyData列表
land_polydatas = []
for region_id in np.unique(region_ids):
    mask = region_ids == region_id
    land = regions.extract_cells(mask)
    land_polydatas.append(land)

# 示例：可视化第一个陆地
land_polydatas[0].plot()

注意事项 #

• 若3D点集不共面，可尝试使用reconstruct_surface()方法，或先拟合平面再处理。
• 多分支线段需额外添加分支判断逻辑，避免错误连接。
• 预处理时可调用feature_edges.clean()清理重复点，减少拓扑分析误差。

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