Extracting regions as PolyData

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?

InputData
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.

Yes, with SetExtractionModeToSpecifiedRegions, InitializeSpecifiedRegionList and AddSpecifiedRegion.

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.

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.

Please upload your vtk file to paste.ee and I will help you.

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.

Thanks,
Srihitha.

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.

Hello. I’m sorry, I had no time for it. Have you solved it yet?

Hi,

I am trying different ways to solve it. But, couldn’t reach the expected result.

Thanks,
Srihitha.

The solution for your problem is this Python script:

#!/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 (87.6 KB)
OutputData.vtk (13.9 KB)

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<String^, List<int>^> 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<vtkIdList> cell;
					pd->GetCellPoints(i, cell.GetPointer());
					List<int>^ 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<int>^ cellIds;
					if (!cells.TryGetValue(key, cellIds))
					{
						cellIds = gcnew List<int>();
						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<vtkPolyData> _pd;
			_pd->DeepCopy(pd);
			List<int> 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<vtkIdList> 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<vtkIdList> 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<int> used;
			vtkNew<vtkIdList> linksA;
			vtkNew<vtkIdList> linksB;

			for each (int l in lines)
			{
				vtkNew<vtkIdList> 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<int> _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<int> poly;

				vtkNew<vtkIdList> link;
				for each (int _l in _linksA)
				{
					if (used.IndexOf(_l) >= 0)
						continue;
					link->Reset();
					_pd->GetCellPoints(_l, link.GetPointer());

					List<int>^ 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<int>^ 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<vtkIdList> _links0;
					_pd->GetPointCells(poly[0], _links0.GetPointer());
					
					vtkNew<vtkIdList> _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<vtkIdList> pl;
										_pd->GetCellPoints(_links0->GetId(i), pl.GetPointer());

										List<int>^ _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<vtkIdList> cellA;
										_pd->GetCellPoints(_links0->GetId(i), cellA.GetPointer());

										vtkNew<vtkIdList> cellB;
										_pd->GetCellPoints(_links1->GetId(j), cellB.GetPointer());

										List<int>^ plA = this->FromIdList(cellA.GetPointer());
										List<int>^ 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<vtkIdList> _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<vtkCleanPolyData> c2;
			c2->SetInputData(_pd.GetPointer());

			vtkNew<vtkIdFilter> _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


Thanks and regards,
Srihitha.

:face_with_monocle: 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.