[vtkusers] Solid Voxelization with VTK

[Berti Krüger]

Am 11.06.2018 um 18:34 schrieb David Gobbi:

Hi Berti,


Hi David.

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?

Yes. And i now understand why the visualization i get from VTK is not what i expected it to be.
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.

Exactly. Thank you very much for this great explanation. I think some of this should really be added to the examples since i seem not to be the only one who seems to be confused by the terminology of VTK. While these abstractions that vtk does are powerful because they can be applied independly of the dimensionality / topology of the underlying data it would be good to have some notes for common entities in 3D graphics because i got similar problems with 3D meshes (like what means vtkCell in this regard, etc.)


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.

Thanks again. I should have seen this by myself, but i coded it in a hurry to just give an example to show to you where my problems lie.

But originally i meant slower regarding the visualization not the loop. Sorry i didn't make myself clearer. But after your explanation it is now obvious to me why my visualization using the glyphs is a bit slower (not by much) than the original visualization with the vtkThreshold approach. It simply does draw more primitives since my definition of a voxel results in more geometry for the graphics card.


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

Thanks for the suggestion. But since you approved my approach with the glyphs and speedwise i am quite happy with it (i only wondered why it was a bit slower than the vtkThreshold approach thinking i was doing something wrong) i think i will stay with it. However if i am getting problems with the rendering speed with larger volumes later on it might be good to keep this in mind. So thanks again.


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.

Ok. I see. These constraints are to strict for me since i don't want to have the voxel image resolution be constrained by the geometry of the mesh.

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.


Thanks again. While i might not be using the vtkPolyDataToImageStencil approach i think i will still need to check if the surface mesh is watertight. So this is really useful.

So after all i will try the Slicer3D code Csaba has pointed me too. If nothing works as expected i have to write my own.
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


While i am not totally fixed at this voxel approach and still be open to alternatives e.g. working with the mesh directly, from a first consideration the discretization of the space into voxels (image samples / basis functions) makes some operations and reasoning easier and hopefully less computationally demanding. At least i think so but i could err on it.

I have taken this image from a thread of the slicer newsgroup (it is attributed to Bill Lorensen) so i hope i didn't do any copyright infringement to give an example:

[amygdala]

The caption says that this should be the left amygdala. I am not an expert in brain anatomy but i think only the lower red ball-like part of this image is the actual amygdala and the rest is part of the striatum, caudate nucleus, lenticular nucleus and more. But anyway it does not matter here at all.

Instead of just using 0 and 1 for the voxel image samples (unset or set, meaning inside or outside of the mesh) i have more options to use (as shown by the different colors in the image) since i use vtkUnsignedChar (vtkBool is not so useful anyway). So i can use these values to label different regions inside of the mesh. So i might easily slice certain regions away or use morphological or general image analysis operations. I can use these labels as landmarks so algorithms know easily in which part they actually operate. The voxel values can also be used to penalize / forbid or encourage certain paths while routing through the mesh (using them as a kind of weight). For multiple paths the voxel values can be used as marks that there is already some other path going through so there might not be a need for some kind of collision detection with AABBs, kd-Trees, Octrees and so on. The instruments i target have a limited resolution / precision anyway so maybe there might also be some mapping from the visualization / representation to the physical reality which could somehow be exploited.

Without a deep analysis these are the first thoughts that came to my mind at the moment. There might be more. But the voxel representation is not yet set in stone. Of course you could all do this on the original surface mesh. But for me it was the first approach to get these requirements (which can still change) in one relatively simple uniform representation. Maybe all of this will change and maybe there might also be much better ways. But first i think you have to start somewhere and see if your ideas work. And i have the hope that vtk helps me trying out ideas without much additional boilerplate code.


Thanks for all.

Regards

Berti
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/vtkusers/attachments/20180614/591ddde9/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: left_amygdala.png
Type: image/png
Size: 146782 bytes
Desc: left_amygdala.png
URL: <https://public.kitware.com/pipermail/vtkusers/attachments/20180614/591ddde9/attachment-0001.png>