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
! 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):
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)
and zoomed here is an example of an extracted region :
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.
My two cents on it is that you need to preprocess geometry data to:
- duplicate vertexes that connect to more than two vertexes;
- 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).
Thank you Paulo for the different possible solutions.
I had imagined a much more simpler from the use of some VTK filters.
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/ .
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.
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.
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()
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
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
Back on this issue that has stayed unresolved.
Any tries or hints to share ?
Did you try to investigate the core dump?
No. I haven’t been able to investigate the core dump.
Honnestly I was expecting that someone confirms that problem.
The solution to change
connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
connectivityFilter = vtk.vtkEdgesPolyDataConnectivityFilter()
is easy to test.
For me it causes a core dumped.
An abend like that does not necessarily happen with others. So, only posting the high-level code that triggers it does not mean that much. Hence, it can be a long time until someone experiences the same problem as yours as these mysterious crashes often have a multitude of causes. Python code depends on piles and piles of lower level APIs. The crash can be starting in any of them. If you wish the community to help you, the least you can do is to get human-readable information on the dump following the steps presented and share them here.
If you don’t feel like getting into the trouble of using a debugger to have an inteligible core dump, I believe your options are:
b) use a work around;
c) hire someone more experienced to do it for you.
Anyway, have you tried a trivial data set? Maybe the crash is caused by some oddity in the data.
I am investigating trying to understand why changing the filter causes an abend.
Have you tried the code I have proposed ? It is quite simple.
I will try to reduce it more and test indeed this vtkPolyDataEdgeConnectivityFilter that causes problem.
In fact the problem I am facing is how to describe anti-clockwise a polygon. The
vtk.vtkPolyDataConnectivityFilter() extract correctly a region. The
stripper.SetJoinContiguousSegments(True) joins each segments as one entity but I need to describe the polygon with an anti-clockwise convention.
How to force this ?
I have posted this question related to my problem.
Use the polygon area method, if the returned area is negative, then it is counter-clockwise (for Cartesian system such as VTK) or clockwise (for inverted Y axis such as screen coordinates): math - How to determine if a list of polygon points are in clockwise order? - Stack Overflow
These methods are reliable in case of self intersecting polygons provided there is a great assimetry in the proportion of clockwise and anti-clockwise portions such as your example just above. One can see that it is “mostly” clockwise or counter-clockwise.
If the computed winding order is not that you want, just invert the list of vertexes.
One can see that it is “mostly” clockwise or counter-clockwise.
This is exactly the problem. See orient function should work on folded polygons · Issue #1705 · shapely/shapely · GitHub
The solution could be to separate the polydata into convex polygons that then can be oriented correctly.
What do you think of this ? And how to do it ?