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) ?

PBrockmann · Accepted Answer · 2023-04-18T09:17:17.070

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

edited Apr 18 '23 at 09:17

answered Apr 17 '23 at 17:25

PBrockmann

303
5
16

You might want to add `import pyvista as pv` to your solution, as you use it later on (`pl = pv.Plotter()`). Also, at least for the 2nd example ("poly_07.vtk") your result is not convex polygons only, so you might want to rephrase your question in this regard (as I understand it from your 2nd example, you are not looking for convex polygons but for the polygons that result from resolving all self-intersections). – simon Apr 18 '23 at 08:57
1

Thank you for the feedback. Added indeed a line for loading pyvista and ok I agree about the title. What would you suggest ? "Separate vtk self intersected polydata into multiple polygons ?" – PBrockmann Apr 18 '23 at 10:15
You're welcome! Regarding the title, it's a bit tricky, I agree. Your suggestion sounds alright, but then again, I would also explain in the question in more detail your problem at hand. In particular: mention that the self-intersections in the given data are marked by duplicate points, because the latter is exactly what your solution (re)solves and I guess it works for you. It is not a general solution for self-intersections (one could, e.g., also imagine self-intersecting polygons without any explicit point at the intersection at all), that's why I would add a bit of detail to the question. – simon Apr 18 '23 at 10:41
What about "Separate vtk self intersected polydata into multiple polygons from duplicate points ?" – PBrockmann Apr 18 '23 at 13:46
Sounds good to me as well! But what I meant is that I would update also the text of your question, not only its heading/title (it mentions the convex polygons several times, and then needs detail on the duplicate points). In any case, just take this as a suggestion, your question is good as it is already! – simon Apr 18 '23 at 13:52

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

1 Answers1

Linked