将vtk自相交多数据从重复点分离成多个多边形？

从 vtk 自相交多数据中，我想将其分成多个多边形。

请注意，可以从形成多边形的点列表中的重复点检测初始多边形中的交点。

从 wget 获取测试文件https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_11.vtk https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_11.vtk

then

import pyvista as pv

a = pv.read('poly_11.vtk')

pl = pv.Plotter()
pl.add_mesh(a)
viewer = pl.show(jupyter_backend='trame', return_viewer=True)
display(viewer)

事实上，我想按照 matplotlib 路径结构的要求逆时针描述它的坐标。

解决方案可能是将多数据分成 2 个凸多边形，然后根据需要使用 shapely orient 函数（https://shapely.readthedocs.io/en/stable/manual.html?highlight=orient#shapely.geometry.polygon.orient https://shapely.readthedocs.io/en/stable/manual.html?highlight=orient#shapely.geometry.polygon.orient).

那么如何从这组线（polydata）中得到2个凸多边形呢？

这是一个解决方案。

对于每个单元格，查找重复点并提取多边形。

import numpy as np
import vtk
import pyvista as pv
import random

input_file_name = "poly_07.vtk"
reader = vtk.vtkPolyDataReader()
reader.SetFileName(input_file_name)
reader.Update()

p = reader.GetOutput()
print("Number of cells: ", p.GetNumberOfCells())

polys = []
for cellIndex in range(p.GetNumberOfCells()):

    c = p.GetCell(cellIndex)
    print("-- Cell #%d Number of points: %d"  %(cellIndex, c.GetNumberOfPoints()))
    
    d = c.GetPointIds()

    ids = []
    for idIndex in range(d.GetNumberOfIds()):
        ids.append(d.GetId(idIndex))
    #print(ids)

    # Find duplicate points
    unique, count = np.unique(ids, return_counts=True)
    dup = unique[count > 1]
    print("---- Duplicate points: ", dup)

    # Extract points between duplicate points
    for id in dup[::-1]:
        select = np.where(ids == id)[0]
        polys.append(ids[select[0]:select[1]+1])
        idsIndices = list(range(select[0],select[1]))
        ids = list(np.delete(ids, idsIndices))

print("Number of polygons: ", len(polys))

# Display the results
pl = pv.Plotter()
for poly in polys:
    points = []
    for id in poly:
        points.append(p.GetPoint(id))
    m1 = pv.MultipleLines(points)
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    pl.add_mesh(m1, color=random_color, line_width=2)
viewer = pl.show(jupyter_backend='trame', return_viewer=True)
display(viewer)

与poly_07.vtk