# create separated regions from polydata
**URL:** https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409
**Category:** Support
**Tags:** python
**Created:** [January 9, 2023, 5:32pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409 "2023-01-09T17:32:28Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [January 9, 2023, 5:32pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/1 "2023-01-09T17:32:28Z")
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](https://discourse.vtk.org/t/get-a-continuous-line-from-a-polydata-structure/9864)
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](https://github.com/pyvista/pyvista/discussions/3791)
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.
```auto
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 :
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.
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [January 10, 2023, 2:33pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/2 "2023-01-10T14:33:37Z")
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](https://doi.org/10.1016/j.tcs.2020.07.040) (research paper)
[algorithms - Divide self-intersecting polygon into simple polygons - Computer Science Stack Exchange](https://cs.stackexchange.com/questions/29376/divide-self-intersecting-polygon-into-simple-polygons)
[algorithm - How divide self-intersecting polygon into simple polygons? - Stack Overflow](https://stackoverflow.com/questions/27314724/how-divide-self-intersecting-polygon-into-simple-polygons)
[algorithm - A decomposition of the polygon with self intersections - Stack Overflow](https://stackoverflow.com/questions/53884696/a-decomposition-of-the-polygon-with-self-intersections)
Libraries you can try:
[Overview](http://www.angusj.com/clipper2/Docs/Overview.htm) (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](https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html) (CGAL, a massive C/C++ library of industry-level, high-performing and accurate computational geometry algorithms that has what you need).
take care,
Paulo
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [January 11, 2023, 8:53am UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/3 "2023-01-11T08:53:28Z")
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
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [January 11, 2023, 1:42pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/4 "2023-01-11T13:42:38Z")
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/](https://kitware.github.io/vtk-examples/site/VTKBook/09Chapter9/) .
best,
Paulo
---
### Author:  [@will.schroeder](https://discourse.vtk.org/u/will.schroeder)
#### Post date: [January 11, 2023, 1:49pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/5 "2023-01-11T13:49:14Z")
This is a long shot but may be worth trying, there is a vtkPolyDataEdgeConnectivityFilter that defines connectivity based on shared edges.
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [January 11, 2023, 5:39pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/6 "2023-01-11T17:39:51Z")
Sounds easy to test.
But replacing `vtk.vtkPolyDataConnectivityFilter()` by `vtk.vtkPolyDataEdgeConnectivityFilter()` makes my kernel restarting without any error message.
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [January 11, 2023, 9:47pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/7 "2023-01-11T21:47:51Z")
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).
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [January 12, 2023, 9:18am UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/8 "2023-01-12T09:18:08Z")
You are right.
Running the script gives me a core dump.
I will investigate more.
---
### Author:  [@will.schroeder](https://discourse.vtk.org/u/will.schroeder)
#### Post date: [January 12, 2023, 1:27pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/9 "2023-01-12T13:27:44Z")
If there is a break in vtkPolyDataEdgeConnectivityFilter please let me know, and if possible provide sample data / script.
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [January 12, 2023, 2:55pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/10 "2023-01-12T14:55:43Z")
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()
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [January 19, 2023, 1:10pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/11 "2023-01-19T13:10:29Z")
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](https://stackoverflow.com/questions/8305866/how-do-i-analyze-a-programs-core-dump-file-with-gdb-when-it-has-command-line-pa)
regards,
Paulo
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [March 20, 2023, 6:19pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/12 "2023-03-20T18:19:50Z")
Hi,
Back on this issue that has stayed unresolved.
Any tries or hints to share ?
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [March 21, 2023, 12:01pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/13 "2023-03-21T12:01:27Z")
Hello,
Did you try to investigate the core dump?
best,
PC
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [March 21, 2023, 12:08pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/14 "2023-03-21T12:08:29Z")
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()`
by
`connectivityFilter = vtk.vtkEdgesPolyDataConnectivityFilter()`
is easy to test.
For me it causes a core dumped.
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [March 21, 2023, 2:38pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/15 "2023-03-21T14:38:20Z")
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:
a) waiting;
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.
take care,
Paulo
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [March 21, 2023, 3:11pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/16 "2023-03-21T15:11:12Z")
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 `vtk.vtkStripper()` and `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 ?
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [March 27, 2023, 5:25pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/17 "2023-03-27T17:25:35Z")
I have posted this question related to my problem.
>
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [March 27, 2023, 9:08pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/18 "2023-03-27T21:08:34Z")
> [@PBrockmann](#):
>
> In fact the problem I am facing is how to describe anti-clockwise a polygon.
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): [https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order](https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order)
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.
---
### Author:  [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)
#### Post date: [March 27, 2023, 9:10pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/19 "2023-03-27T21:10:06Z")
If the computed winding order is not that you want, just invert the list of vertexes.
---
### Author:  [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)
#### Post date: [March 27, 2023, 10:10pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/20 "2023-03-27T22:10:04Z")
> 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](https://github.com/shapely/shapely/issues/1705)
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 ?
[Next page](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409.md?page=2)