This page was generated from: notebooks/examples/how_to_generate_a_tetrahedral_mesh_from_3d_cells_image.ipynb

---

%load_ext autoreload
%autoreload 2
%config InlineBackend.rc = {'figure.figsize': (10,6)}
%matplotlib inline

Generate a tetrahedral mesh from 3D cells data

This notebook shows how to mesh a 3D volume:

1. Load and visualize a volume

2. Apply image filters and segment image

3. Generate a 3D surface mesh

4. Visualize and export the mesh to other formats

import pyvista as pv
from skimage import filters, data
import numpy as np

Load and vizualize the data

This example uses the cells3d sample data from scikit-image.

If you want to use your own data, any numpy array can be passed to a `Image <https://nanomesh.readthedocs.io/en/latest/nanomesh.volume.html#nanomesh.volume.Volume>`__ object. Data stored as .npy can be loaded using Image.load().

from nanomesh import Image
import numpy as np

from skimage.data import cells3d

data = cells3d()[:, 1, :, :]

vol = Image(data)

# swap axes so depth matches z
vol = vol.apply(np.swapaxes, axis1=0, axis2=2)

vol.show_slice()

<nanomesh.image._utils.SliceViewer at 0x213b1ccb3a0>

For this example, select a subvolume using .select_subvolume to keep the cpu times in check.

from skimage.transform import rescale

subvol = vol.select_subvolume(
    ys=(0, 128),
    zs=(128, -1),
)
subvol.show_slice()

<nanomesh.image._utils.SliceViewer at 0x2139bb41a00>

Nanomesh makes use of `itkwidgets <https://github.com/InsightSoftwareConsortium/itkwidgets>`__ to render the volumes.

subvol.show()

Filter and segment the data

Image segmentation is a way to label the pixels of different regions of interest in an image. In this example, we are interested in separating the bulk material (Si) from the nanopores. In the image, the Si is bright, and the pores are dark.

First, we apply a `gaussian filter <https://scikit-image.org/docs/dev/api/skimage.filters.html#skimage.filters.gaussian>`__ to smooth out some of the image noise to get a cleaner segmentation.

subvol_gauss = subvol.gaussian(sigma=1)
subvol_gauss.show_slice()

<nanomesh.image._utils.SliceViewer at 0x213b8ed0d60>

scikit-image contains a useful function to try all threshold finders on the data. These methods analyse the contrast histogram and try to find the optimal value to separate which parts of the image belong to each domain.

Since the function only works on a single slice, we first select a slice using the .select_plane method.

from skimage import filters

plane = subvol_gauss.select_plane(x=30)
plane.try_all_threshold(figsize=(5, 10))

We will use the yen method, because it gives nice separation.

The threshold value is used to segment the image using `np.digitize <https://numpy.org/doc/stable/reference/generated/numpy.digitize.html#numpy-digitize>`__.

subvol_gauss.binary_digitize(threshold='yen')

Volume(shape=(127, 128, 60), range=(0,1), dtype=int64)

subvol_seg = subvol_gauss.binary_digitize(threshold='yen')
subvol_seg.show_slice()

<nanomesh.image._utils.SliceViewer at 0x213b9604d00>

from scipy import ndimage


def fill_holes(image):
    return ndimage.binary_fill_holes(image).astype(int)


subvol_seg = subvol_seg.apply(fill_holes)
subvol_seg.show_slice()

<nanomesh.image._utils.SliceViewer at 0x213b95fb520>

Generate 3d tetragonal mesh

Meshes can be generated using the Mesher class. Meshing consists of two steps:

