Solid Voxelization with VTK

classic Classic list List threaded Threaded
33 messages Options
12
Reply | Threaded
Open this post in threaded view
|

Solid Voxelization with VTK

bertikrueger
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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://vtk.org/mailman/listinfo/vtkusers
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

David Gobbi
Hi Berti,

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 - David



On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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://vtk.org/mailman/listinfo/vtkusers
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger

Hi David,


thank you very much for your help.


I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:


Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).


For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/



I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.


Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?


Thank you very much in advance.



Regards


Berti









Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK
 
Hi Berti,

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:
https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/



If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 - David



On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

CMakeLists.txt (464 bytes) Download Attachment
PolyDataToImageDataVisualization.cxx (9K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

Csaba Pinter

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers <[hidden email]> On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi <[hidden email]>
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 

Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

lorensen.github.io

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

David Gobbi
In reply to this post by bertikrueger
Hi Berti,

The correct pad filter to use for this is "vtkImageConstantPad", not "WrapPad", but padding really isn't necessary.  Neither is the use of cell scalars, it's more appropriate to use point scalars when working with image data since each voxel is represented as a point scalar.

The threshold should be set like this, for a threshold halfway between your "inval" and "outval":

  selector->ThresholdByUpper(0.5*(startLabel + endLabel));

The lookup table wasn't correctly set. Calling SetRange() does not set the number of entries in the lookup table (by default, the lookup table has 256 entries). 

  lut->SetTableRange(0, 1);
  lut->SetNumberOfColors(2);

Also, your color for "1" was transparent, whereas it has to be opaque in order to properly see a scalar-colored cone.



I stripped the visualization part of your pipeline down to the following.  No padding is done, the thresholding uses point scalars, no shifting, and no lookup table (instead, I call ScalarVisibilityOff() and set the color via the actor's property object).  I replaced vtkPolyDataMapper with vtkDataSetMapper, which automatically skins the data.


  unsigned char startLabel = outval;
  unsigned char endLabel = inval;

  vtkNew<vtkThreshold> selector;
  // the call to SetInputArrayToProcess() isn't necessary, the default is fine
  //selector->SetInputArrayToProcess(0, 0, 0,
  //                                 vtkDataObject::FIELD_ASSOCIATION_POINTS,
  //                                 vtkDataSetAttributes::SCALARS);
  selector->SetInputData(voxelImage);
  selector->ThresholdByUpper(0.5*(startLabel + endLabel));
  selector->Update();
    
  vtkNew<vtkDataSetMapper> mapper;
  mapper->SetInputConnection(selector->GetOutputPort());
  mapper->ScalarVisibilityOff();

  vtkNew<vtkActor> actor;
  actor->SetMapper(mapper);
  actor->GetProperty()->SetColor(0.0, 1.0, 0.0);


An image of the resulting cone is attached.

 - David


On Tue, Jun 5, 2018 at 1:26 AM, Berti Krüger <[hidden email]> wrote:

Hi David,


thank you very much for your help.


I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:


Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).


For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/



I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.


Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?


Thank you very much in advance.



Regards


Berti


If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:



Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.





Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK
 
Hi Berti,

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:
https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/
If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:



If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 - David



On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

SolidCone.png (22K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger

Hi David,

thank you very much for all your help and for all the effort.

I am sorry, but i am an absolute beginner with VTK and for me it is really hard to find informations / documentation / explanations beyond the doxygen class documentation about the advanced visualization/manipulation operations in VTK (the VTK User Manual sadly only covers the very basics).

So thanks again for all your clarifications and all your help so far.

I replaced the visualization part of my pipeline with your code which makes it so much clearer, shorter and better in terms of the visualization. And it is also easy enough so that even i get it :-)

Unfortunately when i extend the simple cone source voxelization example i attached the last time by using general STL meshes via the vtkSTLReader, i get missing voxels in the visualization. To be sure that the problem is indeed in the visualization part and not in the voxelization part i used binvox, a mature voxelization program (http://www.patrickmin.com/binvox/ ) to do the voxelization and the code you suggested for the visualization of the vtkImageData via VTK.


Using this setup the Utah Teapot looks like this (the handle and other thin parts of the mesh are missing in the visualization):

Teapot
        missing voxels


For surface-only-voxelized meshes there even more voxels missing (the example here is the Stanford Bunny which has only been surface voxelized):

bunny surface
        voxelization visualization

What i don't get is that the visualization code is as simple as it gets:

    vtkNew<vtkThreshold> selector;
    selector->SetInputData(voxelImage);
    selector->ThresholdByUpper(1.0);
    selector->Update();
     
    vtkNew<vtkDataSetMapper> mapper;
    mapper->SetInputConnection(selector->GetOutputPort());
    mapper->ScalarVisibilityOff();
    mapper->StaticOn();   

When i understand it correctly, then we simply visualize every cell (voxel) whose threshold is equal or greater than 1.0 which means we visualize (filter) all voxels which are set (in my vtkImageData voxels which are set are 1, voxels which are not set 0) and the vtkDataSetMapper maps these to graphic primitives and this does the rest of the visualization. Not much to do wrong here and that basically should do the trick but it doesn't.

Any idea why this does not work correctly all the time and has problems with thin parts?


I therefore tried another approach by using the vtkGlyph3DMapper based on the example code from

https://www.vtk.org/Wiki/VTK/Examples/Cxx/Visualization/Glyph3DMapper :


    vtkNew<vtkPoints> points;
 
    for (int x = 0; x < dims[0]; x++)
    {
      for (int y = 0; y < dims[1]; y++)
      {
        for (int z = 0; z < dims[2]; z++)
        {
          char* pixel = static_cast<char*>(voxelImage->GetScalarPointer(x, y, z));
         
          if(pixel[0] == 1) // if the voxel is set (1) in the voxel image
            points->InsertNextPoint(x, y, z);
        }
      }
    }
   
    vtkNew<vtkPolyData> polydata;
    polydata->SetPoints(points);
 
    vtkNew<vtkPolyData> glyph;
 
    vtkNew<vtkCubeSource> cubeSource;
 
    vtkNew<vtkGlyph3DMapper> glyph3Dmapper;
    glyph3Dmapper->SetSourceConnection(cubeSource->GetOutputPort());
    glyph3Dmapper->SetInputData(polydata);
    glyph3Dmapper->Update();
   
    actor->SetMapper(glyph3Dmapper);


The visualization with this approach is a bit slower (probably because it does no backface culling of the backfaces of the cube polygons) but it seems to work correctly even with surface only voxelized meshes:

bunny okteapot


I attached a simple working vtk demo with both visualizations and integrated voxelized meshes (no external data files or other dependencies) as a simple cmake/vtk compileable project (PolyDataToImageDataVisualizationVoxel). Maybe it might help other vtk users in the future.

Are there any other possibilities (maybe faster) to do the visualization with vtk ?



Now having a visualization that seems to work, i tried this visualization with the suggested voxelization code from the vtk examples here https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

and while the voxelization is lighting-fast and some example meshes did get correctly voxelized, i also had a lot of meshes whose voxelizations had misclassified voxels (e.g. voxels are set where voxels should not have been set and the opposite (holes)) like with the Utah Teapot here:

teapot errorsteapot errors
      2teapot
      errors3

The problem is that for my project i have to find paths through the voxel image. So miscIassified voxels would lead to either wrong paths found (paths which do not really exist like in the first image) or maybe to a situation where no path can be found because of holes.

I attached a compileable cmake/vtk example (PolyDataToImageDataVisualizationSTL) together with some small example stl-mesh-files (teapot etc.) using the example code from the vtk example repository for showing the problem with the given approach there.


Is it possible by tweaking the vtk PolyDataToImageData example code to always get a correct voxelization or is another approach needed?


If you or anybody else has ideas / suggestions i would really be glad.
Thank for reading if you managed to get this far and have not fallen asleep and thank you very much in advance.


Regards

Berti


Am 06.06.2018 um 01:49 schrieb David Gobbi:
Hi Berti,

The correct pad filter to use for this is "vtkImageConstantPad", not "WrapPad", but padding really isn't necessary.  Neither is the use of cell scalars, it's more appropriate to use point scalars when working with image data since each voxel is represented as a point scalar.

The threshold should be set like this, for a threshold halfway between your "inval" and "outval":

  selector->ThresholdByUpper(0.5*(startLabel + endLabel));

The lookup table wasn't correctly set. Calling SetRange() does not set the number of entries in the lookup table (by default, the lookup table has 256 entries). 

  lut->SetTableRange(0, 1);
  lut->SetNumberOfColors(2);

Also, your color for "1" was transparent, whereas it has to be opaque in order to properly see a scalar-colored cone.



I stripped the visualization part of your pipeline down to the following.  No padding is done, the thresholding uses point scalars, no shifting, and no lookup table (instead, I call ScalarVisibilityOff() and set the color via the actor's property object).  I replaced vtkPolyDataMapper with vtkDataSetMapper, which automatically skins the data.


  unsigned char startLabel = outval;
  unsigned char endLabel = inval;

  vtkNew<vtkThreshold> selector;
  // the call to SetInputArrayToProcess() isn't necessary, the default is fine
  //selector->SetInputArrayToProcess(0, 0, 0,
  //                                 vtkDataObject::FIELD_ASSOCIATION_POINTS,
  //                                 vtkDataSetAttributes::SCALARS);
  selector->SetInputData(voxelImage);
  selector->ThresholdByUpper(0.5*(startLabel + endLabel));
  selector->Update();
    
  vtkNew<vtkDataSetMapper> mapper;
  mapper->SetInputConnection(selector->GetOutputPort());
  mapper->ScalarVisibilityOff();

  vtkNew<vtkActor> actor;
  actor->SetMapper(mapper);
  actor->GetProperty()->SetColor(0.0, 1.0, 0.0);


An image of the resulting cone is attached.

 - David


On Tue, Jun 5, 2018 at 1:26 AM, Berti Krüger <[hidden email]> wrote:

Hi David,


thank you very much for your help.


I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:


Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).


For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/



I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.


Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?


Thank you very much in advance.



Regards


Berti


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK
 
Hi Berti,

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:
If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 - David



On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

PolyDataToImageDataVisualizationVoxel.zip (16K) Download Attachment
PolyDataToImageDataVisualizationSTL.zip (283K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger
In reply to this post by Csaba Pinter

Hi Csaba,


thank you very much for you answer and for offering your help.


I already have looked at it and i think if no simpler option will come up, i will give it a try (hoping not to have to change too much to get it working with my pipeline since i am an absolute vtk beginner :-) ).


How is this code licensed ? (i am only doing research work at the moment, the whole project i am working on is internal and just for testing and part of some larger work, it might be open sourced one day in the very far future but could even be not at all  and it will never be used commercially).



Regards


Berti



Am 05.06.2018 um 16:14 schrieb Csaba Pinter:

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers [hidden email] On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 

Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

lorensen.github.io

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

Elvis Stansvik
In reply to this post by bertikrueger
Den mån 11 juni 2018 08:33Berti Krüger <[hidden email]> skrev:

Hi David,

thank you very much for all your help and for all the effort.

I am sorry, but i am an absolute beginner with VTK and for me it is really hard to find informations / documentation / explanations beyond the doxygen class documentation about the advanced visualization/manipulation operations in VTK (the VTK User Manual sadly only covers the very basics).

So thanks again for all your clarifications and all your help so far.

I replaced the visualization part of my pipeline with your code which makes it so much clearer, shorter and better in terms of the visualization. And it is also easy enough so that even i get it :-)

Unfortunately when i extend the simple cone source voxelization example i attached the last time by using general STL meshes via the vtkSTLReader, i get missing voxels in the visualization. To be sure that the problem is indeed in the visualization part and not in the voxelization part i used binvox, a mature voxelization program (http://www.patrickmin.com/binvox/ ) to do the voxelization and the code you suggested for the visualization of the vtkImageData via VTK.


Using this setup the Utah Teapot looks like this (the handle and other thin parts of the mesh are missing in the visualization):


For surface-only-voxelized meshes there even more voxels missing (the example here is the Stanford Bunny which has only been surface voxelized):

What i don't get is that the visualization code is as simple as it gets:

    vtkNew<vtkThreshold> selector;
    selector->SetInputData(voxelImage);
    selector->ThresholdByUpper(1.0);
    selector->Update();
     
    vtkNew<vtkDataSetMapper> mapper;
    mapper->SetInputConnection(selector->GetOutputPort());
    mapper->ScalarVisibilityOff();
    mapper->StaticOn();   

When i understand it correctly, then we simply visualize every cell (voxel) whose threshold is equal or greater than 1.0 which means we visualize (filter) all voxels which are set (in my vtkImageData voxels which are set are 1, voxels which are not set 0) and the vtkDataSetMapper maps these to graphic primitives and this does the rest of the visualization. Not much to do wrong here and that basically should do the trick but it doesn't.

Any idea why this does not work correctly all the time and has problems with thin parts?

I'm just jumping in from the side here, and I'm by no means a VTK expert. But I think what you're seeing is simply due to how VTK does volume rendering. It interpolates the given data points, which means you need at least two data points to get a voxel.

I'm guessing that at the resolution you're doing the voxelization here, your teapot handle is simply too thin to result in even a single voxel.

Think of your data points as a grid of fenceposts. You'll get a voxel in the middle of each 2x2x2 cluster of posts.

Someone correct me if I'm wrong here.

HTH,
Elvis



I therefore tried another approach by using the vtkGlyph3DMapper based on the example code from

https://www.vtk.org/Wiki/VTK/Examples/Cxx/Visualization/Glyph3DMapper :


    vtkNew<vtkPoints> points;
 
    for (int x = 0; x < dims[0]; x++)
    {
      for (int y = 0; y < dims[1]; y++)
      {
        for (int z = 0; z < dims[2]; z++)
        {
          char* pixel = static_cast<char*>(voxelImage->GetScalarPointer(x, y, z));
         
          if(pixel[0] == 1) // if the voxel is set (1) in the voxel image
            points->InsertNextPoint(x, y, z);
        }
      }
    }
   
    vtkNew<vtkPolyData> polydata;
    polydata->SetPoints(points);
 
    vtkNew<vtkPolyData> glyph;
 
    vtkNew<vtkCubeSource> cubeSource;
 
    vtkNew<vtkGlyph3DMapper> glyph3Dmapper;
    glyph3Dmapper->SetSourceConnection(cubeSource->GetOutputPort());
    glyph3Dmapper->SetInputData(polydata);
    glyph3Dmapper->Update();
   
    actor->SetMapper(glyph3Dmapper);


The visualization with this approach is a bit slower (probably because it does no backface culling of the backfaces of the cube polygons) but it seems to work correctly even with surface only voxelized meshes:




I attached a simple working vtk demo with both visualizations and integrated voxelized meshes (no external data files or other dependencies) as a simple cmake/vtk compileable project (PolyDataToImageDataVisualizationVoxel). Maybe it might help other vtk users in the future.

Are there any other possibilities (maybe faster) to do the visualization with vtk ?



Now having a visualization that seems to work, i tried this visualization with the suggested voxelization code from the vtk examples here https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

and while the voxelization is lighting-fast and some example meshes did get correctly voxelized, i also had a lot of meshes whose voxelizations had misclassified voxels (e.g. voxels are set where voxels should not have been set and the opposite (holes)) like with the Utah Teapot here:



The problem is that for my project i have to find paths through the voxel image. So miscIassified voxels would lead to either wrong paths found (paths which do not really exist like in the first image) or maybe to a situation where no path can be found because of holes.

I attached a compileable cmake/vtk example (PolyDataToImageDataVisualizationSTL) together with some small example stl-mesh-files (teapot etc.) using the example code from the vtk example repository for showing the problem with the given approach there.


Is it possible by tweaking the vtk PolyDataToImageData example code to always get a correct voxelization or is another approach needed?


If you or anybody else has ideas / suggestions i would really be glad.
Thank for reading if you managed to get this far and have not fallen asleep and thank you very much in advance.


Regards

Berti


Am 06.06.2018 um 01:49 schrieb David Gobbi:
Hi Berti,

The correct pad filter to use for this is "vtkImageConstantPad", not "WrapPad", but padding really isn't necessary.  Neither is the use of cell scalars, it's more appropriate to use point scalars when working with image data since each voxel is represented as a point scalar.

The threshold should be set like this, for a threshold halfway between your "inval" and "outval":

  selector->ThresholdByUpper(0.5*(startLabel + endLabel));

The lookup table wasn't correctly set. Calling SetRange() does not set the number of entries in the lookup table (by default, the lookup table has 256 entries). 

  lut->SetTableRange(0, 1);
  lut->SetNumberOfColors(2);

Also, your color for "1" was transparent, whereas it has to be opaque in order to properly see a scalar-colored cone.



I stripped the visualization part of your pipeline down to the following.  No padding is done, the thresholding uses point scalars, no shifting, and no lookup table (instead, I call ScalarVisibilityOff() and set the color via the actor's property object).  I replaced vtkPolyDataMapper with vtkDataSetMapper, which automatically skins the data.


  unsigned char startLabel = outval;
  unsigned char endLabel = inval;

  vtkNew<vtkThreshold> selector;
  // the call to SetInputArrayToProcess() isn't necessary, the default is fine
  //selector->SetInputArrayToProcess(0, 0, 0,
  //                                 vtkDataObject::FIELD_ASSOCIATION_POINTS,
  //                                 vtkDataSetAttributes::SCALARS);
  selector->SetInputData(voxelImage);
  selector->ThresholdByUpper(0.5*(startLabel + endLabel));
  selector->Update();
    
  vtkNew<vtkDataSetMapper> mapper;
  mapper->SetInputConnection(selector->GetOutputPort());
  mapper->ScalarVisibilityOff();

  vtkNew<vtkActor> actor;
  actor->SetMapper(mapper);
  actor->GetProperty()->SetColor(0.0, 1.0, 0.0);


An image of the resulting cone is attached.

 - David


On Tue, Jun 5, 2018 at 1:26 AM, Berti Krüger <[hidden email]> wrote:

Hi David,


thank you very much for your help.


I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:



I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).


For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/



I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.


Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?


Thank you very much in advance.



Regards


Berti


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK
 
Hi Berti,

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:
If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 - David



On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

_______________________________________________
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

screenshot.jpg (83K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

Csaba Pinter
In reply to this post by bertikrueger

Hi Berti,

 

That code is under BSD-style license, so you’re free to do basically anything with it. This implementation is already pretty simple; just a concatenation of several VTK filters, so you can just give it a try by copy-pasting the relevant part of that function into your code and see if it works. If you’re interested in using the vtkSegmentation library, let me know though!

 

Best,

csaba

 

From: Berti Krüger <[hidden email]>
Sent: Monday, June 11, 2018 02:46
To: Csaba Pinter <[hidden email]>; David Gobbi <[hidden email]>
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Csaba,

 

thank you very much for you answer and for offering your help.

 

I already have looked at it and i think if no simpler option will come up, i will give it a try (hoping not to have to change too much to get it working with my pipeline since i am an absolute vtk beginner :-) ).

 

How is this code licensed ? (i am only doing research work at the moment, the whole project i am working on is internal and just for testing and part of some larger work, it might be open sourced one day in the very far future but could even be not at all  and it will never be used commercially).

 

 

Regards

 

Berti

 

 

Am 05.06.2018 um 16:14 schrieb Csaba Pinter:

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers [hidden email] On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 

Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

lorensen.github.io

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

David Gobbi
In reply to this post by bertikrueger
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



_______________________________________________
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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger
In reply to this post by Elvis Stansvik

Hi Elvis,

thanks for jumping in. As you can see in David's response you were right or at
least close with your assumptions.

Regards

Berti


Am 11.06.2018 um 08:50 schrieb Elvis Stansvik:
Den mån 11 juni 2018 08:33Berti Krüger <[hidden email]> skrev:

Hi David,

thank you very much for all your help and for all the effort.

I am sorry, but i am an absolute beginner with VTK and for me it is really hard to find informations / documentation / explanations beyond the doxygen class documentation about the advanced visualization/manipulation operations in VTK (the VTK User Manual sadly only covers the very basics).

So thanks again for all your clarifications and all your help so far.

I replaced the visualization part of my pipeline with your code which makes it so much clearer, shorter and better in terms of the visualization. And it is also easy enough so that even i get it :-)

Unfortunately when i extend the simple cone source voxelization example i attached the last time by using general STL meshes via the vtkSTLReader, i get missing voxels in the visualization. To be sure that the problem is indeed in the visualization part and not in the voxelization part i used binvox, a mature voxelization program (http://www.patrickmin.com/binvox/ ) to do the voxelization and the code you suggested for the visualization of the vtkImageData via VTK.


Using this setup the Utah Teapot looks like this (the handle and other thin parts of the mesh are missing in the visualization):


For surface-only-voxelized meshes there even more voxels missing (the example here is the Stanford Bunny which has only been surface voxelized):

What i don't get is that the visualization code is as simple as it gets:

    vtkNew<vtkThreshold> selector;
    selector->SetInputData(voxelImage);
    selector->ThresholdByUpper(1.0);
    selector->Update();
     
    vtkNew<vtkDataSetMapper> mapper;
    mapper->SetInputConnection(selector->GetOutputPort());
    mapper->ScalarVisibilityOff();
    mapper->StaticOn();   

When i understand it correctly, then we simply visualize every cell (voxel) whose threshold is equal or greater than 1.0 which means we visualize (filter) all voxels which are set (in my vtkImageData voxels which are set are 1, voxels which are not set 0) and the vtkDataSetMapper maps these to graphic primitives and this does the rest of the visualization. Not much to do wrong here and that basically should do the trick but it doesn't.

Any idea why this does not work correctly all the time and has problems with thin parts?

I'm just jumping in from the side here, and I'm by no means a VTK expert. But I think what you're seeing is simply due to how VTK does volume rendering. It interpolates the given data points, which means you need at least two data points to get a voxel.

I'm guessing that at the resolution you're doing the voxelization here, your teapot handle is simply too thin to result in even a single voxel.

Think of your data points as a grid of fenceposts. You'll get a voxel in the middle of each 2x2x2 cluster of posts.

Someone correct me if I'm wrong here.

HTH,
Elvis



I therefore tried another approach by using the vtkGlyph3DMapper based on the example code from

https://www.vtk.org/Wiki/VTK/Examples/Cxx/Visualization/Glyph3DMapper :


    vtkNew<vtkPoints> points;
 
    for (int x = 0; x < dims[0]; x++)
    {
      for (int y = 0; y < dims[1]; y++)
      {
        for (int z = 0; z < dims[2]; z++)
        {
          char* pixel = static_cast<char*>(voxelImage->GetScalarPointer(x, y, z));
         
          if(pixel[0] == 1) // if the voxel is set (1) in the voxel image
            points->InsertNextPoint(x, y, z);
        }
      }
    }
   
    vtkNew<vtkPolyData> polydata;
    polydata->SetPoints(points);
 
    vtkNew<vtkPolyData> glyph;
 
    vtkNew<vtkCubeSource> cubeSource;
 
    vtkNew<vtkGlyph3DMapper> glyph3Dmapper;
    glyph3Dmapper->SetSourceConnection(cubeSource->GetOutputPort());
    glyph3Dmapper->SetInputData(polydata);
    glyph3Dmapper->Update();
   
    actor->SetMapper(glyph3Dmapper);


The visualization with this approach is a bit slower (probably because it does no backface culling of the backfaces of the cube polygons) but it seems to work correctly even with surface only voxelized meshes:




I attached a simple working vtk demo with both visualizations and integrated voxelized meshes (no external data files or other dependencies) as a simple cmake/vtk compileable project (PolyDataToImageDataVisualizationVoxel). Maybe it might help other vtk users in the future.

Are there any other possibilities (maybe faster) to do the visualization with vtk ?



Now having a visualization that seems to work, i tried this visualization with the suggested voxelization code from the vtk examples here https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

and while the voxelization is lighting-fast and some example meshes did get correctly voxelized, i also had a lot of meshes whose voxelizations had misclassified voxels (e.g. voxels are set where voxels should not have been set and the opposite (holes)) like with the Utah Teapot here:



The problem is that for my project i have to find paths through the voxel image. So miscIassified voxels would lead to either wrong paths found (paths which do not really exist like in the first image) or maybe to a situation where no path can be found because of holes.

I attached a compileable cmake/vtk example (PolyDataToImageDataVisualizationSTL) together with some small example stl-mesh-files (teapot etc.) using the example code from the vtk example repository for showing the problem with the given approach there.


Is it possible by tweaking the vtk PolyDataToImageData example code to always get a correct voxelization or is another approach needed?


If you or anybody else has ideas / suggestions i would really be glad.
Thank for reading if you managed to get this far and have not fallen asleep and thank you very much in advance.


Regards

Berti


Am 06.06.2018 um 01:49 schrieb David Gobbi:
Hi Berti,

The correct pad filter to use for this is "vtkImageConstantPad", not "WrapPad", but padding really isn't necessary.  Neither is the use of cell scalars, it's more appropriate to use point scalars when working with image data since each voxel is represented as a point scalar.

The threshold should be set like this, for a threshold halfway between your "inval" and "outval":

  selector->ThresholdByUpper(0.5*(startLabel + endLabel));

The lookup table wasn't correctly set. Calling SetRange() does not set the number of entries in the lookup table (by default, the lookup table has 256 entries). 

  lut->SetTableRange(0, 1);
  lut->SetNumberOfColors(2);

Also, your color for "1" was transparent, whereas it has to be opaque in order to properly see a scalar-colored cone.



I stripped the visualization part of your pipeline down to the following.  No padding is done, the thresholding uses point scalars, no shifting, and no lookup table (instead, I call ScalarVisibilityOff() and set the color via the actor's property object).  I replaced vtkPolyDataMapper with vtkDataSetMapper, which automatically skins the data.


  unsigned char startLabel = outval;
  unsigned char endLabel = inval;

  vtkNew<vtkThreshold> selector;
  // the call to SetInputArrayToProcess() isn't necessary, the default is fine
  //selector->SetInputArrayToProcess(0, 0, 0,
  //                                 vtkDataObject::FIELD_ASSOCIATION_POINTS,
  //                                 vtkDataSetAttributes::SCALARS);
  selector->SetInputData(voxelImage);
  selector->ThresholdByUpper(0.5*(startLabel + endLabel));
  selector->Update();
    
  vtkNew<vtkDataSetMapper> mapper;
  mapper->SetInputConnection(selector->GetOutputPort());
  mapper->ScalarVisibilityOff();

  vtkNew<vtkActor> actor;
  actor->SetMapper(mapper);
  actor->GetProperty()->SetColor(0.0, 1.0, 0.0);


An image of the resulting cone is attached.

 - David


On Tue, Jun 5, 2018 at 1:26 AM, Berti Krüger <[hidden email]> wrote:

Hi David,


thank you very much for your help.


I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:



I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).


For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/



I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.


Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?


Thank you very much in advance.



Regards


Berti


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK
 
Hi Berti,

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:
If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 - David



On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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


_______________________________________________
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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger
In reply to this post by Csaba Pinter

Hi Csaba,


thanks again. This is just perfect. I will try it.


Regards


Berti



Am 11.06.2018 um 14:45 schrieb Csaba Pinter:

Hi Berti,

 

That code is under BSD-style license, so you’re free to do basically anything with it. This implementation is already pretty simple; just a concatenation of several VTK filters, so you can just give it a try by copy-pasting the relevant part of that function into your code and see if it works. If you’re interested in using the vtkSegmentation library, let me know though!

 

Best,

csaba

 

From: Berti Krüger [hidden email]
Sent: Monday, June 11, 2018 02:46
To: Csaba Pinter [hidden email]; David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Csaba,

 

thank you very much for you answer and for offering your help.

 

I already have looked at it and i think if no simpler option will come up, i will give it a try (hoping not to have to change too much to get it working with my pipeline since i am an absolute vtk beginner :-) ).

 

How is this code licensed ? (i am only doing research work at the moment, the whole project i am working on is internal and just for testing and part of some larger work, it might be open sourced one day in the very far future but could even be not at all  and it will never be used commercially).

 

 

Regards

 

Berti

 

 

Am 05.06.2018 um 16:14 schrieb Csaba Pinter:

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers [hidden email] On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 

Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

lorensen.github.io

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger
In reply to this post by David Gobbi

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

_______________________________________________
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
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[hidden email]> 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:



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
_______________________________________________
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

_______________________________________________
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

left_amygdala.png (196K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

David Gobbi
In reply to this post by bertikrueger
Hi Berti, I'm glad that I could help.

 - David


On Wed, Jun 13, 2018 at 7:30 PM, Berti Krüger <[hidden email]> wrote:
Thanks for all.

Regards

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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

Elvis Stansvik
In reply to this post by bertikrueger
2018-06-14 3:21 GMT+02:00 Berti Krüger <[hidden email]>:

Hi Elvis,

thanks for jumping in. As you can see in David's response you were right or at
least close with your assumptions.

Yep, you can always trust David to have all the nitty gritty details. He's the friendly VTK ninja around here :)

Elvis


Regards

Berti


Am 11.06.2018 um 08:50 schrieb Elvis Stansvik:
Den mån 11 juni 2018 08:33Berti Krüger <[hidden email]> skrev:

Hi David,

thank you very much for all your help and for all the effort.

I am sorry, but i am an absolute beginner with VTK and for me it is really hard to find informations / documentation / explanations beyond the doxygen class documentation about the advanced visualization/manipulation operations in VTK (the VTK User Manual sadly only covers the very basics).

So thanks again for all your clarifications and all your help so far.

I replaced the visualization part of my pipeline with your code which makes it so much clearer, shorter and better in terms of the visualization. And it is also easy enough so that even i get it :-)

Unfortunately when i extend the simple cone source voxelization example i attached the last time by using general STL meshes via the vtkSTLReader, i get missing voxels in the visualization. To be sure that the problem is indeed in the visualization part and not in the voxelization part i used binvox, a mature voxelization program (http://www.patrickmin.com/binvox/ ) to do the voxelization and the code you suggested for the visualization of the vtkImageData via VTK.


Using this setup the Utah Teapot looks like this (the handle and other thin parts of the mesh are missing in the visualization):

Teapot missing voxels


For surface-only-voxelized meshes there even more voxels missing (the example here is the Stanford Bunny which has only been surface voxelized):

bunny surface voxelization visualization

What i don't get is that the visualization code is as simple as it gets:

    vtkNew<vtkThreshold> selector;
    selector->SetInputData(voxelImage);
    selector->ThresholdByUpper(1.0);
    selector->Update();
     
    vtkNew<vtkDataSetMapper> mapper;
    mapper->SetInputConnection(selector->GetOutputPort());
    mapper->ScalarVisibilityOff();
    mapper->StaticOn();   

When i understand it correctly, then we simply visualize every cell (voxel) whose threshold is equal or greater than 1.0 which means we visualize (filter) all voxels which are set (in my vtkImageData voxels which are set are 1, voxels which are not set 0) and the vtkDataSetMapper maps these to graphic primitives and this does the rest of the visualization. Not much to do wrong here and that basically should do the trick but it doesn't.

Any idea why this does not work correctly all the time and has problems with thin parts?

I'm just jumping in from the side here, and I'm by no means a VTK expert. But I think what you're seeing is simply due to how VTK does volume rendering. It interpolates the given data points, which means you need at least two data points to get a voxel.

I'm guessing that at the resolution you're doing the voxelization here, your teapot handle is simply too thin to result in even a single voxel.

Think of your data points as a grid of fenceposts. You'll get a voxel in the middle of each 2x2x2 cluster of posts.

Someone correct me if I'm wrong here.

HTH,
Elvis



I therefore tried another approach by using the vtkGlyph3DMapper based on the example code from

https://www.vtk.org/Wiki/VTK/Examples/Cxx/Visualization/Glyph3DMapper :


    vtkNew<vtkPoints> points;
 
    for (int x = 0; x < dims[0]; x++)
    {
      for (int y = 0; y < dims[1]; y++)
      {
        for (int z = 0; z < dims[2]; z++)
        {
          char* pixel = static_cast<char*>(voxelImage->GetScalarPointer(x, y, z));
         
          if(pixel[0] == 1) // if the voxel is set (1) in the voxel image
            points->InsertNextPoint(x, y, z);
        }
      }
    }
   
    vtkNew<vtkPolyData> polydata;
    polydata->SetPoints(points);
 
    vtkNew<vtkPolyData> glyph;
 
    vtkNew<vtkCubeSource> cubeSource;
 
    vtkNew<vtkGlyph3DMapper> glyph3Dmapper;
    glyph3Dmapper->SetSourceConnection(cubeSource->GetOutputPort());
    glyph3Dmapper->SetInputData(polydata);
    glyph3Dmapper->Update();
   
    actor->SetMapper(glyph3Dmapper);


The visualization with this approach is a bit slower (probably because it does no backface culling of the backfaces of the cube polygons) but it seems to work correctly even with surface only voxelized meshes:

bunny okteapot


I attached a simple working vtk demo with both visualizations and integrated voxelized meshes (no external data files or other dependencies) as a simple cmake/vtk compileable project (PolyDataToImageDataVisualizationVoxel). Maybe it might help other vtk users in the future.

Are there any other possibilities (maybe faster) to do the visualization with vtk ?



Now having a visualization that seems to work, i tried this visualization with the suggested voxelization code from the vtk examples here https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

and while the voxelization is lighting-fast and some example meshes did get correctly voxelized, i also had a lot of meshes whose voxelizations had misclassified voxels (e.g. voxels are set where voxels should not have been set and the opposite (holes)) like with the Utah Teapot here:

teapot errorsteapot errors 2teapot errors3

The problem is that for my project i have to find paths through the voxel image. So miscIassified voxels would lead to either wrong paths found (paths which do not really exist like in the first image) or maybe to a situation where no path can be found because of holes.

I attached a compileable cmake/vtk example (PolyDataToImageDataVisualizationSTL) together with some small example stl-mesh-files (teapot etc.) using the example code from the vtk example repository for showing the problem with the given approach there.


Is it possible by tweaking the vtk PolyDataToImageData example code to always get a correct voxelization or is another approach needed?


If you or anybody else has ideas / suggestions i would really be glad.
Thank for reading if you managed to get this far and have not fallen asleep and thank you very much in advance.


Regards

Berti


Am 06.06.2018 um 01:49 schrieb David Gobbi:
Hi Berti,

The correct pad filter to use for this is "vtkImageConstantPad", not "WrapPad", but padding really isn't necessary.  Neither is the use of cell scalars, it's more appropriate to use point scalars when working with image data since each voxel is represented as a point scalar.

The threshold should be set like this, for a threshold halfway between your "inval" and "outval":

  selector->ThresholdByUpper(0.5*(startLabel + endLabel));

The lookup table wasn't correctly set. Calling SetRange() does not set the number of entries in the lookup table (by default, the lookup table has 256 entries). 

  lut->SetTableRange(0, 1);
  lut->SetNumberOfColors(2);

Also, your color for "1" was transparent, whereas it has to be opaque in order to properly see a scalar-colored cone.



I stripped the visualization part of your pipeline down to the following.  No padding is done, the thresholding uses point scalars, no shifting, and no lookup table (instead, I call ScalarVisibilityOff() and set the color via the actor's property object).  I replaced vtkPolyDataMapper with vtkDataSetMapper, which automatically skins the data.


  unsigned char startLabel = outval;
  unsigned char endLabel = inval;

  vtkNew<vtkThreshold> selector;
  // the call to SetInputArrayToProcess() isn't necessary, the default is fine
  //selector->SetInputArrayToProcess(0, 0, 0,
  //                                 vtkDataObject::FIELD_ASSOCIATION_POINTS,
  //                                 vtkDataSetAttributes::SCALARS);
  selector->SetInputData(voxelImage);
  selector->ThresholdByUpper(0.5*(startLabel + endLabel));
  selector->Update();
    
  vtkNew<vtkDataSetMapper> mapper;
  mapper->SetInputConnection(selector->GetOutputPort());
  mapper->ScalarVisibilityOff();

  vtkNew<vtkActor> actor;
  actor->SetMapper(mapper);
  actor->GetProperty()->SetColor(0.0, 1.0, 0.0);


An image of the resulting cone is attached.

 - David


On Tue, Jun 5, 2018 at 1:26 AM, Berti Krüger <[hidden email]> wrote:

Hi David,


thank you very much for your help.


I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:


Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).


For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/



I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.


Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?


Thank you very much in advance.



Regards


Berti


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK
 
Hi Berti,

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:
If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 - David



On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:
Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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



_______________________________________________
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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger
In reply to this post by Csaba Pinter

Hi all.


I am very sorry to revive this ancient thread. 


But i thought that maybe it could help other vtk users in the future to have a compilable and working example of the final result. If someone thinks that this code could be useful as an vtk example then feel free to modify / correct / pimp / add this to the vtk examples on github.


I adapted the 3D slicer vtk code Csaba has pointed me to and incorporated the hints from David and everything works like a charm:



The whole procedure is incredibly simple. Basically the vtkPolyData of the mesh is converted via the vtkPolyDataToImageStencil filter to a vtkImageStencil and the rest of the whole magic is only done by the vtkImageStencil filter:


// Convert stencil to image
vtkNew<vtkImageStencil> stencil;
stencil->SetInputData(binaryLabelMap);
stencil->SetStencilConnection(polyDataToImageStencil->GetOutputPort());
stencil->ReverseStencilOn();
stencil->SetBackgroundValue(1); // General foreground value is 1 (background value because of reverse stencil)

I would really like to know very roughly how the vtkImageStencil works internally without digging through all of the source code. Because it is really fast (some kind of rasterization maybe ?). If someone could drop me a line on how it works only in principle  that would be really great.


So thanks again to Bill, Csaba, David, Elvis and everyone else on this thread.



Kind regards,


Berti



PS: Example code, CMakeLists and example STL mesh are attached


Von: Csaba Pinter <[hidden email]>
Gesendet: Montag, 11. Juni 2018 12:45
An: Berti Krüger
Cc: [hidden email]
Betreff: RE: [vtkusers] Solid Voxelization with VTK
 

Hi Berti,

 

That code is under BSD-style license, so you’re free to do basically anything with it. This implementation is already pretty simple; just a concatenation of several VTK filters, so you can just give it a try by copy-pasting the relevant part of that function into your code and see if it works. If you’re interested in using the vtkSegmentation library, let me know though!

 

Best,

csaba

 

From: Berti Krüger <[hidden email]>
Sent: Monday, June 11, 2018 02:46
To: Csaba Pinter <[hidden email]>; David Gobbi <[hidden email]>
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Csaba,

 

thank you very much for you answer and for offering your help.

 

I already have looked at it and i think if no simpler option will come up, i will give it a try (hoping not to have to change too much to get it working with my pipeline since i am an absolute vtk beginner :-) ).

 

How is this code licensed ? (i am only doing research work at the moment, the whole project i am working on is internal and just for testing and part of some larger work, it might be open sourced one day in the very far future but could even be not at all  and it will never be used commercially).

 

 

Regards

 

Berti

 

 

Am 05.06.2018 um 16:14 schrieb Csaba Pinter:

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers [hidden email] On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 

Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

lorensen.github.io

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

lorensen.github.io

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

CMakeLists.txt (482 bytes) Download Attachment
PolyDataToImageDataVisualizationSTL.cxx (10K) Download Attachment
spacefighter.stl (12K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

David Gobbi
Hi Berti,

Thanks for the code, I'm sure that people will find it useful.  The vtkImageStencilData class isn't described in much detail in the VTK documentation, but I can provide a brief description.

The vtkImageStencilData object stores each raster-line of a binary image as a list of (begin, end) pairs.  That is, each raster line has zero or more non-overlapping (begin, end) pairs.  When vtkImageStencil operates on an image, it simply iterates through the (begin, end) pairs to figure out what parts of the image to copy.  One way to think about this is that if the original binary image volume is stored as O(n^3) voxels, the vtkImageStencilData stores O(n^2) values.

When a vtkImageStencilData is created, it allocates enough memory to store one (begin, end) pair per raster line, since this matches the most common use case (i.e. a mask that corresponds to a convex object). Thereafter it performs additional allocations as necessary.

As for vtkPolyDataToImageStencil, it works by cutting your STL object along each XY plane to create polygons (zero or more polygons per slice), and then it rasterizes each polygon to create the (begin, end) pairs for each raster line.


One comment about your code: it would be more efficient to use vtkImageStencilToImage instead of vtkImageStencil and vtkImageCast.  The vtkImageStencilToImage filter requires only one input (the stencil), i.e. it doesn't require 'binaryLabelMap' as input.  Also, you can set the data type of the image that it produces as output, so it doesn't have to be followed by a cast.

Cheers,
 - David

On Sat, Sep 22, 2018 at 10:29 PM Berti Krüger <[hidden email]> wrote:

Hi all.


I am very sorry to revive this ancient thread. 


But i thought that maybe it could help other vtk users in the future to have a compilable and working example of the final result. If someone thinks that this code could be useful as an vtk example then feel free to modify / correct / pimp / add this to the vtk examples on github.


I adapted the 3D slicer vtk code Csaba has pointed me to and incorporated the hints from David and everything works like a charm:



The whole procedure is incredibly simple. Basically the vtkPolyData of the mesh is converted via the vtkPolyDataToImageStencil filter to a vtkImageStencil and the rest of the whole magic is only done by the vtkImageStencil filter:


// Convert stencil to image
vtkNew<vtkImageStencil> stencil;
stencil->SetInputData(binaryLabelMap);
stencil->SetStencilConnection(polyDataToImageStencil->GetOutputPort());
stencil->ReverseStencilOn();
stencil->SetBackgroundValue(1); // General foreground value is 1 (background value because of reverse stencil)

I would really like to know very roughly how the vtkImageStencil works internally without digging through all of the source code. Because it is really fast (some kind of rasterization maybe ?). If someone could drop me a line on how it works only in principle  that would be really great.


So thanks again to Bill, Csaba, David, Elvis and everyone else on this thread.



Kind regards,


Berti



PS: Example code, CMakeLists and example STL mesh are attached


Von: Csaba Pinter <[hidden email]>
Gesendet: Montag, 11. Juni 2018 12:45
An: Berti Krüger
Cc: [hidden email]
Betreff: RE: [vtkusers] Solid Voxelization with VTK
 

Hi Berti,

 

That code is under BSD-style license, so you’re free to do basically anything with it. This implementation is already pretty simple; just a concatenation of several VTK filters, so you can just give it a try by copy-pasting the relevant part of that function into your code and see if it works. If you’re interested in using the vtkSegmentation library, let me know though!

 

Best,

csaba

 

From: Berti Krüger <[hidden email]>
Sent: Monday, June 11, 2018 02:46
To: Csaba Pinter <[hidden email]>; David Gobbi <[hidden email]>
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Csaba,

 

thank you very much for you answer and for offering your help.

 

I already have looked at it and i think if no simpler option will come up, i will give it a try (hoping not to have to change too much to get it working with my pipeline since i am an absolute vtk beginner :-) ).

 

How is this code licensed ? (i am only doing research work at the moment, the whole project i am working on is internal and just for testing and part of some larger work, it might be open sourced one day in the very far future but could even be not at all  and it will never be used commercially).

 

 

Regards

 

Berti

 

 

Am 05.06.2018 um 16:14 schrieb Csaba Pinter:

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers [hidden email] On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 


I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

image001.jpg (83K) Download Attachment
dragon.jpg (88K) Download Attachment
image001.jpg (83K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

Bill Lorensen
Looks like I should add some examples!

On Sun, Sep 23, 2018, 9:15 AM David Gobbi <[hidden email]> wrote:
Hi Berti,

Thanks for the code, I'm sure that people will find it useful.  The vtkImageStencilData class isn't described in much detail in the VTK documentation, but I can provide a brief description.

The vtkImageStencilData object stores each raster-line of a binary image as a list of (begin, end) pairs.  That is, each raster line has zero or more non-overlapping (begin, end) pairs.  When vtkImageStencil operates on an image, it simply iterates through the (begin, end) pairs to figure out what parts of the image to copy.  One way to think about this is that if the original binary image volume is stored as O(n^3) voxels, the vtkImageStencilData stores O(n^2) values.

When a vtkImageStencilData is created, it allocates enough memory to store one (begin, end) pair per raster line, since this matches the most common use case (i.e. a mask that corresponds to a convex object). Thereafter it performs additional allocations as necessary.

As for vtkPolyDataToImageStencil, it works by cutting your STL object along each XY plane to create polygons (zero or more polygons per slice), and then it rasterizes each polygon to create the (begin, end) pairs for each raster line.


One comment about your code: it would be more efficient to use vtkImageStencilToImage instead of vtkImageStencil and vtkImageCast.  The vtkImageStencilToImage filter requires only one input (the stencil), i.e. it doesn't require 'binaryLabelMap' as input.  Also, you can set the data type of the image that it produces as output, so it doesn't have to be followed by a cast.

Cheers,
 - David

On Sat, Sep 22, 2018 at 10:29 PM Berti Krüger <[hidden email]> wrote:

Hi all.


I am very sorry to revive this ancient thread. 


But i thought that maybe it could help other vtk users in the future to have a compilable and working example of the final result. If someone thinks that this code could be useful as an vtk example then feel free to modify / correct / pimp / add this to the vtk examples on github.


I adapted the 3D slicer vtk code Csaba has pointed me to and incorporated the hints from David and everything works like a charm:



The whole procedure is incredibly simple. Basically the vtkPolyData of the mesh is converted via the vtkPolyDataToImageStencil filter to a vtkImageStencil and the rest of the whole magic is only done by the vtkImageStencil filter:


// Convert stencil to image
vtkNew<vtkImageStencil> stencil;
stencil->SetInputData(binaryLabelMap);
stencil->SetStencilConnection(polyDataToImageStencil->GetOutputPort());
stencil->ReverseStencilOn();
stencil->SetBackgroundValue(1); // General foreground value is 1 (background value because of reverse stencil)

I would really like to know very roughly how the vtkImageStencil works internally without digging through all of the source code. Because it is really fast (some kind of rasterization maybe ?). If someone could drop me a line on how it works only in principle  that would be really great.


So thanks again to Bill, Csaba, David, Elvis and everyone else on this thread.



Kind regards,


Berti



PS: Example code, CMakeLists and example STL mesh are attached


Von: Csaba Pinter <[hidden email]>
Gesendet: Montag, 11. Juni 2018 12:45
An: Berti Krüger
Cc: [hidden email]
Betreff: RE: [vtkusers] Solid Voxelization with VTK
 

Hi Berti,

 

That code is under BSD-style license, so you’re free to do basically anything with it. This implementation is already pretty simple; just a concatenation of several VTK filters, so you can just give it a try by copy-pasting the relevant part of that function into your code and see if it works. If you’re interested in using the vtkSegmentation library, let me know though!

 

Best,

csaba

 

From: Berti Krüger <[hidden email]>
Sent: Monday, June 11, 2018 02:46
To: Csaba Pinter <[hidden email]>; David Gobbi <[hidden email]>
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Csaba,

 

thank you very much for you answer and for offering your help.

 

I already have looked at it and i think if no simpler option will come up, i will give it a try (hoping not to have to change too much to get it working with my pipeline since i am an absolute vtk beginner :-) ).

 

How is this code licensed ? (i am only doing research work at the moment, the whole project i am working on is internal and just for testing and part of some larger work, it might be open sourced one day in the very far future but could even be not at all  and it will never be used commercially).

 

 

Regards

 

Berti

 

 

Am 05.06.2018 um 16:14 schrieb Csaba Pinter:

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers [hidden email] On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 

Screenshot result voxelization
I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

_______________________________________________
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
Reply | Threaded
Open this post in threaded view
|

Re: Solid Voxelization with VTK

bertikrueger
In reply to this post by David Gobbi


Hi David,

thank you very much again for your great explanations (i looked at the vtkImageStencilData class description before writing the e-mail and as you have pointed out it is really missing some explanations on how it works internally) and of course your hints to improve the example code!

As you have written, i changed my code from

  // Convert stencil to image
  vtkNew<vtkImageStencil> stencil;
  stencil->SetInputData(binaryLabelMap);
  stencil->SetStencilConnection(polyDataToImageStencil->GetOutputPort());
  stencil->ReverseStencilOn();
  stencil->SetBackgroundValue(1); // General foreground value is 1 (background value because of reverse stencil)

  // Save result to output
  vtkNew<vtkImageCast> imageCast;
  imageCast->SetInputConnection(stencil->GetOutputPort());
  imageCast->SetOutputScalarTypeToUnsignedChar();
  imageCast->Update();
  binaryLabelMap->ShallowCopy(imageCast->GetOutput());



To:

  // Convert stencil to image
 vtkNew<vtkImageStencilToImage> imageStencilToImage;
 imageStencilToImage->SetInputConnection(polyDataToImageStencil->GetOutputPort());
 imageStencilToImage->SetOutsideValue(0);
 imageStencilToImage->SetInsideValue(1);
 imageStencilToImage->Update();

 // Save result to output
 binaryLabelMap->DeepCopy(imageStencilToImage->GetOutput());


which is really much shorter and better.


I also found out that i can change the vtkImage initialization code from

// fill the image with background voxels:
vtkIdType count = voxelImage->GetNumberOfPoints();
for (vtkIdType i = 0; i < count; ++i)
{
    voxelImage->GetPointData()->GetScalars()->SetTuple1(i, outval);
}

to simply:

voxelImage->GetPointData()->GetScalars()->Fill(outval);

I attached the improved version and the CMakeLists file.


So far it works really great and is quite fast:





But one thing is that i still get problems with some kind of meshes (i guess meshes which are not closed or watertight):


Alien Queen (minor problems):






Eiffel Tower (completely wrong):





Maybe this could be solved by rasterizing the mesh in different directions and doing some kind of majority voting if a voxel is set or not.

I solved these cases by implementing a raycasting solution using a generalization of the point in polygon algorithm and 

  // load STL file
  vtkNew<vtkSTLReader> reader;
  reader->SetFileName("eiffel_tower_small.stl"); 
  reader->Update();

  vtkNew<vtkPolyData> pd;
  pd->DeepCopy(reader->GetOutput());

  // compute bounds for stl mesh polydata
  double bounds[6];
  pd->GetBounds(bounds);

/*
 // Create the locator (obb tree)
  vtkNew<vtkOBBTree> obbTree;
  tree->SetDataSet(pd);
  tree->BuildLocator();
*/

  // Create the locator (modified bsp tree)
  vtkNew<vtkModifiedBSPTree> bspTree;
  bspTree->SetDataSet(pd);
  bspTree->BuildLocator();
  
  constexpr int numSamples = 32;
            
  
  //#pragma omp parallel for // not for vtkModifiedBSPTree !!
  for(int x = 0; x < imageResolutionX; x++)
  {
    for(int y = 0; y < imageResolutionY; y++)
    {
        for(int z = 0; z < imageResolutionZ; z++)
        {
    // pick an arbitrary sampling direction
    vtkVector3f ray_direction(0.0, 0.0, 0.0);
    vtkVector3f ray_direction2(0.0, 0.0, 0.0);

            // find world space position of the voxel
            constexpr int X_MIN = 0;
            constexpr int Y_MIN = 2;
            constexpr int Z_MIN = 4;
  
            vtkVector3f ray_origin( bounds[X_MIN] + spacing[0] * x,
                                    bounds[Y_MIN] + spacing[1] * y,
                                    bounds[Z_MIN] + spacing[2] * z);


    // randomly sample 3D space by sending rays in randomized directions
            unsigned int numIntersections = 0;

            int id = 0;
            
    for (int j = 0; j < numSamples; ++j)
    {
    // compute the random direction. Convert from polar to 
          // euclidean to get an even distribution
               float theta = 2 * M_PI * random(genenerator);
               float phi = acos(1 - 2 * random(generator));
            
               ray_direction[0] = sin(phi) * cos(theta);
               ray_direction[1] = sin(phi) * sin(theta);
               ray_direction[2] = cos(phi);


            // check if the voxel is inside of the mesh.
            vtkNew<vtkPoints> intersectPoints;
      
               double lineP0[3] = {ray_origin[0], ray_origin[1], ray_origin[2]};
                
               double lineP1[3] = {ray_origin[0] + ray_direction[0] * 10000.0, 
                                   ray_origin[1] + ray_direction[1] * 10000.0, 
                                   ray_origin[2] + ray_direction[2] * 10000.0};
                
                
               //obbTree->IntersectWithLine(lineP0, lineP1, intersectPoints, NULL);
                
               bspTree->IntersectWithLine(lineP0, lineP1, 0.0001, intersectPoints, NULL);       

               if (intersectPoints->GetNumberOfPoints() % 2 == 1)
                   numIntersections++;      
    }
    
              // save answer
            unsigned int index = x + y * imageResolutionX 
                                   + z * imageResolutionX * imageResolutionY;

            
            float hitPercentage = static_cast<float>(numIntersections) / static_cast<float>(numSamples);
            
if(hitPercentage >= 0.5)
                voxelImage[index] = inval;
    else
        voxelImage[index] = outval;
}
}
  }





Von: David Gobbi <[hidden email]>
Gesendet: Sonntag, 23. September 2018 16:15
An: Berti Krüger
Cc: Csaba Pinter; VTK Users
Betreff: Re: [vtkusers] Solid Voxelization with VTK
 
Hi Berti,

Thanks for the code, I'm sure that people will find it useful.  The vtkImageStencilData class isn't described in much detail in the VTK documentation, but I can provide a brief description.

The vtkImageStencilData object stores each raster-line of a binary image as a list of (begin, end) pairs.  That is, each raster line has zero or more non-overlapping (begin, end) pairs.  When vtkImageStencil operates on an image, it simply iterates through the (begin, end) pairs to figure out what parts of the image to copy.  One way to think about this is that if the original binary image volume is stored as O(n^3) voxels, the vtkImageStencilData stores O(n^2) values.

When a vtkImageStencilData is created, it allocates enough memory to store one (begin, end) pair per raster line, since this matches the most common use case (i.e. a mask that corresponds to a convex object). Thereafter it performs additional allocations as necessary.

As for vtkPolyDataToImageStencil, it works by cutting your STL object along each XY plane to create polygons (zero or more polygons per slice), and then it rasterizes each polygon to create the (begin, end) pairs for each raster line.


One comment about your code: it would be more efficient to use vtkImageStencilToImage instead of vtkImageStencil and vtkImageCast.  The vtkImageStencilToImage filter requires only one input (the stencil), i.e. it doesn't require 'binaryLabelMap' as input.  Also, you can set the data type of the image that it produces as output, so it doesn't have to be followed by a cast.

Cheers,
 - David

On Sat, Sep 22, 2018 at 10:29 PM Berti Krüger <[hidden email]> wrote:

Hi all.


I am very sorry to revive this ancient thread. 


But i thought that maybe it could help other vtk users in the future to have a compilable and working example of the final result. If someone thinks that this code could be useful as an vtk example then feel free to modify / correct / pimp / add this to the vtk examples on github.


I adapted the 3D slicer vtk code Csaba has pointed me to and incorporated the hints from David and everything works like a charm:



The whole procedure is incredibly simple. Basically the vtkPolyData of the mesh is converted via the vtkPolyDataToImageStencil filter to a vtkImageStencil and the rest of the whole magic is only done by the vtkImageStencil filter:


// Convert stencil to image
vtkNew<vtkImageStencil> stencil;
stencil->SetInputData(binaryLabelMap);
stencil->SetStencilConnection(polyDataToImageStencil->GetOutputPort());
stencil->ReverseStencilOn();
stencil->SetBackgroundValue(1); // General foreground value is 1 (background value because of reverse stencil)

I would really like to know very roughly how the vtkImageStencil works internally without digging through all of the source code. Because it is really fast (some kind of rasterization maybe ?). If someone could drop me a line on how it works only in principle  that would be really great.


So thanks again to Bill, Csaba, David, Elvis and everyone else on this thread.



Kind regards,


Berti



PS: Example code, CMakeLists and example STL mesh are attached


Von: Csaba Pinter <[hidden email]>
Gesendet: Montag, 11. Juni 2018 12:45
An: Berti Krüger
Cc: [hidden email]
Betreff: RE: [vtkusers] Solid Voxelization with VTK
 

Hi Berti,

 

That code is under BSD-style license, so you’re free to do basically anything with it. This implementation is already pretty simple; just a concatenation of several VTK filters, so you can just give it a try by copy-pasting the relevant part of that function into your code and see if it works. If you’re interested in using the vtkSegmentation library, let me know though!

 

Best,

csaba

 

From: Berti Krüger <[hidden email]>
Sent: Monday, June 11, 2018 02:46
To: Csaba Pinter <[hidden email]>; David Gobbi <[hidden email]>
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Csaba,

 

thank you very much for you answer and for offering your help.

 

I already have looked at it and i think if no simpler option will come up, i will give it a try (hoping not to have to change too much to get it working with my pipeline since i am an absolute vtk beginner :-) ).

 

