UK-Biobank atlas#
The UK-Biobank provides an atlas for BiV meshes (see https://www.cardiacatlas.org/biventricular-modes/). We also provide a separate package to generate meshes from this atlas, see https://computationalphysiology.github.io/ukb-atlas.
However, if you want to use these meshes for cardiac simulations you probably need to generate fiber orientations using `fenicsx-ldrb as well. Here we will show how to do this programmatically from python
Select the mode you want to use (-1 will you the mean mode)
mode = -1
# Standard deviation of the mode
std = 1.5
Choose output directory
# Characteristic length of the mesh (smaller values will give finer meshes). Let us also pick the case from end-diastole (ED)
char_length = 10.0
case = "ED"
Generate surfaces
ukb.cli.main(["surf", str(outdir), "--mode", str(mode), "--std", str(std), "--case", case])
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 ukb_mesh/EPI_ED.stl
INFO:ukb.surface:Saved ukb_mesh/MV_ED.stl
INFO:ukb.surface:Saved ukb_mesh/AV_ED.stl
INFO:ukb.surface:Saved ukb_mesh/TV_ED.stl
INFO:ukb.surface:Saved ukb_mesh/PV_ED.stl
INFO:ukb.surface:Saved ukb_mesh/LV_ED.stl
INFO:ukb.surface:Saved ukb_mesh/RV_ED.stl
INFO:ukb.surface:Saved ukb_mesh/RVFW_ED.stl
0
Generate mesh
ukb.cli.main([
"mesh",
str(outdir),
"--case",
case,
"--char_length_max",
str(char_length),
"--char_length_min",
str(char_length),
])
INFO:ukb.mesh:Creating mesh for ED with char_length_max=10.0, char_length_min=10.0
INFO:ukb.mesh:Created mesh ukb_mesh/ED.msh
0
# Choose the ED mesh
mesh_name = outdir / "ED.msh"
# Convert mesh to dolfinx
comm = MPI.COMM_WORLD
geometry = cg.utils.gmsh2dolfin(comm=comm, msh_file=mesh_name)
2025-10-06 08:42:54 [debug ] Convert file ukb_mesh/ED.msh to dolfin
Info : Reading 'ukb_mesh/ED.msh'...
Info : 27 entities
Info : 2953 nodes
Info : 14635 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 'ukb_mesh/ED.msh'
Now we can plot the mesh
vtk_mesh = dolfinx.plot.vtk_mesh(geometry.mesh, geometry.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("ukb_mesh.png")
2025-10-06 08:42:55.348 ( 11.694s) [ 7F5262D87140]vtkXOpenGLRenderWindow.:1458 WARN| bad X server connection. DISPLAY=
INFO:trame_server.utils.namespace:Translator(prefix=None)
INFO:wslink.backends.aiohttp:awaiting runner setup
INFO:wslink.backends.aiohttp:awaiting site startup
INFO:wslink.backends.aiohttp:Print WSLINK_READY_MSG
INFO:wslink.backends.aiohttp:Schedule auto shutdown with timout 0
INFO:wslink.backends.aiohttp:awaiting running future
We als save the markers to json
print(geometry.markers)
if comm.rank == 0:
(outdir / "markers.json").write_text(
json.dumps(geometry.markers, default=cg.utils.json_serial)
)
comm.barrier()
{'LV': (1, 2), 'RV': (2, 2), 'MV': (3, 2), 'AV': (4, 2), 'PV': (5, 2), 'TV': (6, 2), 'EPI': (7, 2), 'Wall': (9, 3)}
And plot the facet tags
assert geometry.ffun is not None
vtk_bmesh = dolfinx.plot.vtk_mesh(
geometry.mesh, geometry.ffun.dim, geometry.ffun.indices
)
bgrid = pyvista.UnstructuredGrid(*vtk_bmesh)
bgrid.cell_data["Facet tags"] = geometry.ffun.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_ukb.png")
We could also combine all the outflow tracts into one entity, i.e the BASE
entities = [
geometry.ffun.find(1), # LV
geometry.ffun.find(2), # RV
geometry.ffun.find(7), # EPI
np.hstack(
[
geometry.ffun.find(3),
geometry.ffun.find(4),
geometry.ffun.find(5),
geometry.ffun.find(6),
]
), # BASE
]
entity_values = [
np.full(e.shape, i + 1, dtype=np.int32) for i, e in enumerate(entities)
]
and create new mesh tags
ffun = dolfinx.mesh.meshtags(
geometry.mesh, 2, np.hstack(entities), np.hstack(entity_values)
)
and plot the new tags
vtk_bmesh_new = dolfinx.plot.vtk_mesh(geometry.mesh, ffun.dim, ffun.indices)
bgrid_new = pyvista.UnstructuredGrid(*vtk_bmesh_new)
bgrid_new.cell_data["Facet tags"] = ffun.values
bgrid_new.set_active_scalars("Facet tags")
p = pyvista.Plotter(window_size=[800, 800])
p.add_mesh(bgrid_new, show_edges=True)
if not pyvista.OFF_SCREEN:
p.show()
else:
figure = p.screenshot("new_facet_tags_ukb.png")
Let us also save the markers so that we can inspect them later
with dolfinx.io.XDMFFile(geometry.mesh.comm, outdir / "new_ffun.xdmf", "w") as xdmf:
xdmf.write_mesh(geometry.mesh)
xdmf.write_meshtags(ffun, geometry.mesh.geometry)
Now let ut combine the markers as well so that they align with the markers needed in fenicsx-ldrb
markers = {}
keymap = {"LV": "lv", "RV": "rv", "EPI": "epi"}
for k, v in geometry.markers.items():
markers[keymap.get(k, k)] = [v[0]]
markers["base"] = [geometry.markers[m][0] for m in ["MV", "PV", "TV", "AV"]]
And finally generate the fibers
fiber_space = "P_1"
fiber_angle_endo = 60
fiber_angle_epi = -60
system = ldrb.dolfinx_ldrb(
mesh=geometry.mesh,
ffun=geometry.ffun,
markers=markers,
alpha_endo_lv=fiber_angle_endo,
alpha_epi_lv=fiber_angle_epi,
beta_endo_lv=0,
beta_epi_lv=0,
fiber_space=fiber_space,
)
INFO:ldrb.ldrb:Calculating scalar fields
WARNING:root:Unknown marker MV. Expected one of {'lv', 'rv', 'mv', 'av', 'epi', 'base'}
WARNING:root:Unknown marker AV. Expected one of {'lv', 'rv', 'mv', 'av', 'epi', 'base'}
WARNING:root:Unknown marker PV. Expected one of {'lv', 'rv', 'mv', 'av', 'epi', 'base'}
WARNING:root:Unknown marker TV. Expected one of {'lv', 'rv', 'mv', 'av', 'epi', 'base'}
WARNING:root:Unknown marker Wall. Expected one of {'lv', 'rv', 'mv', 'av', 'epi', 'base'}
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
rv: [2]
MV: [3]
AV: [4]
PV: [5]
TV: [6]
epi: [7]
Wall: [9]
base: [3, 5, 6, 4]
INFO:ldrb.ldrb: Num vertices: 2953
INFO:ldrb.ldrb: Num cells: 8977
INFO:ldrb.ldrb: Apex coord: (45.82, -21.12, -28.18)
INFO:ldrb.ldrb:
Calculating gradients
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 60
epi_lv: -60
endo_septum: 60
epi_septum: -60
endo_rv: 60
epi_rv: -60
INFO:ldrb.ldrb:beta:
endo_lv: 0
epi_lv: 0
endo_septum: 0
epi_septum: 0
endo_rv: 0
epi_rv: 0
and save them
cg.fibers.utils.save_microstructure(
mesh=geometry.mesh, functions=(system.f0, system.s0, system.n0), outdir=outdir
)
# Let us also plot the fibers
topology_f0, cell_types_f0, geometry_f0 = dolfinx.plot.vtk_mesh(
system.f0.function_space
)
values = np.zeros((geometry_f0.shape[0], 3), dtype=np.float64)
values[:, : len(system.f0)] = system.f0.x.array.real.reshape(
(geometry_f0.shape[0], len(system.f0))
)
function_grid = pyvista.UnstructuredGrid(topology_f0, cell_types_f0, geometry_f0)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=1.0)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe", color="r")
plotter.add_mesh(glyphs)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
fig_as_array = plotter.screenshot("fiber_ukb.png")
You should now be able to load the geometry and the fibers using the following code
geo = cg.geometry.Geometry.from_folder(comm=comm, folder=outdir)
INFO:cardiac_geometries.geometry:Reading geometry from ukb_mesh
We have also implemented a wrapper around the ukb-atlas package that allows you to generate the fibers directly. In this case you can simply de
INFO:ukb.atlas:Generating points from /github/home/.ukb/UKBRVLV.h5 using mode -1 and std 1.5
INFO:ukb.surface:Saved ukb_mesh2/EPI_ED.stl
INFO:ukb.surface:Saved ukb_mesh2/MV_ED.stl
INFO:ukb.surface:Saved ukb_mesh2/AV_ED.stl
INFO:ukb.surface:Saved ukb_mesh2/TV_ED.stl
INFO:ukb.surface:Saved ukb_mesh2/PV_ED.stl
INFO:ukb.surface:Saved ukb_mesh2/LV_ED.stl
INFO:ukb.surface:Saved ukb_mesh2/RV_ED.stl
INFO:ukb.surface:Saved ukb_mesh2/RVFW_ED.stl
INFO:ukb.mesh:Creating mesh for ED with char_length_max=2.0, char_length_min=2.0
0
INFO:ukb.mesh:Created mesh ukb_mesh2/ED.msh
2025-10-06 08:43:06 [debug ] Convert file ukb_mesh2/ED.msh to dolfin
Info : Reading 'ukb_mesh2/ED.msh'...
Info : 27 entities
Info : 24693 nodes
Info : 128345 elements
Info : [ 0%] Reading elements
Info : [ 10%] Reading elements
Info : [ 20%] Reading elements
Info : [ 30%] Reading elements
Info : [ 40%] Reading elements
Info : [ 50%] Reading elements
Info : [ 60%] Reading elements
Info : [ 70%] Reading elements
Info : [ 80%] Reading elements
Info : [ 90%] Reading 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 'ukb_mesh2/ED.msh'
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
rv: [2]
epi: [7]
base: [5, 6, 4, 3]
INFO:ldrb.ldrb: Num vertices: 24693
INFO:ldrb.ldrb: Num cells: 89853
INFO:ldrb.ldrb: Apex coord: (47.23, -19.64, -28.24)
INFO:ldrb.ldrb:
Calculating gradients
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 60
epi_lv: -60
endo_septum: 60
epi_septum: -60
endo_rv: 60
epi_rv: -60
INFO:ldrb.ldrb:beta:
endo_lv: 0
epi_lv: 0
endo_septum: 0
epi_septum: 0
endo_rv: 0
epi_rv: 0
2025-10-06 08:43:11 [debug ] Write f0: f0
2025-10-06 08:43:11 [debug ] Write s0: s0
2025-10-06 08:43:11 [debug ] Write n0: n0
2025-10-06 08:43:11 [debug ] Write lv: f
2025-10-06 08:43:11 [debug ] Write rv: f
2025-10-06 08:43:12 [debug ] Write epi: f
2025-10-06 08:43:12 [debug ] Write lv_rv: f
2025-10-06 08:43:12 [debug ] Write apex: f
2025-10-06 08:43:12 [debug ] Write lv_scalar: f
2025-10-06 08:43:12 [debug ] Write rv_scalar: f
2025-10-06 08:43:12 [debug ] Write epi_scalar: f
2025-10-06 08:43:12 [debug ] Write lv_rv_scalar: f
2025-10-06 08:43:12 [debug ] Write lv_gradient: f
2025-10-06 08:43:12 [debug ] Write rv_gradient: f
2025-10-06 08:43:12 [debug ] Write epi_gradient: f
2025-10-06 08:43:12 [debug ] Write apex_gradient: f
2025-10-06 08:43:12 [debug ] Write markers_scalar: f
INFO:cardiac_geometries.geometry:Reading geometry from ukb_mesh2