[vtkusers] Solid Voxelization with VTK

[David Gobbi]

Hi Berti,

Currently you have two distinct parts to your pipeline:

1) generate the image (via vtkPolyDataToImageStencil)

2) visualize the result (via vtkThreshold)


Let's start with vtkThreshold, since it seems to be the filter whose
behavior is furthest from your expectations.

The vtkThreshold filter is designed to extract cells.  In your case, each
"cell" that it extracts is a cube.  But it's crucial to understand how a
vtk "cell" is defined: each corner of the cell is a point in the data set.
In this case, each of these points is one of the image samples, i.e. each
point is what image analysis folks commonly refer to as a "voxel".

So the cubes that vtkThreshold produces are not the image voxels.  Instead,
each cube is in fact a finite element whose corners are positioned at the
image voxels.  To make things even more confusing, VTK calls each of these
elements a vtkVoxel.  But a vtkVoxel is not an image voxel.  A vtkVoxel is
a cube that has an image voxel at each corner.  You with me so far?

What vtkThreshold output, by default, is every cube for which every corner
of the cube is above the threshold.  You can use
threshold->SetAllScalars(0) to make it output every cube for which _any_
corner is above threshold.  This makes it much more inclusive.  But this is
probably still not what you want.


Your approach of using glyphs is good, since what you really want to do is
render a small cube for every point in the image.  Your loop will run
faster if you call voxelImage->GetScalarPointer(x, y, z) just once instead
of calling it for every iteration.  It returns a pointer, so all you have
to do is increment the pointer at every iteration.  Or you could replace
the loop with the vtkThresholdPoints filter.


An alternative that might be faster (I say "might" because I'm not sure) is
to use the vtkDiscreteMarchingCubes class.


Regarding vtkPolyDataToImageStencil, this class selects all voxels whose
center points are within volume enclosed by the surface.  It is not enough
for the surface to merely pass through the voxel, this class wasn't
designed to do that.   So it is necessary for the voxel spacing to be
smaller than the thinnest part of the model.

Also, vtkPolyDataToImageStencil expects the input data to be truly closed
(it isn't enough if the surface merely "looks" like it is watertight).
Otherwise, it can produce artifacts like the ones you have shown.  You can
test whether a surface is closed with the vtkFeatureEdges class.  If you
update vtkFeatureEdges for your original polydata, it should produce no
output cells when used like this:

vtkNew<vtkFeatureEdges> features;
features->BoundaryEdgesOn();
features->NonManifoldEdgesOn();
features->FeatureEdgesOff();
features->SetInputData(....);
features->Update();

If features->GetOutput()->GetNumberOfCells() is not zero, then the polydata
is not truly watertight.

I have a question for you: if your goal is to find paths through the
object, then why convert it into voxels?  Why not work directly with the
original surface data?

 - David
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/vtkusers/attachments/20180611/d6fb7afb/attachment.html>