Intersection between model and volume to select voxels

I have a generated a surface model by applying vtkMarchingCubes on a loaded volume. The next step I would like to perform is to identify which voxels intersect with the triangles that were generated. I was thinking of starting with iterating the vertices to see which voxels overlap and do some triangle intersect from there but I am not sure if that is the best method. Doe anybody know if there is a better way to select voxels that are inside, outside or on the contour of a model? I am using python.

What are you trying to achieve in general? I have the feeling that there should be an easier way for that in Slicer compared to low-level VTK programming.

Yes, you might be right. I would like to dilate the outer layer of the scan to replace this gradient-artifact in the scan with the color of the inside voxels. So only the replacing of this thin outer layer and leave the inside voxels untouched.

I found this Module:

So my guess is that in slicer there must be a way to identify which voxels are inside a model?

There is a complete infrastructure for this in the past few years called Segmentations. If you convert your model to a segmentation then set the reference volume to your volume in Segmentations module, then you can do all kinds of computations and editing.

How I’d do it:

  • Use latest 4.9.0 nightly
  • Create new segmentation node in Segmentations module
  • In Representations section make Closed surface as master (because you want to have control over the geometry of the labelmap afterwards)
  • Import model in Import/Export section
  • In Representations section long-click Create and choose Advanced create
  • Click the only path in the top part
  • Click Specify geometry button in the bottom part
  • Select your volume

Now you have a segmentation node which contains the labelmap representation of your model in the exact grid as your volume. Now you can do all kinds of things with it in the Segment Editor module (click yes when it asks if you want to change master to labelmap), such as dilation in the Margin effect. If you only need the dilated part then you can use Logical operators to copy the original segment, dilate the copy, and subtract the original from the dilated. You can calculate statistics in the Segment Statistics module.

Hi cpinter, thanks for the quick reply. I am not familiar with segmentations, I will try the steps you specify. If I understand correctly a segment model is created from the model that can be assigned to the volume to perform segment operations? In the end I would like to automate this process in a script and use the marchingcubes isovalue as a guide for automatically regenerating this segmentation.

You can do the whole segmentation in Segment Editor, and then you don’t have to even start from a model. You can do all that from python.

Here’s some information about Segmentations:
And about driving Segment Editor from python: (this is just a discourse search, please use this resource, a lot of questions have been answered already)
I think for your use case you can do the segmentation as a simple threshold, it gives you the same geometry as an isosurface (both are directly based on voxel intensity).