如何从多数据中提取边缘作为连接特征？

Question

我有一个多数据结构及其提取的边，但使用

extract_feature_edges

函数计算为未连接的单元格（分隔线）。

是否有可能将这些细胞（线）从它们的共同点连接起来，然后获得不同的特征（陆地、岛屿，例如您在图像中看到的——南极洲、澳大利亚……——顺便说一句，它们是古大陆）？

在简历中，我想从我的网格及其边缘中提取不同的土地部分作为单独的多数据。我已经尝试使用 python 模块 shapely 和 polygonize 函数，它可以工作但不能使用 3D 坐标（https://shapely.readthedocs.io/en/latest/reference/shapely.polygonize.html）。

import pyvista as pv

! wget -q -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/pyvista/mesh.vtk
mesh = pv.PolyData('mesh.vtk')
edges = mesh.extract_feature_edges(boundary_edges=True)

pl = pv.Plotter()

pl.add_mesh(pv.Sphere(radius=0.999, theta_resolution=360, phi_resolution=180))
pl.add_mesh(mesh, show_edges=True, edge_color="gray")
pl.add_mesh(edges, color="red", line_width=2)

viewer = pl.show(jupyter_backend='pythreejs', return_viewer=True)
display(viewer)

有什么想法吗？

Answer 1

这是一个使用 vtk.vtkStripper() 将连续线段连接成折线的解决方案。 请参阅 https://discourse.vtk.org/t/get-a-continuous-line-from-a-polydata-structure/9864

中的线程

import pyvista as pv
import vtk
import random

! wget -q -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/pyvista/mesh3D.vtk
mesh = pv.PolyData('mesh3D.vtk')
edges = mesh.extract_feature_edges(boundary_edges=True)

pl = pv.Plotter()

pl.add_mesh(pv.Sphere(radius=0.999, theta_resolution=360, phi_resolution=180))
pl.add_mesh(mesh, show_edges=True, edge_color="gray")

regions = edges.connectivity()
regCount = len(set(regions.get_array("RegionId")))

connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
stripper = vtk.vtkStripper()

for r in range(regCount):
    connectivityFilter.SetInputData(edges)
    connectivityFilter.SetExtractionModeToSpecifiedRegions()
    connectivityFilter.InitializeSpecifiedRegionList()
    connectivityFilter.AddSpecifiedRegion(r)
    connectivityFilter.Update()

    if r == 11:
        p = pv.PolyData(connectivityFilter.GetOutput())
        p.save('poly1.vtk')
    
    stripper.SetInputData(connectivityFilter.GetOutput())
    stripper.SetJoinContiguousSegments(True)
    stripper.Update()
    reg = stripper.GetOutput()
    
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    pl.add_mesh(reg, color=random_color, line_width=4)

viewer = pl.show(jupyter_backend='trame', return_viewer=True)
display(viewer)

Answer 2

这在 github discussions 之前出现过。结论是 PyVista 没有内置任何东西来重新排序边缘，但可能有第三方库可以做到这一点（这个答案提到

libigl

，但我没有这方面的经验）。

我对如何解决这个问题有一些想法，但有人担心这种助手在一般情况下的适用性。但是，在您的特定情况下，我们知道每条边都是一个闭环，并且它们不是很多，因此我们不必担心性能（尤其是内存占用）那么多。

这是一种手动方法，通过构建邻接图并一直走到我们结束每个循环的起点来重新排序边缘：

from collections import defaultdict

import pyvista as pv

# load example mesh
mesh = pv.read('mesh.vtk')

# get edges
edges = mesh.extract_feature_edges(boundary_edges=True)

# build undirected adjacency graph from edges (2-length lines)
# (potential performance improvement: use connectivity to only do this for each closed loop)
# (potentially via calling edges.split_bodies())
lines = edges.lines.reshape(-1, 3)[:, 1:]
adjacency = defaultdict(set)  # {2: {1, 3}, ...} if there are lines from point 2 to point 1 and 3
for first, second in lines:
    adjacency[first].add(second)
    adjacency[second].add(first)

# start looping from whichever point, keep going until we run out of adjacent points
points_left = set(range(edges.n_points))
loops = []
while points_left:
    point = points_left.pop()  # starting point for next loop
    loop = [point]
    loops.append(loop)
    while True:
        # keep walking the loop
        neighb = adjacency[point].pop()
        loop.append(neighb)
        if neighb == loop[0]:
            # this loop is done
            break
        # make sure we never backtrack
        adjacency[neighb].remove(point)
        # bookkeeping
        points_left.discard(neighb)
        point = neighb

# assemble new lines based on the existing ones, flatten
lines = sum(([len(loop)] + loop for loop in loops), [])

# overwrite the lines in the original edges; optionally we could create a copy here
edges.lines = lines

# edges are long, closed loops by construction, so it's probably correct
# plot each curve with an individual colour just to be safe
plotter = pv.Plotter()
plotter.add_mesh(pv.Sphere(radius=0.999))
plotter.add_mesh(edges, scalars=range(edges.n_cells), line_width=3, show_scalar_bar=False)
plotter.enable_anti_aliasing('msaa')
plotter.show()

这段代码用定义每个循环的 14 条更大的线替换了原来的 1760 条 2 长线。不过，您必须要小心一点：在澳大利亚北部，您有一个自相交的循环：

交点出现 4 次而不是 2 次。这意味着我的强力求解器没有给出明确定义的结果：它会随机选择交点，如果（运气不好）我们从交点算法可能会失败。让它更健壮留给读者作为练习（我关于将边缘拆分为单个边缘的评论可能有助于解决这个问题）。