create separated regions from polydata

Hi,

I describe here a problem I am facing to. It follows a question that I had asked here as: get a continuous line from a polydata structure

Today, I would like to produce maps with filled polygons (that may have holes).
The problem is described here: Make polygons valid and ready to be plotted as matplotlib patches · Discussion #3791 · pyvista/pyvista · GitHub

But the problem could be solved from a VTK point of view by creating separate regions even if they are connected by a single vertex (extract_feature_edges). I have not been able to do that
and this is the purpose of my question.

Here is snippet for testing.

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

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

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

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

for r in range(regCount):
    connectivityFilter.SetInputData(edges)
    connectivityFilter.SetExtractionModeToSpecifiedRegions()
    connectivityFilter.InitializeSpecifiedRegionList()
    connectivityFilter.AddSpecifiedRegion(r)
    connectivityFilter.Update()
    
    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=2)

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

It produces:

and zoomed here is an example of an extracted region :
image
that I would like to separate and be seen as 2 distinct regions. I will then be able to describe for each polygon their coordinates anti-clockwise as requested by the matplotlib path structure.

Hello,

My two cents on it is that you need to preprocess geometry data to:

  1. duplicate vertexes that connect to more than two vertexes;
  2. for each of the resulting pairs, set the connectivity of the two collocated vertexes such that you have two closed polygons and not a single self intersecting one.

Not a trivial task though, so I suggest refering to these:
https://doi.org/10.1016/j.tcs.2020.07.040 (research paper)
algorithms - Divide self-intersecting polygon into simple polygons - Computer Science Stack Exchange
algorithm - How divide self-intersecting polygon into simple polygons? - Stack Overflow
algorithm - A decomposition of the polygon with self intersections - Stack Overflow

Libraries you can try:
Overview (Clipper2, a small library with an algorithm that does what you need for C++, C# and Delphi).
CGAL 5.5.1 - 2D Arrangements: User Manual (CGAL, a massive C/C++ library of industry-level, high-performing and accurate computational geometry algorithms that has what you need).

take care,

Paulo

1 Like

Hi,
Thank you Paulo for the different possible solutions.
I had imagined a much more simpler from the use of some VTK filters.
Regards
Patrick

Hello,

Well, it’s likely doable in VTK only by combining two or more filters, but off the top of my head I don’t know how. Please, take a look at this chapter about VTK’s advanced CG algorithms: https://kitware.github.io/vtk-examples/site/VTKBook/09Chapter9/ .

best,

Paulo

This is a long shot but may be worth trying, there is a vtkPolyDataEdgeConnectivityFilter that defines connectivity based on shared edges.

Sounds easy to test.

But replacing vtk.vtkPolyDataConnectivityFilter() by vtk.vtkPolyDataEdgeConnectivityFilter() makes my kernel restarting without any error message.

I assume you’re running in a Python notebook. Can you try running it as a Python script? Running as notebooks often results in important feedback being ignored (e.g. error messages).

You are right.

Running the script gives me a core dump.
I will investigate more.

If there is a break in vtkPolyDataEdgeConnectivityFilter please let me know, and if possible provide sample data / script.

See above.
Please download the mesh3D.vtk file with a:
wget -q -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/pyvista/mesh3D.vtk

The script runs with vtkPolyDataConnectivityFilter(), not with vtkPolyDataEdgeConnectivityFilter()

Hello,

So, a core dump means you had an abend. Besides the core dump, was there any higher-level messages (e.g. Python stack trace or OpenGL error messages)? If not, even a core dump can provide some insight to what went wrong depending on where exactly the crash ensued.

Can you identify which executable crashed? Depending on the answer, it’ll be easy or difficult to analyze the dump. For example, if the crash occured inside python or libvtkCommonCore-9.1.dll you can get debug versions of these with relative ease. If the crash took place inside one of the OpenGL backend libraries, you can have a hard time getting the debug version of them because they are supplied by your graphics card vendor.

Once you have the debug version of the binary that crashed, you can analyze the core dump with gdb: linux - How do I analyze a program's core dump file with GDB when it has command-line parameters? - Stack Overflow

regards,

Paulo