Create 3D volume mesh

I have a 3D numpy array containing loads of processed binary images. I’d like to turn this into a 3D volume mesh. I’d rather keep this to python rather than using external software.

I have used this so far:

verts, faces, normals, values = marching_cubes_lewiner(all_trabecular_slices, 0.0) 
mlab.triangular_mesh([vert[0] for vert in verts],
                 [vert[1] for vert in verts],
                 [vert[2] for vert in verts],
                 normals)
mlab.savefig("Output\\Trabecular_Mask\\trabecular_model.obj")

but this only gives me the mesh of the “perimeter” or surface mesh of the shape:

3D Mesh

From a stack of images like this:

Sample of single image

How can I get this to be a “solid” 3D mesh?

1 Like

Hey @steersteer,

could you try to append a black image at the beginning and the end of your stack? The algorithm may then be able to see that your volume is closed…

Cheers,
Robert

2 Likes

Yes! That was it.
I’ve spent about 20 hours trying to work this out, can’t believe it was so simple.
Thanks

1 Like

You’re welcome. I think I spent 20 hours on thinking about that 3 years ago :wink:

2 Likes

Some other tips:

  1. you can use np.pad to make sure pad the cube’s every face.
  2. if you use the marching_cubes_lewiner’s step_size parameter, you may need a thick zero pad.
  3. you can do a gaussian filter to make the result surface smooth (but you may need a thick pad).
  4. you can use mesh_surface_area to count the surface area, (sum of all triangle’s area)
  5. if your surface is a closed volume, you can count the object’s volume. (sum of mixed moment of every triangle to the (0,0,0))
2 Likes

Do you have any suggestions for a package that could create a tetrahedral mesh? The triangle mesh only creates a “shell” but I’d like elements throughout. Something like PyMesh but this doesn’t work on windows yet

what 's your goal? you want to count the volume?

I do not think so. a isolate triangle is a surface, but all triangles build a volume, and it is more friendly to the gl engine. you can count surface area and volume. (if you need, I can tell you how to count). I think you can get all the geometry features by the triangles. except you want do a finite element analysis

and if you need tetrahedral indeed, you can follow the steps below:

  1. use all the triangle vertex to build a Delaunay net. (scipy.spatial.Delaunay) then every cell is a tetrahedral.
  2. but some tetrahedral is in the volume, others is out. you can choose the four vertex’s mean as a inner point, and use the original array mask to check if the mean point is true. you can remove all the out tetrahedral
1 Like

Yes, finite element analysis is my final goal for this. So I’ll need to use scipy.spatial.Delaunay for this?

I give you a demo in 2D.

from skimage.io import imread, imsave
from skimage.measure import find_contours, approximate_polygon
from scipy.spatial import Delaunay
import numpy as np
import matplotlib.pyplot as plt

img = imread('img.png')
ls = find_contours(img, 0)
ls = [approximate_polygon(i, 1) for i in ls]
vert = np.vstack(ls)
tri = Delaunay(vert).simplices

ax = plt.subplot(411)
ax.imshow(img, cmap='gray')

ax = plt.subplot(412)
ax.imshow(img, cmap='gray')
for i in ls: ax.plot(*i.T[::-1], 'red')

ax = plt.subplot(413)
ax.imshow(img, cmap='gray')
plt.triplot(*vert.T[::-1], tri)

ax = plt.subplot(414)
center = vert[tri[:,0]] + vert[tri[:,1]] + vert[tri[:,2]]
msk = img[tuple((center/3).astype(int).T)] == 0
ax.imshow(img, cmap='gray')
plt.triplot(*vert.T[::-1], tri[msk])
plt.show()


Is this what you want?
(there is a little mistake in e, that because it is too sparse when approximate, I think it is not exists in your case)

@sofroniewn , @jni Does napari support text? As I know treating text in 3d engine is not easy, (using texture).

An inspiration from this question. This could be a general method to show text in 3d engine.

  1. using pil draw in 2d array
  2. find_contours, and approximate. (maybe we need limit the approximate method must leave a vertex in specific distance, nomatter how straight the line is)
  3. Delaunay, and check every triangle in or out.
  4. send the inner triangles to the engine.

Hi @yxdragon, that approach sounds interesting.

Adding text to napari is in progress here, including 3D rendering and editing, resizing etc. We’re just using the default text rendering provided by the vispy library napari that uses under the hood, and I’m not sure about the details of their implementation, but you can find some of it here.

1 Like