vtkMarchingCubes / vtkDiscreteMarchingCubes does not produce a closed mesh surface

Hello Everyone.

I have a simple three dimensional vtkImageData object representing a voxelgrid:

If a voxel at an index (i, j, k) is set, the entry is 1 ( voxelImageData[idx(x, y, z)] = 1), if the voxel at this index is not set, then the entry is 0.

I want to convert this voxelgrid (3D-volume) to a regular vtkPolyData-Object (mesh surface) for STL-Export with the vtkStlWriter.

I try to archive this using the vtkMarchingCubes / vtkDiscreteMarchingCubes / vtkFlyingEdges3D - Algorithm:

vtkNew<vtkDiscreteMarchingCubes> surface;
//vtkNew<vtkMarchingCubes> surface;
//vtkNew<vtkFlyingEdges3D> surface;
surface->SetValue(0, 1.0);
//surface->GenerateValues(1, 1, 1);

While it succeeds in generating a mesh-surface which represents the surface of the occupied voxels, the resulting mesh-surfaces has holes (mostly axis-parallel surfaces are missing).

For example, given this voxel grid as input:

The output should be simply a cube with all sides closed.

However, this is what i get:

A cube surface, where only half of the sides are closed.

I tried this with a variety of voxel representations and mostly axis parallel surfaces are almost always missing (but there are also additional holes at times.)

What can i do to get a closed surface mesh from a binary voxel grid (labelmap) using vtk?

Thank you very much in advance.

Make sure you have a margin of at least one empty voxel around the filled voxels.

1 Like

Hello Andras.

Thank you very much for your reply. I will try it.

Can i use the vtkImageConstantPad for this or should i do it on my own?


Hello again Andras.

I did what you suggested.
And it works like charm!
Thank you!

Here is what i did:

vtkNew<vtkImageChangeInformation> translator;
translator->SetExtentTranslation (1, 1, 1);
translator->SetOriginTranslation(-spacing[0], -spacing[1], -spacing[2]);

vtkNew<vtkImageConstantPad> pad;
pad->SetInputConnection(translator->GetOutputPort ());
pad->SetOutputWholeExtent(0, resolution + 1,
                          0, resolution + 1,
                          0, resolution + 1);

And that did the trick!

Not sure if I have missed something, but to transform your image data into a surface you can also
use the vtkGeometryFilter.

Hello Charles.

Thank you very much for the suggestion.

But i have tried it with the vtkGeometryFilter and it does not work. When i export the result with the vtkSTLWriter i don’t get a valid output. I think that since the vtkGeometryFilter has no option for specifing the isosurface (that is the only thing i got here, since it is in VTK-terms a “binary labelmap”) it simply tries to output anything.

With the above solution (using the marching cubes or flying edges algorithm and having a margin of an empty voxel) everything works as it should.

But thanks again for the suggestion!