[vtkusers] Solid Voxelization with VTK

[Bill Lorensen]

No copyright issuses. Use however you wish.

Enjoy,

Bill

On Wed, Jun 13, 2018, 6:31 PM Berti Krüger <berti_krueger at hotmail.com>
wrote:

> 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:*
>
> [image: 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.*
>
> *Regard**s*
>
> *Berti*
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
>
> Search the list archives at: http://markmail.org/search/?q=vtkusers
>
> Follow this link to subscribe/unsubscribe:
> https://public.kitware.com/mailman/listinfo/vtkusers
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/vtkusers/attachments/20180613/0524e3de/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: left_amygdala.png
Type: image/png
Size: 146782 bytes
Desc: not available
URL: <https://public.kitware.com/pipermail/vtkusers/attachments/20180613/0524e3de/attachment-0001.png>