# 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:** 7\
**Page:** 2

<div class="post-metadata">

**Author:** ![marcomusy](https://discourse.vtk.org/user_avatar/discourse.vtk.org/marcomusy/32/95_2.png) [@marcomusy](https://discourse.vtk.org/u/marcomusy)\
**Post date:** [March 28, 2023, 12:56pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/21 "2023-03-28T12:56:11Z")

</div>

The `vtkPolyDataEdgeConnectivityFilter` works fine, it is crashing because it needs polygons, not just lines.  
Here is a solution in vedo, (uncomment `reverse()` to have it clockwise):

```python
from vedo import *
m = Mesh("poly1.vtk").clean().triangulate()
ms = m.split(must_share_edge=True) # return a list of meshes
s1 = ms[0].boundaries().join_segments()[0]#.reverse()
print("First polyline", s1.points())
show(ms, s1, s1.labels2d(), axes=1)

```

```auto
First polyline 
[[1562187.9 -5830366. 2038637.4]
 [1764714.6 -5772294. 2038637.4]
 [1965086. -5707191.5 2038637.4]
 [2163063.2 -5635136. 2038637.4]
 [2358410.5 -5556212.5 2038637.4]
 [2331017.8 -5491532.5 2235859.8]
 [2137944.8 -5569539. 2235859.8]
 [1942269.2 -5640759. 2235859.8]
 [1744224.8 -5705107.5 2235859.8]
 [1544058. -5762504.5 2235859.8]]

```

![Screenshot from 2023-03-28 14-28-58](https://discourse.vtk.org/uploads/default/original/2X/1/1e6956e99edcf49d02c35ec8e5c08bbf68d12849.png)

```auto
pip install -U git+https://github.com/marcomusy/vedo.git

```

---

<div class="post-metadata">

**Author:** ![Paulo\_Carvalho](https://discourse.vtk.org/user_avatar/discourse.vtk.org/paulo_carvalho/32/370_2.png) [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)\
**Post date:** [March 28, 2023, 1:58pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/22 "2023-03-28T13:58:50Z")

</div>

> [@marcomusy](#):
>
> it is crashing because it needs polygons, not just lines.

This is what I meant by some oddity in the dataset causing the crash. Hence, I believe @PBrockmann could join the lines with coincident vertexes to form a single polyline and then process it with `vtkPolyDataEdgeConnectivityFilter `: [vtkStripper not joining contiguous lines to form a single vtkPolyLine.](https://discourse.vtk.org/t/vtkstripper-not-joining-contiguous-lines-to-form-a-single-vtkpolyline/2402)

---

<div class="post-metadata">

**Author:** ![PBrockmann](https://discourse.vtk.org/user_avatar/discourse.vtk.org/pbrockmann/32/524_2.png) [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)\
**Post date:** [March 29, 2023, 9:36pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/23 "2023-03-29T21:36:21Z")

</div>

May I ask you how you “join the lines with coincident vertexes to form a single polyline” ?  
With vtk python API.

@marcomusy : thank you very much for this real test. I didn’t know vedo. I would try to stay witk vtk/pyvista.

---

<div class="post-metadata">

**Author:** ![Paulo\_Carvalho](https://discourse.vtk.org/user_avatar/discourse.vtk.org/paulo_carvalho/32/370_2.png) [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)\
**Post date:** [March 30, 2023, 10:14pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/24 "2023-03-30T22:14:43Z")

</div>

> [@PBrockmann](#):
>
> May I ask you how you “join the lines with coincident vertexes to form a single polyline” ?  
> With vtk python API.

Well, port it. It’s not difficult. The name of the classes and methods are much like the same.

---

<div class="post-metadata">

**Author:** ![Paulo\_Carvalho](https://discourse.vtk.org/user_avatar/discourse.vtk.org/paulo_carvalho/32/370_2.png) [@Paulo\_Carvalho](https://discourse.vtk.org/u/Paulo_Carvalho)\
**Post date:** [March 30, 2023, 10:22pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/25 "2023-03-30T22:22:36Z")

</div>

For example:

```cpp
   vtkSmartPointer<vtkPolyData> poly = vtkSmartPointer<vtkPolyData>::New();
    poly->SetPoints( points );
    poly->SetLines( segments );
    poly->GetCellData()->SetScalars( dataIndexes );

```

Ports to:

```python
   poly = vtk.vtkPolyData()
   poly.SetPoints( points )
   poly.SetLines( segments )
   poly.GetCellData().SetScalars( dataIndexes )

```

---

<div class="post-metadata">

**Author:** ![PBrockmann](https://discourse.vtk.org/user_avatar/discourse.vtk.org/pbrockmann/32/524_2.png) [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)\
**Post date:** [April 3, 2023, 5:47pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/26 "2023-04-03T17:47:50Z")

</div>

I can reproduce what the vedo code from @marcomusy has proposed with the VTK API.  
But the split function from the vedo code gives anyway an uncorrect number of regions.

I have prepared 3 vtk files to play with.  
[https://thredds-su.ipsl.fr/thredds/fileServer/ipsl\_thredds/brocksce/tmp/poly\_11.vtk](https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_11.vtk)  
 ![image](https://discourse.vtk.org/uploads/default/original/2X/5/557bb30ddb8a0349e2f82ef8a83b4273d0247eea.png)

[https://thredds-su.ipsl.fr/thredds/fileServer/ipsl\_thredds/brocksce/tmp/poly\_07.vtk](https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_07.vtk)  
 ![image](https://discourse.vtk.org/uploads/default/original/2X/f/f5b1d63e1aea015368bd73ebe174bcb5029c4b84.png)

[https://thredds-su.ipsl.fr/thredds/fileServer/ipsl\_thredds/brocksce/tmp/poly\_03.vtk](https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_03.vtk)  
 ![image](https://discourse.vtk.org/uploads/default/original/2X/6/6f852e547450d233ee1d1a4cc5e66b030c6aa8b9.png)

In the poly\_11.vtk, I expect to find 2 regions.  
In the poly\_07.vtk, 3 regions.  
In the poly\_03.vtk, 1 region.

```auto
from vedo import *
m = Mesh("poly_07.vtk").clean().triangulate()
ms = m.split(must_share_edge=True) # return a list of meshes
print(len(ms))

```

Gives 1 instead of 3.

---

<div class="post-metadata">

**Author:** ![PBrockmann](https://discourse.vtk.org/user_avatar/discourse.vtk.org/pbrockmann/32/524_2.png) [@PBrockmann](https://discourse.vtk.org/u/PBrockmann)\
**Post date:** [April 20, 2023, 3:00pm UTC](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409/27 "2023-04-20T15:00:08Z")

</div>

Here is the solution I can proposed.

Described in this answer as well.

> <https://stackoverflow.com/questions/75858629/separate-vtk-self-intersected-polydata-into-multiple-polygons-from-duplicate-poi>

Notice that all the different self intersected polygons are separated from their duplicate points into multiple polygons.

![image](https://discourse.vtk.org/uploads/default/original/2X/6/6e844209cc6a9fbf83f353a1988c7f8f51ac3664.jpeg)

```auto
import pyvista as pv
import vtk
import random
import numpy as np

! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/mesh3D_03.vtk
mesh = pv.PolyData('mesh3D_03.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(regions.get_array("RegionId")))

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

allPolys = []
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() 
    
    cleaner = vtk.vtkCleanPolyData()
    cleaner.SetInputData(stripper.GetOutput())
    cleaner.Update()
    
    reg = cleaner.GetOutput()
    print(r, reg.GetNumberOfCells(), reg.GetNumberOfPoints())
    
    #==============================
    polys = []
    for cellIndex in range(reg.GetNumberOfCells()):

        c = reg.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))

        unique, count = np.unique(ids, return_counts=True)
        dup = unique[count > 1]
    
        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(" ", len(polys))
    for p in polys:
        points = []
        for id in p:
            points.append(reg.GetPoint(id))
        m1 = pv.MultipleLines(points)
        print(" add poly")
        allPolys.append(points)
      
        random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
        pl.add_mesh(m1, color=random_color, line_width=4)

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

```

Then my intention was to produce maps with filled polygons in various projections.  
have a look to

> <https://stackoverflow.com/questions/74710358/how-to-detect-inner-polygons-from-a-multipolygon-shapely-object>

The polygon file polys\_03.pickle can be produced from

```auto
import shapely
from pyproj import Transformer

#=======================================================
# From a cartesian 3D CS (geocentric) to a global latitude/longitude CS
trans_GPS_to_XYZ = Transformer.from_crs(4978, 4326, always_xy=True)

#=======================================================
# Define a function to make continuous the polygon
# because the projection transformation produces output
# between -180 and 180 

def makePolygonContinuous(arrayLonLat):
    for i in range(1,len(arrayLonLat)):
        lon1 = arrayLonLat[i-1][0]
        lon2 = arrayLonLat[i][0]
        if (abs(lon2-360-lon1) < abs(lon2-lon1)):
            arrayLonLat[i][0] = lon2-360
        elif (abs(lon2+360-lon1) < abs(lon2-lon1)):
            arrayLonLat[i][0] = lon2+360
        
    return arrayLonLat

#=======================================================
allPolys2 = []
for polyIndex in range(len(allPolys)):
    print(polyIndex, len(allPolys[polyIndex]))
    points = []
    for point in allPolys[polyIndex]:              
        x,y,z = point
        lon, lat, Z = trans_GPS_to_XYZ.transform(x,y,z)
        points.append([lon, lat])

    points = makePolygonContinuous(points)
        
    polygonShapely = shapely.Polygon(points)
    polygonShapelyOriented = shapely.geometry.polygon.orient(polygonShapely)
    
    allPolys2.append(polygonShapelyOriented)

#=======================================================
import pickle
pickle.dump(allPolys2, open("polys_03.pickle","wb"))

```

 ![image](https://discourse.vtk.org/uploads/default/original/2X/0/0e5a694d3ddf0b0bc2c2490471c735cdaf2b466c.png)

Thank you for the different help messages received.

[Previous page](https://discourse.vtk.org/t/create-separated-regions-from-polydata/10409.md?page=1)