Separate vtk self intersected polydata into multiple polygons from duplicate points?

Question

From a vtk self intersected polydata, I would like to separate it into multiple polygons.

Note that intersections in the initial polygon can be detected from duplicate points in the list of points that form the polygon.

Get the test file from a wget 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)

In fact I would like to describe its coordinates anti-clockwise as requested by the matplotlib path structure.

A solution could be to separate the polydata into 2 convex polygons then use the shapely orient function as needed (https://shapely.readthedocs.io/en/stable/manual.html?highlight=orient#shapely.geometry.polygon.orient).

So how to have 2 convex polygons from this set of lines (polydata) ?

Answer 1

Here is a solution.

For each cell, find duplicate points, and extract polygons.

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)

with poly_07.vtk