# Extracting regions as PolyData
**URL:** https://discourse.vtk.org/t/extracting-regions-as-polydata/1788
**Category:** Support
**Created:** [September 19, 2019, 2:41pm UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788 "2019-09-19T14:41:17Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@Srihitha](https://discourse.vtk.org/u/Srihitha)
#### Post date: [September 19, 2019, 2:41pm UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/1 "2019-09-19T14:41:17Z")
Hi,
I am having a polydata containing multiple regions. I want to extract these regions as individual polydata.
Using “vtkPolyDataConnectivityFilter”, I am able to get the total count of extracted regions and able to visualize them as separate regions by setting ColorRegions to ‘On’. Is there any way of accessing the polydata for extracted regions by region id or so?
From the attached image, you can find the poly data having three regions(differentiated by colors).
I need to extract these three regions as three different polydata.
Thanks,
Srihitha.
---
### Author:  [@Ron84](https://discourse.vtk.org/u/Ron84)
#### Post date: [September 19, 2019, 4:00pm UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/2 "2019-09-19T16:00:07Z")
Yes, with SetExtractionModeToSpecifiedRegions, InitializeSpecifiedRegionList and AddSpecifiedRegion.
---
### Author:  [@Srihitha](https://discourse.vtk.org/u/Srihitha)
#### Post date: [September 20, 2019, 5:32am UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/3 "2019-09-20T05:32:18Z")
Thanks, Ronald for the help!
I tried using all the methods you suggested. I am getting each of the individual lines as a region, which is not expected.
---
### Author:  [@Srihitha](https://discourse.vtk.org/u/Srihitha)
#### Post date: [September 20, 2019, 7:18am UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/4 "2019-09-20T07:18:01Z")
Hi,
This is the piece of code I am using to obtain the extracted region’s polydata.
Code snippet:
vtkPolyDataConnectivityFilter\* connectivityFilter = vtkPolyDataConnectivityFilter::New() ;
//The reader’s output will the input polydata attached
connectivityFilter-\>SetInputData(this-\>reader-\>GetOutput());
connectivityFilter-\>SetExtractionModeToAllRegions();
connectivityFilter-\>ScalarConnectivityOff();
connectivityFilter-\>Update();
int regCount = connectivityFilter-\>GetNumberOfExtractedRegions();
//Iterating all the extracted regions to get the region’s polydata
for (int i = 0; i \< regCount; i++)
{
connectivityFilter-\>SetInputData(this-\>reader-\>GetOutput());
connectivityFilter-\>SetExtractionModeToSpecifiedRegions();
connectivityFilter-\>InitializeSpecifiedRegionList();
connectivityFilter-\>AddSpecifiedRegion(i);
connectivityFilter-\>ScalarConnectivityOff();
connectivityFilter-\>Update();
}
InputData:
PS: I am not able to attach .vtk format of InputData. So, just attached an image.
From the above image, I need to obtain the polydata of all the extracted regions differentiated by colors(29 regions).
Thanks,
Srihitha.
---
### Author:  [@Ron84](https://discourse.vtk.org/u/Ron84)
#### Post date: [September 20, 2019, 9:31am UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/5 "2019-09-20T09:31:18Z")
Please upload your vtk file to paste.ee and I will help you.
---
### Author:  [@Srihitha](https://discourse.vtk.org/u/Srihitha)
#### Post date: [September 20, 2019, 12:31pm UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/6 "2019-09-20T12:31:01Z")
Thank you, Ronald for your time!
Sadly, I am not allowed to upload any attachments except images.
So, please find the drive link for the input data.
> **[InputData.vtk](https://drive.google.com/file/d/1LIQfJHZR59-KPskX3GMi7sj0AwvMBb5-/view?usp=sharing)**
>
> Google Drive file.
Thanks,
Srihitha.
---
### Author:  [@Srihitha](https://discourse.vtk.org/u/Srihitha)
#### Post date: [September 24, 2019, 5:16am UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/7 "2019-09-24T05:16:03Z")
Hi,
I hope you are able to download the vtk file, which I shared. If you are having any issues downloading the file, please let me know.
Thanks,
Srihitha.
---
### Author:  [@Ron84](https://discourse.vtk.org/u/Ron84)
#### Post date: [September 26, 2019, 11:34pm UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/8 "2019-09-26T23:34:05Z")
Hello. I’m sorry, I had no time for it. Have you solved it yet?
---
### Author:  [@Srihitha](https://discourse.vtk.org/u/Srihitha)
#### Post date: [September 27, 2019, 11:48am UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/9 "2019-09-27T11:48:09Z")
Hi,
I am trying different ways to solve it. But, couldn’t reach the expected result.
Thanks,
Srihitha.
---
### Author:  [@Ron84](https://discourse.vtk.org/u/Ron84)
#### Post date: [October 4, 2019, 6:32pm UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/10 "2019-10-04T18:32:39Z")
The solution for your problem is this Python script:
```python
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from collections import defaultdict
import math
import vtk
r = vtk.vtkPolyDataReader()
r.SetFileName('InputData.vtk')
c = vtk.vtkCleanPolyData()
c.SetInputConnection(r.GetOutputPort())
c.ToleranceIsAbsoluteOn()
c.SetAbsoluteTolerance(1e-5)
c.Update()
pd = c.GetOutput()
cells = defaultdict(list)
# remove duplicates
for i in range(pd.GetNumberOfCells()):
cell = vtk.vtkIdList()
pd.GetCellPoints(i, cell)
if pd.GetCellType(i) in [vtk.VTK_LINE, vtk.VTK_POLY_LINE]:
ids = [cell.GetId(j) for j in range(cell.GetNumberOfIds())]
if ids[0] > ids[-1]:
ids.reverse()
cells[tuple(ids)].append(i)
to_keep = list(zip(*cells.values()))[0]
_pd = vtk.vtkPolyData()
_pd.DeepCopy(pd)
lines = []
for i in range(_pd.GetNumberOfCells()):
if i not in to_keep:
_pd.DeleteCell(i)
else:
# vertically?
if _pd.GetCellType(i) == vtk.VTK_LINE:
line = vtk.vtkIdList()
_pd.GetCellPoints(i, line)
a = _pd.GetPoint(line.GetId(0))
b = _pd.GetPoint(line.GetId(1))
rev = vtk.vtkIdList()
rev.SetNumberOfIds(2)
rev.SetId(0, line.GetId(1))
rev.SetId(1, line.GetId(0))
if abs(a[0]-b[0]) < 1e-5:
lines.extend([i, _pd.InsertNextCell(vtk.VTK_LINE, rev)])
_pd.BuildLinks()
used = []
for l in lines:
line = vtk.vtkIdList()
_pd.GetCellPoints(l, line)
linksA = vtk.vtkIdList()
linksB = vtk.vtkIdList()
_pd.GetPointCells(line.GetId(0), linksA)
_pd.GetPointCells(line.GetId(1), linksB)
_linksA = [linksA.GetId(i) for i in range(linksA.GetNumberOfIds()) if linksA.GetId(i) not in lines]
_linksB = [linksB.GetId(i) for i in range(linksB.GetNumberOfIds()) if linksB.GetId(i) not in lines]
pA = _pd.GetPoint(line.GetId(0))
pB = _pd.GetPoint(line.GetId(1))
v = [pB[0]-pA[0], pB[1]-pA[1]]
n = -v[1]/abs(v[1])
d = n*pA[0]
poly = []
for _l in _linksA:
if _l in used:
continue
link = vtk.vtkIdList()
_pd.GetCellPoints(_l, link)
pts = [link.GetId(i) for i in range(link.GetNumberOfIds())]
if pts[-1] == line.GetId(0):
pts.reverse()
pt = _pd.GetPoint(pts[-1])
_d = n*pt[0]-d
assert(abs(_d) > 1e-5)
if _d > 0:
poly.extend(reversed(pts[1:]))
used.append(_l)
break
poly.extend([line.GetId(0), line.GetId(1)])
for _l in _linksB:
if _l in used:
continue
link = vtk.vtkIdList()
_pd.GetCellPoints(_l, link)
pts = [link.GetId(i) for i in range(link.GetNumberOfIds())]
if pts[-1] == line.GetId(1):
pts.reverse()
pt = _pd.GetPoint(pts[-1])
_d = n*pt[0]-d
assert(abs(_d) > 1e-5)
if _d > 0:
poly.extend(pts[1:])
used.append(_l)
break
if len(poly) == 2:
continue
if poly[0] == poly[-1]:
del poly[-1]
else:
_links0 = vtk.vtkIdList()
_pd.GetPointCells(poly[0], _links0)
_links1 = vtk.vtkIdList()
_pd.GetPointCells(poly[-1], _links1)
try:
for i in range(_links0.GetNumberOfIds()):
for j in range(_links1.GetNumberOfIds()):
if _pd.GetCellType(_links0.GetId(i)) == vtk.VTK_POLY_LINE \
and _pd.GetCellType(_links1.GetId(j)) == vtk.VTK_POLY_LINE:
if _links0.GetId(i) == _links1.GetId(j):
pl = vtk.vtkIdList()
_pd.GetCellPoints(_links0.GetId(i), pl)
_pl = [pl.GetId(k) for k in range(pl.GetNumberOfIds())]
if _pl[1] == poly[-1]:
_pl.reverse()
poly.extend(_pl[1:-1])
raise StopIteration
elif _links0.GetId(i) != _links1.GetId(j):
cellA = vtk.vtkIdList()
_pd.GetCellPoints(_links0.GetId(i), cellA)
cellB = vtk.vtkIdList()
_pd.GetCellPoints(_links1.GetId(j), cellB)
plA = [cellA.GetId(k) for k in range(cellA.GetNumberOfIds())]
plB = [cellB.GetId(k) for k in range(cellB.GetNumberOfIds())]
if plA[-1] == poly[0]:
plA.reverse()
if plB[-1] == poly[-1]:
plB.reverse()
if plA[-1] != poly[-1] \
and plB[-1] != poly[0] \
and plA[-1] == plB[-1]:
poly.extend(plB[1:])
poly.extend(plA[1:-1][::-1])
raise StopIteration
except StopIteration:
pass
_poly = vtk.vtkIdList()
[_poly.InsertNextId(p) for p in poly]
_pd.InsertNextCell(vtk.VTK_POLYGON, _poly)
for i in range(_pd.GetNumberOfCells()):
if _pd.GetCellType(i) != vtk.VTK_POLYGON:
_pd.DeleteCell(i)
_pd.RemoveDeletedCells()
print(_pd.GetNumberOfCells())
c2 = vtk.vtkCleanPolyData()
c2.SetInputData(_pd)
_if = vtk.vtkIdFilter()
_if.SetInputConnection(c2.GetOutputPort())
_if.CellIdsOn()
_if.PointIdsOff()
w = vtk.vtkPolyDataWriter()
w.SetInputConnection(_if.GetOutputPort())
w.SetFileName('OutputData.vtk')
w.Update()
```
This is the result:
[InputData.vtk](https://discourse.vtk.org/uploads/short-url/opzEUPhI3exOSBGnOKHwnFstgyX.vtk) (87.6 KB)
[OutputData.vtk](https://discourse.vtk.org/uploads/short-url/5wvo5dH3ifrRYqXv7j3iZmPIr0p.vtk) (13.9 KB)
---
### Author:  [@Srihitha](https://discourse.vtk.org/u/Srihitha)
#### Post date: [November 14, 2019, 8:01am UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/11 "2019-11-14T08:01:52Z")
Hi Ronald,
Thank you so much for your time and support.
Looking at your Python script, we converted it into CLI/C++ for our usage. We are able to get NoOfExtractedRegions(polyData-\>GetNumberOfCells()) as 27, but the expected count is 29. We will be glad if you can have a look at our piece of code.
```
vtkPolyDataReader* reader = vtkPolyDataReader::New();
reader->SetFileName("InputData.vtk");
reader->Update();
vtkCleanPolyData* c = vtkCleanPolyData::New();
c->SetInputData(reader->GetOutput());
c->ToleranceIsAbsoluteOn();
c->SetAbsoluteTolerance(1e-5);
c->Update();
vtkPolyData* pd = c->GetOutput();
Dictionary^> cells;
auto pdNumberOfCells = pd->GetNumberOfCells();
for (int i = 0; i < pdNumberOfCells; i++)
{
int cellType = pd->GetCellType(i);
//to identify dupplicate lines shared by neighbouring regions
if (VTK_LINE == cellType || VTK_POLY_LINE == cellType)
{
vtkNew cell;
pd->GetCellPoints(i, cell.GetPointer());
List^ ids = this->FromIdList(cell.GetPointer());
if (ids[0] > ids[ids->Count - 1])
ids->Reverse();//this is only to get unique key to identify duplication
auto key = String::Join("#", ids);
List^ cellIds;
if (!cells.TryGetValue(key, cellIds))
{
cellIds = gcnew List();
cells[key] = cellIds;
}
cellIds->Add(i);//track duplicate index for further processing
}
}
auto zip_list = this->Zip(%cells);
auto to_keep = zip_list[0];
vtkNew _pd;
_pd->DeepCopy(pd);
List lines;
auto _pdNumberOfCells = _pd->GetNumberOfCells();
for (int i = 0; i < _pdNumberOfCells; i++)
{
if (to_keep->IndexOf(i) < 0)
_pd->DeleteCell(i);
else
{
int cellType = _pd->GetCellType(i);
if (VTK_LINE == cellType)
{
vtkNew line;
_pd->GetCellPoints(i, line.GetPointer());
double a[3], b[3];
_pd->GetPoint(line->GetId(0), a);
_pd->GetPoint(line->GetId(1), b);
if (Math::Abs(a[0] - b[0]) < 1e-5)
{
vtkNew rev;
rev->SetNumberOfIds(2);
rev->SetId(0, line->GetId(1));
rev->SetId(1, line->GetId(0));
lines.Add(i);
lines.Add(_pd->InsertNextCell(VTK_LINE, rev.GetPointer()));
}
}
}
}
_pd->BuildLinks();
List used;
vtkNew linksA;
vtkNew linksB;
for each (int l in lines)
{
vtkNew line;
_pd->GetCellPoints(l, line.GetPointer());
linksA->Reset();
linksB->Reset();
_pd->GetPointCells(line->GetId(0), linksA.GetPointer());
_pd->GetPointCells(line->GetId(1), linksB.GetPointer());
List _linksA, _linksB;
auto linksANumberOfIds = linksA->GetNumberOfIds();
for (int i = 0; i < linksANumberOfIds; i++)
{
int id = linksA->GetId(i);
if (lines.IndexOf(id) < 0)
_linksA.Add(id);
}
auto linksBNumberOfIds = linksB->GetNumberOfIds();
for (int i = 0; i < linksBNumberOfIds; i++)
{
int id = linksB->GetId(i);
if (lines.IndexOf(id) < 0)
_linksB.Add(id);
}
double pA[3], pB[3];
_pd->GetPoint(line->GetId(0), pA);
_pd->GetPoint(line->GetId(1), pB);
double v[] = { pB[0] - pA[0], pB[1] - pA[1] };
double n = -v[1] / Math::Abs(v[1]);
double d = n * pA[0];
List poly;
vtkNew link;
for each (int _l in _linksA)
{
if (used.IndexOf(_l) >= 0)
continue;
link->Reset();
_pd->GetCellPoints(_l, link.GetPointer());
List^ pts = this->FromIdList(link.GetPointer());
if (pts[pts->Count - 1] == line->GetId(0))
pts->Reverse();
double pt[3];
_pd->GetPoint(pts[pts->Count - 1], pt);
double _d = n * pt[0] - d;
assert(abs(-d) > 1e-5);
if (_d > 0)
{
for (int k = pts->Count - 1; k >= 1; k--)
poly.Add(pts[k]);//add in revrese order except first item
used.Add(_l);
break;
}
}
poly.Add(line->GetId(0));
poly.Add(line->GetId(1));
for each (int _l in _linksB)
{
if (used.IndexOf(_l) >= 0)
continue;
link->Reset();
_pd->GetCellPoints(_l, link.GetPointer());
List^ pts = this->FromIdList(link.GetPointer());
if (pts[pts->Count - 1] == line->GetId(1))
pts->Reverse();
double pt[3];
_pd->GetPoint(pts[pts->Count - 1], pt);
double _d = n * pt[0] - d;
assert(abs(-d) > 1e-5);
if (_d > 0)
{
for (int k = pts->Count - 1; k >= 1; k --)
poly.Add(pts[k]);
used.Add(_l);
break;
}
}
if (poly.Count == 2)
continue;
if (poly[0] == poly[poly.Count - 1])
poly.RemoveAt(poly.Count - 1);
else
{
vtkNew _links0;
_pd->GetPointCells(poly[0], _links0.GetPointer());
vtkNew _links1;
_pd->GetPointCells(poly[poly.Count - 1], _links1.GetPointer());
try
{
auto _links0NumberOfIds = _links0->GetNumberOfIds();
auto _links1NumberOfIds = _links1->GetNumberOfIds();
for (int i = 0; i < _links0NumberOfIds; i++)
{
for (int j = 0; j < _links1NumberOfIds; j++)
{
int cellType1 = _pd->GetCellType(_links0->GetId(i));
int cellType2 = _pd->GetCellType(_links1->GetId(j));
if (VTK_POLY_LINE == cellType1 && VTK_POLY_LINE == cellType2)
{
if (_links0->GetId(i) == _links1->GetId(j))
{
vtkNew pl;
_pd->GetCellPoints(_links0->GetId(i), pl.GetPointer());
List^ _pl = this->FromIdList(pl.GetPointer());
if (_pl[1] == poly[poly.Count - 1])
_pl->Reverse();
for (int k = 1; k < _pl->Count - 1; k++)
poly.Add(_pl[k]);
throw gcnew Exception();
}
else
{
vtkNew cellA;
_pd->GetCellPoints(_links0->GetId(i), cellA.GetPointer());
vtkNew cellB;
_pd->GetCellPoints(_links1->GetId(j), cellB.GetPointer());
List^ plA = this->FromIdList(cellA.GetPointer());
List^ plB = this->FromIdList(cellB.GetPointer());
if (plA[plA->Count - 1] == poly[0])
plA->Reverse();
if (plB[plB->Count - 1] == poly[poly.Count - 1])
plB->Reverse();
if( (plA[plA->Count -1] != poly[poly.Count-1]) &&
(plB[plB->Count - 1] != poly[0]) &&
(plA[plA->Count - 1] == plB[plB->Count - 1])
)
{
for (int x = 1; x < plB->Count; x ++)
poly.Add(plB[x]);
for (int x = plA->Count -1; x >=1 ; x --)
poly.Add(plA[x]);
throw gcnew Exception();
}
}
}
}
}
}
catch(...)
{
}
vtkNew _poly;
for each (int p in poly)
{
_poly->InsertNextId(p);
}
_pd->InsertNextCell(VTK_POLYGON, _poly.GetPointer());
}
}
_pdNumberOfCells = _pd->GetNumberOfCells();
for (int i = 0; i < _pdNumberOfCells; i++)
{
int cellType = _pd->GetCellType(i);
if (VTK_POLYGON != cellType)
_pd->DeleteCell(cellType);
}
_pd->RemoveDeletedCells();
_pdNumberOfCells = _pd->GetNumberOfCells();
vtkNew c2;
c2->SetInputData(_pd.GetPointer());
vtkNew _if;
_if->SetInputConnection(c2->GetOutputPort());
_if->CellIdsOn();
_if->PointIdsOff();
_if->Update();
auto cpd = vtkPolyData::SafeDownCast(_if->GetOutput());
_pdNumberOfCells = cpd->GetNumberOfCells();
vtkPolyDataWriter* writer = vtkPolyDataWriter::New();
writer->SetFileName("OutputData.vtk");
writer->SetFileTypeToASCII();
writer->SetInputData(_if->GetOutput());
writer->Write();
```
This is the result we are getting in which two regions(Top Right and Bottom Left) are missing.
Please find the drive links for Input and Output data
> **[InputData.vtk](https://drive.google.com/file/d/1LIQfJHZR59-KPskX3GMi7sj0AwvMBb5-/view)**
>
> Google Drive file.
> **[OutputData.vtk](https://drive.google.com/file/d/1fa8aQQC6tUoYk6qZxnlYFTu6NB5xXPgj/view?usp=sharing)**
>
> Google Drive file.
Thanks and regards,
Srihitha.
---
### Author:  [@Ron84](https://discourse.vtk.org/u/Ron84)
#### Post date: [November 15, 2019, 11:54pm UTC](https://discourse.vtk.org/t/extracting-regions-as-polydata/1788/12 "2019-11-15T23:54:34Z")
🧐 It’s difficult to say what’s wrong. You skipped this line `elif _links0.GetId(i) != _links1.GetId(j):`. That’s the only difference I can see. Maybe this is the problem.