Storing Additional Data in the Geometry#
This demo illustrates how to extend the standard Geometry storage format to include
custom fields. We use a Bi-Ventricular (BiV) mesh from the UK Biobank Atlas as an example.
We often need more than just the standard microstructure (\(\mathbf{f}_0, \mathbf{s}_0, \mathbf{n}_0\)). Common requirements include:
AHA Segments: The 17-segment model for regional analysis.
Auxiliary Vectors: Purely circumferential (\(\mathbf{c}_0\)) or longitudinal (\(\mathbf{l}_0\)) directions for strain analysis.
Custom Markers: Patient-specific regions of interest (e.g., scars).
We will generate these fields using cardiac-geometries and fenicsx-ldrb, save them
alongside the mesh, and demonstrate how to reload them.
Imports#
import shutil
from pathlib import Path
import numpy as np
from mpi4py import MPI
import dolfinx
import ldrb
import scifem
import cardiac_geometries as cg
1. Geometry Generation#
We generate a mesh from the UK Biobank (UKB) statistical atlas. See UKB Atlas Documentation for details.
comm = MPI.COMM_WORLD
mode = -1 # Use the mean mode
std = 0
char_length = 5.0
outdir = Path("ukb_atlas_additional_data")
geodir = outdir / "geometry"
# Clean up previous geometry if it exists
if geodir.exists():
shutil.rmtree(geodir)
geo = cg.mesh.ukb(
outdir=geodir,
comm=comm,
mode=mode,
std=std,
case="ED",
create_fibers=True,
char_length_max=char_length,
char_length_min=char_length,
clipped=True,
fiber_angle_endo=60,
fiber_angle_epi=-60,
fiber_space="Quadrature_6",
)
INFO:ukb.atlas:Downloading https://github.com/ComputationalPhysiology/ukb-atlas/releases/download/atlas-data/UKBRVLV.h5 to /github/home/.ukb/UKBRVLV.h5. This may take a while.
INFO:ukb.atlas:Done downloading.
INFO:ukb.atlas:Generating points from /github/home/.ukb/UKBRVLV.h5
INFO:ukb.atlas:Using mode -1 and std 0.0
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/LV_ENDO_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/RV_SEPTUM_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/RV_FREEWALL_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/LV_EPI_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/RV_EPI_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/MV_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/AV_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/TV_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/PV_ED.stl
INFO:ukb.surface:Saved ukb_atlas_additional_data/geometry/landmarks_ED.json
INFO:ukb.clip:Folder: ukb_atlas_additional_data/geometry
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 ukb_atlas_additional_data/geometry/LV_ENDO_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 ukb_atlas_additional_data/geometry/lv_clipped.ply
INFO:ukb.clip:Reading ukb_atlas_additional_data/geometry/RV_SEPTUM_ED.stl, ukb_atlas_additional_data/geometry/RV_FREEWALL_ED.stl
INFO:ukb.clip:Smoothing RV
INFO:ukb.clip:Saving ukb_atlas_additional_data/geometry/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 ukb_atlas_additional_data/geometry/LV_EPI_ED.stl, ukb_atlas_additional_data/geometry/RV_EPI_ED.stl
INFO:ukb.clip:Saving ukb_atlas_additional_data/geometry/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.
INFO:ukb.mesh:Creating clipped mesh for ED with char_length_max=5.0, char_length_min=5.0, order=1
0
INFO:ukb.mesh:Created mesh ukb_atlas_additional_data/geometry/ED_clipped.msh
2026-11-12 21:36:34 [debug ] Convert file ukb_atlas_additional_data/geometry/ED_clipped.msh to dolfin
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
epi: [3]
base: [4]
rv: [2]
INFO:ldrb.ldrb: Num vertices: 2593
INFO:ldrb.ldrb: Num cells: 8054
Info : Reading 'ukb_atlas_additional_data/geometry/ED_clipped.msh'...
Info : 11 entities
Info : 2593 nodes
Info : 12948 elements
Info : 3 parametrizations
Info : [ 0%] Processing parametrizations
Info : [ 10%] Processing parametrizations
Info : [ 40%] Processing parametrizations
Info : Done reading 'ukb_atlas_additional_data/geometry/ED_clipped.msh'
INFO:ldrb.ldrb: Apex coord: (47.76, -18.06, -29.37)
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
# Visualize the standard fiber field
try:
import pyvista
except ImportError:
print("Pyvista is not installed")
else:
plotter = pyvista.Plotter()
# Only show every 60th point to avoid clutter (and keep the size of the docs down). The
# quadrature space has 24 points per cell, six times as many as the DG 1 space below
N = 60
assert geo.f0 is not None
points = geo.f0.function_space.tabulate_dof_coordinates()
point_cloud = pyvista.PolyData(points[::N, :])
f0_arr = geo.f0.x.array.reshape((points.shape[0], 3))
point_cloud["fibers"] = f0_arr[::N, :]
fibers = point_cloud.glyph(
orient="fibers",
scale=False,
factor=5.0,
)
plotter.add_mesh(fibers, color="red")
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(outdir / "f0.png")
INFO:trame_server.utils.namespace:Translator(prefix=None)
INFO:trame_rca.encoders.image_encoder:Turbo JPEG - NOT AVAILABLE (Python package)
INFO:trame_rca.encoders:VTKStreaming Video encoding is NOT AVAILABLE (missing Python package)
INFO:trame_rca.schedulers:cannot import name 'RcaVideoEncoder' from 'trame_rca.encoders' (/dolfinx-env/lib/python3.12/site-packages/trame_rca/encoders/__init__.py)
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 timeout 0
INFO:wslink.backends.aiohttp:awaiting running future
2. Generating AHA Segments#
We use fenicsx-ldrb to generate the 17-segment American Heart Association (AHA) model.
This assigns an integer ID to each cell in the mesh corresponding to its region (Basal Anterior, Septal, etc.).
# Generate AHA function (DG0 space)
aha = ldrb.aha.gernerate_aha_biv(
mesh=geo.mesh,
ffun=geo.ffun,
markers=cg.mesh.transform_markers(geo.markers, clipped=True),
function_space="DG_0",
base_max=0.75,
mid_base=0.70,
apex_mid=0.65,
)
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
epi: [3]
base: [4]
rv: [2]
INFO:ldrb.ldrb: Num vertices: 2593
INFO:ldrb.ldrb: Num cells: 8054
INFO:ldrb.ldrb: Apex coord: (47.76, -18.06, -29.37)
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: 90
epi_rv: -25
INFO:ldrb.ldrb:beta:
endo_lv: -20
epi_lv: 20
endo_septum: -20
epi_septum: 20
endo_rv: -20
epi_rv: 20
INFO:ldrb.aha:Basal region: 4196 inds between 0.938978 and 1.006047
INFO:ldrb.aha:Mid region: 2167 inds between 0.871908 and 0.938978
# Convert the DG0 function to dolfinx.mesh.MeshTags
# This is often more convenient for post-processing and sub-domain integration.
entities = np.hstack(comm.allgather(np.arange(*geo.mesh.topology.index_map(3).local_range, dtype=np.int32)))
values = np.hstack(comm.allgather(aha.x.array.astype(np.int32)))
aha_mt = dolfinx.mesh.meshtags(
geo.mesh, 3, entities, values,
)
# Visualize AHA segments
try:
import pyvista
except ImportError:
pass
else:
vtk_mesh = dolfinx.plot.vtk_mesh(geo.mesh, geo.mesh.topology.dim)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
grid.cell_data["AHA"] = aha_mt.values
grid.set_active_scalars("AHA")
p = pyvista.Plotter(window_size=[800, 800])
p.add_mesh(grid, show_edges=True, cmap="tab20")
if not pyvista.OFF_SCREEN:
p.show()
else:
p.screenshot(outdir / "aha_segments.png")
Note however that the markers near the base is a little bit off in this case.
3. Generating Custom Fiber Vectors#
For strain analysis, it is useful to have vector fields representing the purely Circumferential (\(\mathbf{c}_0\)) and Longitudinal (\(\mathbf{l}_0\)) directions, separate from the helical myofibers.
We can generate these by running ldrb with specific angles (\(0^\circ\) and \(90^\circ\)).
fiber_space = "DG_1"
Circumferential Vectors (0 degrees)#
system_c = ldrb.dolfinx_ldrb(
mesh=geo.mesh,
ffun=geo.ffun,
markers=cg.mesh.transform_markers(geo.markers, clipped=True),
alpha_endo_lv=0,
alpha_epi_lv=0,
fiber_space=fiber_space,
)
c0 = system_c.f0
c0.name = "c0"
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
epi: [3]
base: [4]
rv: [2]
INFO:ldrb.ldrb: Num vertices: 2593
INFO:ldrb.ldrb: Num cells: 8054
INFO:ldrb.ldrb: Apex coord: (47.76, -18.06, -29.37)
INFO:ldrb.ldrb:
Calculating gradients
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 0
epi_lv: 0
endo_septum: 0
epi_septum: 0
endo_rv: 0
epi_rv: 0
INFO:ldrb.ldrb:beta:
endo_lv: -65
epi_lv: 25
endo_septum: -65
epi_septum: 25
endo_rv: -65
epi_rv: 25
# Visualize Circumferential Vectors
try:
import pyvista
except ImportError:
pass
else:
plotter = pyvista.Plotter()
N = 10 # Only show every 10th point to avoid clutter
points = c0.function_space.tabulate_dof_coordinates()
point_cloud = pyvista.PolyData(points[::N, :])
c0_arr = c0.x.array.reshape((points.shape[0], 3))
point_cloud["fibers"] = c0_arr[::N, :]
fibers = point_cloud.glyph(orient="fibers", scale=False, factor=5.0)
plotter.add_mesh(fibers, color="blue")
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(outdir / "c0.png")
Longitudinal Vectors (90 degrees)#
system_l = ldrb.dolfinx_ldrb(
mesh=geo.mesh,
ffun=geo.ffun,
markers=cg.mesh.transform_markers(geo.markers, clipped=True),
alpha_endo_lv=90,
alpha_epi_lv=90,
epi_only=True, # Optimization: Longitudinal direction is roughly constant transmurally
fiber_space=fiber_space,
)
l0 = system_l.f0
l0.name = "l0"
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
epi: [3]
base: [4]
rv: [2]
INFO:ldrb.ldrb: Num vertices: 2593
INFO:ldrb.ldrb: Num cells: 8054
INFO:ldrb.ldrb: Apex coord: (47.76, -18.06, -29.37)
INFO:ldrb.ldrb:
Calculating gradients
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 90
epi_lv: 90
endo_septum: 90
epi_septum: 90
endo_rv: 90
epi_rv: 90
INFO:ldrb.ldrb:beta:
endo_lv: -65
epi_lv: 25
endo_septum: -65
epi_septum: 25
endo_rv: -65
epi_rv: 25
# Visualize Longitudinal Vectors
try:
import pyvista
except ImportError:
pass
else:
plotter = pyvista.Plotter()
N = 10 # Only show every 10th point to avoid clutter
points = l0.function_space.tabulate_dof_coordinates()
point_cloud = pyvista.PolyData(points[::N, :])
l0_arr = l0.x.array.reshape((points.shape[0], 3))
point_cloud["fibers"] = l0_arr[::N, :]
fibers = point_cloud.glyph(orient="fibers", scale=False, factor=5.0)
plotter.add_mesh(fibers, color="blue")
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(outdir / "l0.png")
4. Saving Geometry with Additional Data#
We now bundle our custom fields into a dictionary and pass it to save_geometry.
This function serializes the mesh, standard markers, standard fibers, and our additional data
into a consistent folder structure (using HDF5/XDMF for mesh tags and ADIOS2 for functions).
# Overwrite the existing geometry.bp folder to include the new data
# Note: We use the existing mesh/markers to ensure consistency and delete the old data first
shutil.rmtree(geodir / "geometry.bp")
cg.geometry.save_geometry(
path=geodir / "geometry.bp",
mesh=geo.mesh,
ffun=geo.ffun,
markers=geo.markers,
info=geo.info,
f0=geo.f0,
s0=geo.s0,
n0=geo.n0,
additional_data=additional_data,
)
5. Reloading and Verification#
We load the geometry from the folder. The Geometry object will now populate its
additional_data attribute with our custom fields.
geo_loaded = cg.geometry.Geometry.from_folder(
comm=comm,
folder=geodir,
)
INFO:cardiac_geometries.geometry:Reading geometry from ukb_atlas_additional_data/geometry
print("Loaded keys:", geo_loaded.additional_data.keys())
Loaded keys: dict_keys(['aha_segments', 'c0', 'l0'])
# Verify AHA Segments
try:
import pyvista
except ImportError:
pass
else:
# Access the loaded MeshTags
loaded_aha = geo_loaded.additional_data["aha_segments"]
vtk_mesh = dolfinx.plot.vtk_mesh(geo_loaded.mesh, geo_loaded.mesh.topology.dim)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
grid.cell_data["AHA"] = loaded_aha.values
grid.set_active_scalars("AHA")
p = pyvista.Plotter(window_size=[800, 800])
p.add_mesh(grid, show_edges=True, cmap="tab20")
if not pyvista.OFF_SCREEN:
p.show()
else:
p.screenshot(outdir / "loaded_aha_segments.png")
# Verify Circumferential Vectors
try:
import pyvista
except ImportError:
pass
else:
# Access the loaded Function
loaded_c0 = geo_loaded.additional_data["c0"]
assert isinstance(loaded_c0, dolfinx.fem.Function)
plotter = pyvista.Plotter()
N = 10 # Only show every 10th point to avoid clutter
points = loaded_c0.function_space.tabulate_dof_coordinates()
point_cloud = pyvista.PolyData(points[::N, :])
c0_arr = loaded_c0.x.array.reshape((points.shape[0], 3))
point_cloud["fibers"] = c0_arr[::N, :]
fibers = point_cloud.glyph(orient="fibers", scale=False, factor=5.0)
plotter.add_mesh(fibers, color="green")
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(outdir / "loaded_c0.png")
# Verify Longitudinal Vectors
try:
import pyvista
except ImportError:
pass
else:
# Access the loaded Function
loaded_l0 = geo_loaded.additional_data["l0"]
assert isinstance(loaded_l0, dolfinx.fem.Function)
plotter = pyvista.Plotter()
N = 10 # Only show every 10th point to avoid clutter
points = loaded_l0.function_space.tabulate_dof_coordinates()
point_cloud = pyvista.PolyData(points[::N, :])
f0_arr = loaded_l0.x.array.reshape((points.shape[0], 3))
point_cloud["fibers"] = f0_arr[::N, :]
fibers = point_cloud.glyph(orient="fibers", scale=False, factor=5.0)
plotter.add_mesh(fibers, color="green")
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(outdir / "loaded_l0.png")
Also make sure the fiber fields still looks good
# Visualize the standard fiber field
try:
import pyvista
except ImportError:
print("Pyvista is not installed")
else:
plotter = pyvista.Plotter()
# Only show every 60th point to avoid clutter (and keep the size of the docs down). The
# quadrature space has 24 points per cell, six times as many as the DG 1 space below
N = 60
assert geo_loaded.f0 is not None
points = geo_loaded.f0.function_space.tabulate_dof_coordinates()
point_cloud = pyvista.PolyData(points[::N, :])
f0_arr = geo_loaded.f0.x.array.reshape((points.shape[0], 3))
point_cloud["fibers"] = f0_arr[::N, :]
fibers = point_cloud.glyph(
orient="fibers",
scale=False,
factor=5.0,
)
plotter.add_mesh(fibers, color="red")
if not pyvista.OFF_SCREEN:
plotter.show()
else:
plotter.screenshot(outdir / "f0.png")