How is this code licensed ? (i am only doing research work at the moment, the whole project i am working on is internal and just for testing and part of some larger work, it might be open sourced one day in the very far future but could even be not at all  and it will never be used commercially).

 

 

Regards

 

Berti

 

 

Am 05.06.2018 um 16:14 schrieb Csaba Pinter:

Hi Berti,

 

There is a pretty robust conversion algorithm implemented within 3D Slicer but using solely VTK:

https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

It consists of a short pipeline of standard VTK filters (some preprocessing then vtkPolyDataToImageStencil and vtkImageStencil).

 

Let me know if you have any questions.

 

csaba

 

 

From: vtkusers [hidden email] On Behalf Of Berti Krüger
Sent: Tuesday, June 5, 2018 03:26
To: David Gobbi [hidden email]
Cc: [hidden email]
Subject: Re: [vtkusers] Solid Voxelization with VTK

 

Hi David,

 

thank you very much for your help.

 

I tried the example from your link. Unfortunately the result is same as using the vtkVoxelModeller:

 


I get only a voxelization of the surface. The interior of the voxelized cone is still hollow and not filled with voxels (at least this is implied by the visualization).

 

For the visualization i use Bill Lorensen's example code from here:

https://lorensen.github.io/VTKExamples/site/Cxx/Medical/GenerateCubesFromLabels/

 

 

