Example with dolfinx#
import subprocess
from pathlib import Path
from mpi4py import MPI
import pyvista
import gmsh
import dolfinx
Let us first list the help menu
subprocess.run(["ukb-atlas", "--help"])
usage: ukb-atlas [-h] {surf,clip,mesh} ...
UKB-atlas This is a command line interface for extracting surfaces and
generating Bi-ventricular meshes from the UK Biobank atlas:
https://www.cardiacatlas.org/biventricular-modes/
positional arguments:
{surf,clip,mesh}
surf Extract surfaces from the atlas
clip Clip the surfaces
mesh Generate mesh from the surfaces
options:
-h, --help show this help message and exit
CompletedProcess(args=['ukb-atlas', '--help'], returncode=0)
We will first use the surf command
subprocess.run(["ukb-atlas", "surf", "--help"])
usage: ukb-atlas surf [-h] [-a] [-m MODE] [-s STD] [-v]
[--cache-dir CACHE_DIR] [-c {ED,ES,both}] [--use-burns]
[--burns-path BURNS_PATH]
folder
positional arguments:
folder Directory to save the generated surfaces.
options:
-h, --help show this help message and exit
-a, --all Use the the PCA atlas derived from all 4,329 subjects
from the UK Biobank Study. By default we use the PCA
atlas derived from 630 healthy reference subjects from
the UK Biobank Study (default: False)
-m MODE, --mode MODE Mode to generate points from. If -1, generate points
from the mean shape. If between 0 and the number of
modes, generate points from the specified mode. By
default -1 (default: -1)
-s STD, --std STD Standard deviation to scale the mode by, by default
1.5 (default: 1.5)
-v, --verbose Print verbose output. (default: False)
--cache-dir CACHE_DIR
Directory to save the downloaded atlas. Can also be
set with the UKB_CACHE_DIR environment variable. By
default ~/.ukb (default: /github/home/.ukb)
-c {ED,ES,both}, --case {ED,ES,both}
Case to generate surfaces for. (default: ED)
--use-burns Use the atlas from Richard Burns to generate the
surfaces. (default: False)
--burns-path BURNS_PATH
Path to the burns atlas file. This will be a .mat file
which will be loaded using scipy.io.loadmat. This
needs to be specified if --use-burns is set. (default:
None)
CompletedProcess(args=['ukb-atlas', 'surf', '--help'], returncode=0)
First we try to generate the surfaces for mean shape from the UK Biobank atlas end-diastole
folder = Path("data-mean")
mode = "-1" # -1 refers to mean shape (otherwise you can choice a value between 0 and 200)
std = "1.5" # Standard deviation to scale the mode by, by default 1.5
subprocess.run(["ukb-atlas", "surf", str(folder), "--mode", mode, "--std", std, "--case", "ED"])
INFO:ukb.atlas:Downloading https://www.cardiacatlas.org/share/download.php?id=60&token=AR3JSoaxJ9Ev9n8QAkvV4BHJUniyttqm&download to /github/home/.ukb/UKBRVLV.zip. This may take a while.
INFO:ukb.atlas:Done downloading.
INFO:ukb.atlas:Generating points from /github/home/.ukb/UKBRVLV.h5 using mode -1 and std 1.5
INFO:ukb.surface:Saved data-mean/EPI_ED.stl
INFO:ukb.surface:Saved data-mean/MV_ED.stl
INFO:ukb.surface:Saved data-mean/AV_ED.stl
INFO:ukb.surface:Saved data-mean/TV_ED.stl
INFO:ukb.surface:Saved data-mean/PV_ED.stl
INFO:ukb.surface:Saved data-mean/LV_ED.stl
INFO:ukb.surface:Saved data-mean/RV_ED.stl
INFO:ukb.surface:Saved data-mean/RVFW_ED.stl
CompletedProcess(args=['ukb-atlas', 'surf', 'data-mean', '--mode', '-1', '--std', '1.5', '--case', 'ED'], returncode=0)
Let us now see which files that was generated
subprocess.run(["ls", str(folder)])
AV_ED.stl
EPI_ED.stl
LV_ED.stl
MV_ED.stl
PV_ED.stl
RVFW_ED.stl
RV_ED.stl
TV_ED.stl
parameters.json
CompletedProcess(args=['ls', 'data-mean'], returncode=0)
We will now generate the mesh for the mean shape, but let us first look at the help menu
subprocess.run(["ukb-atlas", "mesh", "--help"])
usage: ukb-atlas mesh [-h] [--char_length_max CHAR_LENGTH_MAX]
[--char_length_min CHAR_LENGTH_MIN] [-v]
[-c {ED,ES,both}] [--clipped]
folder
positional arguments:
folder Directory to save the generated meshes.
options:
-h, --help show this help message and exit
--char_length_max CHAR_LENGTH_MAX
Maximum characteristic length of the mesh elements.
(default: 5.0)
--char_length_min CHAR_LENGTH_MIN
Minimum characteristic length of the mesh elements.
(default: 5.0)
-v, --verbose Print verbose output. (default: False)
-c {ED,ES,both}, --case {ED,ES,both}
Case to generate surfaces for. (default: ED)
--clipped Create a clipped mesh. (default: False)
CompletedProcess(args=['ukb-atlas', 'mesh', '--help'], returncode=0)
Note that this command uses gmsh, which is not part of the default dependencies so this needs to be installed. You can do this using either
python3 -m pip install gmsh
or use
python3 -m pip install ukb-atlas[gmsh]
Now we will generate the mesh for the mean shape (let us pick the rest of the parameters as default)
subprocess.run(["ukb-atlas", "mesh", str(folder)])
INFO:ukb.mesh:Creating mesh for ED with char_length_max=5.0, char_length_min=5.0
INFO:ukb.mesh:Created mesh data-mean/ED.msh
CompletedProcess(args=['ukb-atlas', 'mesh', 'data-mean'], returncode=0)
Note that by decreasing the values of char_length_maxand char_length_min you will get a finer mesh. We can now load the mesh with dolfinx and plot it with pyvista
comm = MPI.COMM_WORLD
msh_file = folder / "ED.msh"
gmsh.initialize()
gmsh.model.add("Mesh from file")
gmsh.merge(str(msh_file))
mesh, ct, ft = dolfinx.io.gmshio.model_to_mesh(gmsh.model, comm, 0)
markers = {
gmsh.model.getPhysicalName(*v): tuple(reversed(v)) for v in gmsh.model.getPhysicalGroups()
}
print(markers)
gmsh.finalize()
Info : Reading 'data-mean/ED.msh'...
Info : 27 entities
Info : 3788 nodes
Info : 18941 elements
Info : 8 parametrizations
Info : [ 0%] Processing parametrizations
Info : [ 10%] Processing parametrizations
Info : [ 20%] Processing parametrizations
Info : [ 30%] Processing parametrizations
Info : [ 40%] Processing parametrizations
Info : [ 60%] Processing parametrizations
Info : [ 70%] Processing parametrizations
Info : [ 80%] Processing parametrizations
Info : Done reading 'data-mean/ED.msh'
{'LV': (1, 2), 'RV': (2, 2), 'MV': (3, 2), 'AV': (4, 2), 'PV': (5, 2), 'TV': (6, 2), 'EPI': (7, 2), 'Wall': (9, 3)}
pyvista.start_xvfb()
vtk_mesh = dolfinx.plot.vtk_mesh(mesh, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
figure = plotter.screenshot("mesh.png")
vtk_bmesh = dolfinx.plot.vtk_mesh(mesh, ft.dim, ft.indices)
bgrid = pyvista.UnstructuredGrid(*vtk_bmesh)
bgrid.cell_data["Facet tags"] = ft.values
bgrid.set_active_scalars("Facet tags")
p = pyvista.Plotter(window_size=[800, 800])
p.add_mesh(bgrid, show_edges=True)
if not pyvista.OFF_SCREEN:
p.show()
else:
figure = p.screenshot("facet_tags_biv.png")
Now lets us perturb the mean shape with the second mode
folder_mode_1 = Path("data-mode-1")
subprocess.run(
["ukb-atlas", "surf", str(folder_mode_1), "--mode", "1", "--std", "1.5", "--case", "ED"]
)
INFO:ukb.atlas:Generating points from /github/home/.ukb/UKBRVLV.h5 using mode 1 and std 1.5
INFO:ukb.surface:Saved data-mode-1/EPI_ED.stl
INFO:ukb.surface:Saved data-mode-1/MV_ED.stl
INFO:ukb.surface:Saved data-mode-1/AV_ED.stl
INFO:ukb.surface:Saved data-mode-1/TV_ED.stl
INFO:ukb.surface:Saved data-mode-1/PV_ED.stl
INFO:ukb.surface:Saved data-mode-1/LV_ED.stl
INFO:ukb.surface:Saved data-mode-1/RV_ED.stl
INFO:ukb.surface:Saved data-mode-1/RVFW_ED.stl
CompletedProcess(args=['ukb-atlas', 'surf', 'data-mode-1', '--mode', '1', '--std', '1.5', '--case', 'ED'], returncode=0)
And create the mesh
subprocess.run(["ukb-atlas", "mesh", str(folder_mode_1)])
INFO:ukb.mesh:Creating mesh for ED with char_length_max=5.0, char_length_min=5.0
INFO:ukb.mesh:Created mesh data-mode-1/ED.msh
CompletedProcess(args=['ukb-atlas', 'mesh', 'data-mode-1'], returncode=0)
comm = MPI.COMM_WORLD
msh_file = folder_mode_1 / "ED.msh"
gmsh.initialize()
gmsh.model.add("Mesh from file")
gmsh.merge(str(msh_file))
mesh2, ct2, ft2 = dolfinx.io.gmshio.model_to_mesh(gmsh.model, comm, 0)
markers = {
gmsh.model.getPhysicalName(*v): tuple(reversed(v)) for v in gmsh.model.getPhysicalGroups()
}
print(markers)
gmsh.finalize()
Info : Reading 'data-mode-1/ED.msh'...
Info : 27 entities
Info : 3959 nodes
Info : 19809 elements
Info : 8 parametrizations
Info : [ 0%] Processing parametrizations
Info : [ 10%] Processing parametrizations
Info : [ 20%] Processing parametrizations
Info : [ 30%] Processing parametrizations
Info : [ 40%] Processing parametrizations
Info : [ 60%] Processing parametrizations
Info : [ 70%] Processing parametrizations
Info : [ 80%] Processing parametrizations
Info : Done reading 'data-mode-1/ED.msh'
{'LV': (1, 2), 'RV': (2, 2), 'MV': (3, 2), 'AV': (4, 2), 'PV': (5, 2), 'TV': (6, 2), 'EPI': (7, 2), 'Wall': (9, 3)}
pyvista.start_xvfb()
vtk_mesh2 = dolfinx.plot.vtk_mesh(mesh2, mesh2.topology.dim)
grid2 = pyvista.UnstructuredGrid(*vtk_mesh2)
plotter2 = pyvista.Plotter()
plotter2.add_mesh(grid2, show_edges=False)
plotter2.add_mesh(grid, style="wireframe", color="k")
plotter2.view_xy()
if not pyvista.OFF_SCREEN:
plotter2.show()
else:
figure = plotter.screenshot("mesh.png")
vtk_bmesh2 = dolfinx.plot.vtk_mesh(mesh2, ft2.dim, ft2.indices)
bgrid2 = pyvista.UnstructuredGrid(*vtk_bmesh2)
bgrid2.cell_data["Facet tags"] = ft2.values
bgrid2.set_active_scalars("Facet tags")
p2 = pyvista.Plotter(window_size=[800, 800])
p2.add_mesh(bgrid2, show_edges=True)
if not pyvista.OFF_SCREEN:
p2.show()
else:
figure = p2.screenshot("facet_tags_biv.png")
Clipping the mesh#
Let us now show how to clip the mesh. First let us just clip the first mesh using pyvista to see how we can the mesh to look like. This is also a good way to figure out what should be the origin and normal for the clip plane. (Note that the default values are the ones given below which works well for the mean shape of the ED mesh)
origin = (-13.612554383622273, 18.55767189380559, 15.135103714006394)
normal = (-0.7160843664428893, 0.544394641424108, 0.4368725838557541)
p3 = pyvista.Plotter()
p3.add_mesh_clip_plane(grid, normal=normal, origin=origin, invert=True)
if not pyvista.OFF_SCREEN:
p3.show()
else:
figure = p3.screenshot("mesh_clip.png")
There is a command called clip, so let us first look at the help section
subprocess.run(["ukb-atlas", "clip", "--help"])
usage: ukb-atlas clip [-h] [-c {ED,ES,both}] [-ox ORIGIN_X] [-oy ORIGIN_Y]
[-oz ORIGIN_Z] [-nx NORMAL_X] [-ny NORMAL_Y]
[-nz NORMAL_Z] [-s] [-si SMOOTH_ITER]
[-sr SMOOTH_RELAXATION] [-v]
folder
positional arguments:
folder Directory to save the generated surfaces.
options:
-h, --help show this help message and exit
-c {ED,ES,both}, --case {ED,ES,both}
Case to generate surfaces for. (default: ED)
-ox ORIGIN_X, --origin-x ORIGIN_X
Origin of the clipping plane in x direction. (default:
-13.612554383622273)
-oy ORIGIN_Y, --origin-y ORIGIN_Y
Origin of the clipping plane in y direction. (default:
18.55767189380559)
-oz ORIGIN_Z, --origin-z ORIGIN_Z
Origin of the clipping plane in z direction. (default:
15.135103714006394)
-nx NORMAL_X, --normal-x NORMAL_X
Normal of the clipping plane in x direction. (default:
-0.7160843664428893)
-ny NORMAL_Y, --normal-y NORMAL_Y
Normal of the clipping plane in y direction. (default:
0.544394641424108)
-nz NORMAL_Z, --normal-z NORMAL_Z
Normal of the clipping plane in z direction. (default:
0.4368725838557541)
-s, --smooth Smooth the RV surface. (default: False)
-si SMOOTH_ITER, --smooth-iter SMOOTH_ITER
Number of iterations to smooth the RV surface.
(default: 100)
-sr SMOOTH_RELAXATION, --smooth-relaxation SMOOTH_RELAXATION
Relaxation factor to smooth the RV surface. (default:
0.1)
-v, --verbose Print verbose output. (default: False)
CompletedProcess(args=['ukb-atlas', 'clip', '--help'], returncode=0)
Note that this command uses pyvista, which is not part of the default dependencies so this needs to be installed. You can do this using either
python3 -m pip install pyvista
or use
python3 -m pip install ukb-atlas[pyvista]
To create the mesh we pass the origin and normal, we will also smooth the RV endocardium because the original surface has some very sharp transistions
subprocess.run(
[
"ukb-atlas",
"clip",
str(folder),
"--smooth",
"--case",
"ED",
"-ox",
str(origin[0]),
"-oy",
str(origin[1]),
"-oz",
str(origin[2]),
"-nx",
str(normal[0]),
"-ny",
str(normal[1]),
"-nz",
str(normal[2]),
]
)
INFO:ukb.clip:Folder: data-mean
INFO:ukb.clip:Case: ED
INFO:ukb.clip:Origin: [-13.612554383622273, 18.55767189380559, 15.135103714006394]
INFO:ukb.clip:Normal: [-0.7160843664428893, 0.544394641424108, 0.4368725838557541]
INFO:ukb.clip:Reading data-mean/LV_ED.stl
Warning: PLY writer doesn't support multidimensional point data yet. Skipping
Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.clip:Saved data-mean/lv_clipped.ply
INFO:ukb.clip:Reading data-mean/RV_ED.stl
INFO:ukb.clip:Reading data-mean/RVFW_ED.stl
INFO:ukb.clip:Merging RV and RVFW
INFO:ukb.clip:Smoothing RV
INFO:ukb.clip:Saving data-mean/rv_clipped.ply
Warning: PLY writer doesn't support multidimensional point data yet. Skipping
Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.clip:Reading data-mean/EPI_ED.stl
INFO:ukb.clip:Saving data-mean/epi_clipped.ply
Warning: PLY writer doesn't support multidimensional point data yet. Skipping
Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
CompletedProcess(args=['ukb-atlas', 'clip', 'data-mean', '--smooth', '--case', 'ED', '-ox', '-13.612554383622273', '-oy', '18.55767189380559', '-oz', '15.135103714006394', '-nx', '-0.7160843664428893', '-ny', '0.544394641424108', '-nz', '0.4368725838557541'], returncode=0)
Let us again see which files we have
subprocess.run(["ls", str(folder)])
AV_ED.stl
ED.msh
EPI_ED.stl
LV_ED.stl
MV_ED.stl
PV_ED.stl
RVFW_ED.stl
RV_ED.stl
TV_ED.stl
epi_clipped.ply
lv_clipped.ply
parameters.json
rv_clipped.ply
CompletedProcess(args=['ls', 'data-mean'], returncode=0)
There are now some new .ply files that represents the clipped surfaces. We can for example take a look at one of them
lv_clipped = pyvista.read(folder / "lv_clipped.ply")
p4 = pyvista.Plotter()
p4.add_mesh(lv_clipped)
if not pyvista.OFF_SCREEN:
p4.show()
else:
figure = p4.screenshot("lv_clipped.png")
To mesh to three surfaces togeher we repeat the mesh command but pass the --clipped flag
subprocess.run(["ukb-atlas", "mesh", str(folder), "--clipped"])
INFO:ukb.mesh:Creating clipped mesh for ED with char_length_max=5.0, char_length_min=5.0
INFO:ukb.mesh:Created mesh data-mean/ED_clipped.msh
CompletedProcess(args=['ukb-atlas', 'mesh', 'data-mean', '--clipped'], returncode=0)
We can now load the new clipped mesh into dolfinx
comm = MPI.COMM_WORLD
msh_file = folder / "ED_clipped.msh"
gmsh.initialize()
gmsh.model.add("Mesh from file")
gmsh.merge(str(msh_file))
mesh_clipped, ct_clipped, ft_clipped = dolfinx.io.gmshio.model_to_mesh(gmsh.model, comm, 0)
markers_clipped = {
gmsh.model.getPhysicalName(*v): tuple(reversed(v)) for v in gmsh.model.getPhysicalGroups()
}
print(markers_clipped)
gmsh.finalize()
Info : Reading 'data-mean/ED_clipped.msh'...
Info : 11 entities
Info : 2603 nodes
Info : 13008 elements
Info : 3 parametrizations
Info : [ 0%] Processing parametrizations
Info : [ 10%] Processing parametrizations
Info : [ 40%] Processing parametrizations
Info : Done reading 'data-mean/ED_clipped.msh'
{'LV': (1, 2), 'RV': (2, 2), 'EPI': (3, 2), 'BASE': (4, 2), 'Wall': (5, 3)}
and visualize it
pyvista.start_xvfb()
vtk_mesh_clipped = dolfinx.plot.vtk_mesh(mesh_clipped, mesh_clipped.topology.dim)
grid_clipped = pyvista.UnstructuredGrid(*vtk_mesh_clipped)
p5 = pyvista.Plotter()
p5.add_mesh(grid_clipped, show_edges=True)
# p5.add_mesh(grid, style="wireframe", color="k")
p5.view_xy()
if not pyvista.OFF_SCREEN:
p5.show()
else:
figure = p5.screenshot("mesh_clipped.png")
vtk_bmesh_clipped = dolfinx.plot.vtk_mesh(mesh_clipped, ft_clipped.dim, ft_clipped.indices)
bgrid_clipped = pyvista.UnstructuredGrid(*vtk_bmesh_clipped)
bgrid_clipped.cell_data["Facet tags"] = ft_clipped.values
bgrid_clipped.set_active_scalars("Facet tags")
p6 = pyvista.Plotter(window_size=[800, 800])
p6.add_mesh(bgrid_clipped, show_edges=True)
if not pyvista.OFF_SCREEN:
p6.show()
else:
figure = p6.screenshot("facet_tags_biv_clipped.png")