vtkMarchingCubes / vtkDiscreteMarchingCubes does not produce a closed mesh surface

Support
1

Hello Everyone.

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

voxelgrid
voxelgrid757Ă—768 45.6 KB

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->SetInputData(voxelImage);
surface->ComputeNormalsOn();
surface->SetValue(0, 1.0);
//surface->GenerateValues(1, 1, 1);
surface->Update();

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:

voxel_cube
voxel_cube749Ă—755 60.9 KB

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

However, this is what i get:

marching_cube_result
marching_cube_result644Ă—598 8.25 KB

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.

2

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

1 Like
3

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?

Regards.

4

Hello again Andras.

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

Here is what i did:

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

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

And that did the trick!

5

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

6

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!