Determining the voxels between two boundary surfaces

Question

Issue description

I am working on human brain tACS simulations where I have the models of the skin, skull, csf, brain and ventricles in STL format. The shape does not matter and there are no intersections. I want to run FE analysis and use SfePy but the only problem I have so far is how to define the volume regions. Using pymesh I am able to generate a pretty detailed and good mesh using the following code:

import pymesh

skin_stl = pymesh.load_mesh('skin.stl')
skull_stl = pymesh.load_mesh('skull.stl')
csf_stl = pymesh.load_mesh('csf.stl')
brain_stl = pymesh.load_mesh('brain.stl')

model = pymesh.merge_meshes( (skin_stl, skull_stl, csf_stl, brain_stl) )
model_tet = pymesh.tetrahedralize(model, 10)

The meshing I get from the code above can be seen in the following picture.

As you can see the boundaries are clear and the meshing algorithm took that into consideration. If I get the number of components of this mesh it is 4 as it should be since I have 4 bounding surfaces. Using the following code I am able to separate and get the vertices and faces of the above tessellated mesh.

new_mesh = pymesh.form_mesh(model_tet.vertices, model_tet.faces)
sep_mesh = pymesh.separate_mesh(new_mesh)

I can not in any way separate the mesh and get the voxels included in between two consecutive bounding surfaces. Is there a way to do that using pymesh or something else?

Answer 1

I finally did it! You can check the code on my GitLab repo. The function is to do this is the domain_extract which identifies the meshed boundaries to separate. To check a working example, you can check the sphere_meshing.py file, under the Scripts/Meshing directory.