I have attached this small example which compiles with VTK 8.1 to show the issue. It basicly only consist of the example code from your link (https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/) and the visualization code from the  Bill Lorensen's example from the link above. I have not changed much.

 

Any idea what is going wrong or what can i do get also voxelization of the interior (solid voxelization)?

 

Thank you very much in advance.

 

 

Regards

 

Berti

 

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

Usage: GenerateCubesFromLabels FullHead.mhd StartLabel EndLabel where InputVolume is a meta file containing a 3 volume of discrete labels.

 

 


Von: David Gobbi <[hidden email]>
Gesendet: Montag, 28. Mai 2018 23:13
An: Berti Krüger
Cc: [hidden email]
Betreff: Re: [vtkusers] Solid Voxelization with VTK

 

Hi Berti,

 

If its a triangulated surface that you want to fill with voxels, use vtkPolyDataToImageStencil:

https://lorensen.github.io/VTKExamples/site/Cxx/PolyData/PolyDataToImageData/

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

 

 

If its a mesh of 3D elements and you want to sample the elements to create voxels, try vtkResampleToImage.

 

 - David

 

 

 

On Mon, May 28, 2018 at 4:41 PM, Berti Krüger <[hidden email]> wrote:

Hello Everyone.

For my Project i have to voxelize a 3D-Triangle-Mesh. I already found the
vtkVoxelModeller which, while somewhat slow (only around 200 Triangles per
second), works, but only voxelizes the outer shell where the mesh boundary
triangles are. The inner part of the mesh stays hollow.

Is there some way to get solid voxelization of 3D-Triangle-Meshes out of the
box with VTK ?


Thank you very much in advance,

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

PolyDataToImageDataVisualizationSTL.cxx (10K) Download Attachment
CMakeLists.txt (482 bytes) Download Attachment
eiffel_tower_small.stl (1014K) Download Attachment
12