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.
It’s difficult to say what’s wrong. You skipped this line