# Possible bug on vtkPolyDataConnectivityFilter

**URL:** https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413
**Category:** Development
**Created:** [March 6, 2019, 3:14am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413 "2019-03-06T03:14:19Z")
**Posts on this page:** 10
**Page:** 1

<div class="post-metadata">

### Author: ![kojikoby](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/k/67e7ee/32.png) [@kojikoby](https://discourse.vtk.org/u/kojikoby)
#### Post date: [March 6, 2019, 3:14am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/1 "2019-03-06T03:14:19Z")

</div>

Hi,  
I am using VTK for medical applications, and encounter a possible bug on vtkPolyConnectivityFilter. I use the function to generate an STL from only the largest connected region. I doubt the spacing parameter of VTK data causes this phenomenon. Your help is greatly appreciated.  
Following code defines a cylinder on 3D space, and generates two STL files from the cylinder. First STL with ordinary axis pitches has only a part of a cylinder, and the second STL with axis pitch of small deviations, difference between float value and double value, has full cylinder.  
The program uses VTK-8.1.1 and runs on Ubuntu 18.04.02 LTS.

![Screenshot%20from%202019-03-06%2011-43-03](https://discourse.vtk.org/uploads/default/original/1X/40cbc84201cee9d423e403c6372f6b6d868930fc.png)  
 ![Screenshot%20from%202019-03-06%2011-43-12](https://discourse.vtk.org/uploads/default/original/1X/e5403fe2c829fbc696886f65b16f32061a526d24.png)

#include “vtkAutoInit.h”  
VTK\_MODULE\_INIT(vtkRenderingOpenGL2); // VTK was built with vtkRenderingOpenGL2  
VTK\_MODULE\_INIT(vtkInteractionStyle);

// vtk headers  
#include \<vtkImageImport.h\>  
#include \<vtkMarchingCubes.h\>  
#include \<vtkPolyDataConnectivityFilter.h\>  
#include \<vtkProperty.h\>  
#include \<vtkSTLReader.h\>  
#include \<vtkSTLWriter.h\>  
// local window display  
#include \<vtkActor.h\>  
#include \<vtkPolyDataMapper.h\>  
#include \<vtkRenderer.h\>  
#include \<vtkRenderWindow.h\>  
#include \<vtkRenderWindowInteractor.h\>

#include   
#include \<unistd.h\>

using namespace std;

//----------------------------------------------------//  
// show stl file //  
//----------------------------------------------------//  
void showStl(string fnstl) {  
// define renderer  
vtkRenderer \*renderer = vtkRenderer::New();  
renderer-\>SetBackground(0.5, 0.5, 0.5); // Background color  
vtkRenderWindow \*renderWindow = vtkRenderWindow::New();  
renderWindow-\>AddRenderer(renderer);  
vtkRenderWindowInteractor \*renderWindowInteractor = vtkRenderWindowInteractor::New();  
renderWindowInteractor-\>SetRenderWindow(renderWindow);

```
// show stl
vtkSTLReader *stlreader = vtkSTLReader::New();
stlreader->SetFileName(fnstl.c_str());
stlreader->Update();
// setup the visualization pipeline
vtkPolyDataMapper *stlmapper = vtkPolyDataMapper::New();
stlmapper->SetInputConnection(stlreader->GetOutputPort());
// actor
vtkActor *stlactor = vtkActor::New();
stlactor->GetProperty()->SetColor(0.5,0.5,0.5);
stlactor->GetProperty()->SetOpacity(0.5);
stlactor->SetMapper(stlmapper);
// set stl to rendere
renderer->AddActor(stlactor);
stlreader->Delete();
stlactor->Delete();

renderWindow->Render();
renderWindowInteractor->Start();

// delete vtk objects
renderer->Delete();
renderWindowInteractor->Delete();
renderWindow->Delete();
stlreader->Delete();
stlmapper->Delete();
stlactor->Delete();

```

}

//----------------------------------------------------//  
// make stl file from voxel data //  
//----------------------------------------------------//  
void makeStl(float \*fpvox, string fnstl, int pix[3], double spc[3]) {  
// Convert the c-style image to a vtkImageData  
vtkImageImport \*imp = vtkImageImport::New();  
imp-\>SetDataSpacing(spc[0], spc[1], spc[2]);  
imp-\>SetDataOrigin(0, 0, 0);  
imp-\>SetWholeExtent(0, pix[0]-1, 0, pix[1]-1, 0, pix[2]-1);  
imp-\>SetDataExtentToWholeExtent();  
imp-\>SetDataScalarTypeToFloat();  
imp-\>SetNumberOfScalarComponents(1);  
imp-\>SetImportVoidPointer(fpvox);  
imp-\>Update();

```
// apply marching cubes
vtkMarchingCubes *sfc = vtkMarchingCubes::New();
sfc->SetInputConnection(imp->GetOutputPort());
sfc->GenerateValues(1, 1, 1);
sfc->Update();

// to remain largest region
vtkPolyDataConnectivityFilter *flt = vtkPolyDataConnectivityFilter::New();
flt->SetInputConnection(sfc->GetOutputPort());
flt->SetExtractionModeToLargestRegion();

// generate stl file (binary)
vtkSTLWriter *stw = vtkSTLWriter::New();
stw->SetInputConnection(flt->GetOutputPort());
stw->SetFileName(fnstl.c_str());
stw->SetFileTypeToBinary();
stw->Write();

// delete vtk classes
imp->Delete();
sfc->Delete();
flt->Delete();
stw->Delete();

```

}

//----------------------------------------------------//  
// entry point //  
//----------------------------------------------------//  
int main(int argc, char \*argv[]) {  
// define sizes  
int pix[3] = {400, 400, 150};  
double spc[3] = {0.626953, 0.626953, 3};

```
// allocate memory
int n = pix[2] * pix[1] * pix[0];	// sizes on axes
float *fpvox = (float *)malloc(n * sizeof(float));	// pitches on axes

// make test data: cylinder
int cntr[2] = {pix[0] / 2, pix[1] / 2};	// cylinder center on x-y plane
int rds = 100;	// radius
for (int k=0; k<pix[2]; k++) {
	for (int j=0; j<pix[1]; j++) {
		for (int i=0; i<pix[0]; i++) {
			int inx = (k * pix[1] + j) * pix[0] + i;
			float r = sqrt((i - cntr[0]) * (i - cntr[0]) + (j - cntr[1]) * (j - cntr[1]));
			if (r <= rds) {
				fpvox[inx] = 1;
			} else {
				fpvox[inx] = 0;
			}
		}
	}
}

// make stl: erroneous
string fnstl0 = "stl0.stl";
makeStl(fpvox, fnstl0, pix, spc);

// show stl
if (fork() == 0) {
	showStl(fnstl0);
}

// make small error on axis pitch
spc[0] = (float)spc[0];
spc[1] = (float)spc[1];
spc[2] = (float)spc[2];

// make stl: correct
string fnstl1 = "stl1.stl";
makeStl(fpvox, fnstl1, pix, spc);

// release memory
free(fpvox);

// show stl
showStl(fnstl1);

return EXIT_SUCCESS;

```

}

---

<div class="post-metadata">

### Author: ![Kenichiro-Yoshimi](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/k/ecc23a/32.png) [@Kenichiro-Yoshimi](https://discourse.vtk.org/u/Kenichiro-Yoshimi)
#### Post date: [March 7, 2019, 1:58am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/2 "2019-03-07T01:58:13Z")

</div>

I think this numerical sensitivity is due to the vtkMarchingCubes rather than the vtkPolyDataConnectivityFilter, since the connectivity filter extracts cells that are only topologically connected and doesn’t know anything about the points’ coordinates.

To safely circumvent this numerical instability, you should use the vtkCleanPolyData which merges duplicate points after applying vtkMarchingCubes.

---

<div class="post-metadata">

### Author: ![lorensen](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lorensen/32/20_2.png) [@lorensen](https://discourse.vtk.org/u/lorensen)
#### Post date: [March 7, 2019, 5:27am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/3 "2019-03-07T05:27:23Z")

</div>

Vtk’s Marching Cubes does not produce duplicate points.

---

<div class="post-metadata">

### Author: ![Kenichiro-Yoshimi](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/k/ecc23a/32.png) [@Kenichiro-Yoshimi](https://discourse.vtk.org/u/Kenichiro-Yoshimi)
#### Post date: [March 7, 2019, 7:49am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/4 "2019-03-07T07:49:13Z")

</div>

Thanks Bill,

In this program, vtkMarchingCubes indeed introduces duplicate points in the output. Although I don’t know the reason why this happens, perhaps the precision for the output points may be single-precision.

---

<div class="post-metadata">

### Author: ![lorensen](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lorensen/32/20_2.png) [@lorensen](https://discourse.vtk.org/u/lorensen)
#### Post date: [March 7, 2019, 8:19am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/5 "2019-03-07T08:19:16Z")

</div>

I’d be surprised if it does. Can you share you program and data?

---

<div class="post-metadata">

### Author: ![Kenichiro-Yoshimi](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/k/ecc23a/32.png) [@Kenichiro-Yoshimi](https://discourse.vtk.org/u/Kenichiro-Yoshimi)
#### Post date: [March 7, 2019, 9:34am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/6 "2019-03-07T09:34:37Z")

</div>

Hi Bill,

I’ve attached CMakeLists.txt and code to replicate the behavior.  
[connect\_test.zip](https://discourse.vtk.org/uploads/default/original/1X/66d6fb5cdd62d4787428c07cf9aecd869083ed08.zip) (2.5 KB)

---

<div class="post-metadata">

### Author: ![lorensen](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lorensen/32/20_2.png) [@lorensen](https://discourse.vtk.org/u/lorensen)
#### Post date: [March 7, 2019, 2:57pm UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/7 "2019-03-07T14:57:19Z")

</div>

Thanks. I can duplicate the problem. Perhaps MC is generating 0 area triangles. I’ll investigate.

Bill

---

<div class="post-metadata">

### Author: ![lorensen](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lorensen/32/20_2.png) [@lorensen](https://discourse.vtk.org/u/lorensen)
#### Post date: [March 7, 2019, 3:34pm UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/8 "2019-03-07T15:34:32Z")

</div>

There may be extra points that are generated but not used. If you print the number of cells, they are the same for both cases.

vtkCleanPolyData will eliminate points that are not used.

Add:  
std::cout \<\< "vtkMarchingCubes: number of cells: " \<\< sfc-\>GetOutput()-\>GetNumberOfCells() \<\< std::endl;

and  
std::cout \<\< "vtkCleanPolyData: number of cells: " \<\< cleaner-\>GetOutput()-\>GetNumberOfCells() \<\< std::e\

ndl;

Result is:  
vtkMarchingCubes: number of points: 121200

vtkMarchingCubes: number of cells: 239592

vtkCleanPolyData: number of points: 120600

vtkCleanPolyData: number of cells: 239592

vtkMarchingCubes: number of points: 120600

vtkMarchingCubes: number of cells: 239592

vtkCleanPolyData: number of points: 120600

vtkCleanPolyData: number of cells: 239592

---

<div class="post-metadata">

### Author: ![Kenichiro-Yoshimi](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/k/ecc23a/32.png) [@Kenichiro-Yoshimi](https://discourse.vtk.org/u/Kenichiro-Yoshimi)
#### Post date: [March 7, 2019, 8:44pm UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/9 "2019-03-07T20:44:19Z")

</div>

Thanks a lot, Bill!  
I can understand this situation well.

---

<div class="post-metadata">

### Author: ![kojikoby](https://discourse.vtk.org/letter_avatar_proxy/v4/letter/k/67e7ee/32.png) [@kojikoby](https://discourse.vtk.org/u/kojikoby)
#### Post date: [March 8, 2019, 12:18am UTC](https://discourse.vtk.org/t/possible-bug-on-vtkpolydataconnectivityfilter/413/10 "2019-03-08T00:18:50Z")

</div>

Many thanks to Kenichiro! My problem was solved by adding vtkCleanPolyData after vtkMarchingCubes exactly as you pointed out. You helped me a lot.  
I was also impressed with high-level discussions exchanged between you and Bill.