1. Contour finding (using the `marching_cubes <https://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.marching_cubes>`__ function

2. Triangulation (using the `tetgen <https://tetgen.pyvista.org/>`__ library)

Mesher requires a segmented image. generate_contour() wraps all domains of the image corresponding to that label. Here, 1 corresponds to the cells.

Meshing options are defined in the tetgen documentation. These can be specified using the opts parameter. The default options are opts='-pAq1.2:

• -A: Assigns attributes to tetrahedra in different regions.

• -p: Tetrahedralizes a piecewise linear complex (PLC).

• -q: Refines mesh (to improve mesh quality).

Also useful:

• -a: Applies a maximum tetrahedron volume constraint. Don’t make -a too small, or the algorithm will take a very long time to complete. If this parameter is left out, the triangles will keep growing without limit.

%%time

from nanomesh import Mesher

mesher = Mesher(subvol_seg)
mesher.generate_contour()
mesh = mesher.tetrahedralize(opts='-pAq')

Wall time: 16 s

Tetrahedralization returns a TetraMesh dataclass that can be used for various operations, for example showing the result using itkwidgets:

mesh.plot_pyvista(jupyter_backend='static', show_edges=True)

Using region markers

By default, the region attributes are assigned automatically by tetgen. Tetrahedra in each enclosed region will be assigned a new label sequentially.

Region markers are used to assign attributes to tetrahedra in different regions. After tetrahedralization, the region markers will ‘flood’ the regions up to the defined boundaries. The elements of the resulting mesh are marked according to the region they belong to (tetras.metadata['tetgenRef'].

You can view the existing region markers by looking at the .region_markers attribute on the contour.

mesher.contour.region_markers

RegionMarkerList(
    RegionMarker(label=1, point=(44, 76, 48), name='background', constraint=0),
    RegionMarker(label=2, point=(47, 46, 33), name='feature', constraint=0),
    RegionMarker(label=2, point=(2, 77, 38), name='feature', constraint=0),
    RegionMarker(label=2, point=(0, 90, 50), name='feature', constraint=0),
    RegionMarker(label=2, point=(31, 98, 34), name='feature', constraint=0),
    RegionMarker(label=2, point=(9, 59, 46), name='feature', constraint=0),
    RegionMarker(label=2, point=(15, 76, 48), name='feature', constraint=0),
    RegionMarker(label=2, point=(23, 69, 51), name='feature', constraint=0),
    RegionMarker(label=2, point=(69, 123, 35), name='feature', constraint=0),
    RegionMarker(label=2, point=(43, 31, 59), name='feature', constraint=0),
    RegionMarker(label=2, point=(48, 50, 59), name='feature', constraint=0),
    RegionMarker(label=2, point=(51.15384615384615, 1.6153846153846154, 50.23076923076923), name='feature', constraint=0),
    RegionMarker(label=2, point=(87, 6, 38), name='feature', constraint=0),
    RegionMarker(label=2, point=(68, 77, 38), name='feature', constraint=0),
    RegionMarker(label=2, point=(73, 92, 47), name='feature', constraint=0),
    RegionMarker(label=1, point=(88, 2, 38), name='background', constraint=0),
    RegionMarker(label=2, point=(119, 107, 41), name='feature', constraint=0),
    RegionMarker(label=2, point=(99, 21, 47), name='feature', constraint=0),
    RegionMarker(label=2, point=(110, 100, 23), name='feature', constraint=0),
    RegionMarker(label=2, point=(113.78571428571429, 78.64285714285714, 47.285714285714285), name='feature', constraint=0),
    RegionMarker(label=2, point=(124, 2, 34), name='feature', constraint=0),
    RegionMarker(label=2, point=(125, 24, 35), name='feature', constraint=0),
    RegionMarker(label=2, point=(126, 76, 45), name='feature', constraint=0),
    RegionMarker(label=2, point=(126, 74, 49), name='feature', constraint=0)
)

Mesh evaluation

The mesh can be evaluated using the metrics module. This example shows how to calculate all metrics and plot them on a section through the generated mesh.

from nanomesh import metrics

tetra_mesh = mesh.get('tetra')

metrics_dict = metrics.calculate_all_metrics(tetra_mesh, inplace=True)
metrics_dict

{'area': array([-1., -1., -1., ..., -1., -1., -1.]),
 'aspect_frobenius': array([2.31677723, 1.19575747, 1.44029751, ..., 1.75947687, 1.4392653 ,
        1.24123145]),
 'aspect_ratio': array([3.80024454, 1.42080963, 1.84731578, ..., 2.71672234, 1.87260525,
        1.49231281]),
 'condition': array([3.05341788, 1.21351406, 1.50935969, ..., 1.97026528, 1.38174864,
        1.28635338]),
 'max_angle': array([-1., -1., -1., ..., -1., -1., -1.]),
 'min_angle': array([14.60862095, 43.95133677, 37.62991895, ..., 29.54078766,
        49.42751675, 44.49077778]),
 'radius_ratio': array([3.25667802, 1.23811768, 1.58989383, ..., 2.63296245, 1.62480334,
        1.33489297]),
 'scaled_jacobian': array([0.25730214, 0.50833528, 0.44610067, ..., 0.39530155, 0.33271911,
        0.60030157]),
 'shape': array([0.43163408, 0.83628999, 0.694301  , ..., 0.56835075, 0.69479894,
        0.80565152]),
 'relative_size_squared': array([0., 0., 0., ..., 0., 0., 0.]),
 'shape_and_size': array([0., 0., 0., ..., 0., 0., 0.]),
 'distortion': array([1., 1., 1., ..., 1., 1., 1.]),
 'max_min_edge_ratio': array([1.41421356, 1.99347457, 1.92068054, ..., 1.58555878, 1.76211173,
        1.54166293])}

Using the .plot_submesh() method, any array that is present in the metadata can be plotted. plot_submesh() is flexible, in that it can show a slice through the mesh as defined using index, along, and invert. Extra keyword arguments, such as show_edges and lighting are passed on to `Plotter.add_mesh() <https://docs.pyvista.org/api/plotting/_autosummary/pyvista.Plotter.add_mesh.html?highlight=add_mesh>`__.

tetra_mesh.plot_submesh(
    along='x',
    index=30,
    scalars='min_angle',
    show_edges=True,
    lighting=True,
    backend='static',
)

<pyvista.plotting.plotting.Plotter at 0x213b30d26d0>

Interoperability

The TetraMesh object can also be used to convert to various other library formats, such as:

• `trimesh.open3d <http://www.open3d.org/docs/release/python_api/open3d.geometry.TetraMesh.html#open3d.geometry.TetraMesh>`__

• `pyvista.UnstructuredGrid <https://docs.pyvista.org/examples/00-load/create-unstructured-surface.html>`__

• `meshio.Mesh <https://docs.pyvista.org/examples/00-load/create-unstructured-surface.html>`__

To save the data, use the .write method. This is essentially a very thin wrapper around meshio, equivalent to meshio_mesh.write(...).

tetra_mesh.write('cells.vtk')

---

Generated by nbsphinx from a Jupyter notebook.