Contracting LV-mesh

This notebook demonstrates the convenience functionality of organizing, handling and preparing dynamic meshes as simulation input.

Note: To keep the notebook memory-size reasonable, to animations are reduced to only a short part of the contracting cylce.

Imports

import os
import math
from typing import List, Union, Tuple, Sequence
import sys
import base64
from IPython.display import display, HTML

import pyvista
import numpy as np
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
import vtk
from pint import Quantity

sys.path.insert(0, "../../")
import cmrsim

Load mesh from list of files

The following cell assumes, that the motion states are stored as separate vtk files in a separate folder and the corresponding time-stamps for each motion state are listed in a text-file. Furthermore, the LV-mesh is assumed to build in shells (4), rings (40 points per ring) and slices(30) with continous node-ordering according to:

    [apex, n_points_per_ring(1st slice, 1st shell), n_points_per_ring(2st slice, 1st shell),
    ..., n_points_per_ring(1st slice, 2nd shell)]

shell --> radial
slice --> longitudinal

files = [f'../example_resources/mesh_displacements/Displ{i:03}.vtk' for i in range(11)]
timing = Quantity(np.loadtxt("../example_resources/mesh_displacements/PV_loop_reordered.txt")[:, 1], "ms")
mesh_dataset = cmrsim.datasets.CardiacMeshDataset.from_list_of_meshes(files, timing[:11], time_precision_ms=3, mesh_numbers=(4, 40, 30))

Refine dataset

To avoid inverse crime during simulation, the mesh should be refined. This means the mesh nodes are linearly interpolated and the cell connectivity is evaluated for all new mesh-nodes

refined_mesh = mesh_dataset.refine(longitudinal_factor=2, circumferential_factor=2, radial_factor=2)
refined_mesh.mesh.save(f"../example_resources/refined_mesh_{tuple(refined_mesh.mesh_numbers[i] for i in (2, 0, 1))}.vtk")

Load from single file

When the displacements are already stored correctly in a single vtk file you call also load it from there

fn = "../example_resources/refined_mesh_(7, 81, 60).vtk"
mesh_numbers = [int(s) for s in os.path.basename(fn).replace("refined_mesh_(", "").replace(").vtk", "").split(",")]
mesh_dataset_refined = cmrsim.datasets.CardiacMeshDataset.from_single_vtk(f"../resources/{fn}", time_precision_ms=3, mesh_numbers=mesh_numbers)

Render Animation

To check the mesh, and its stored field arrays by rendering an animation of the contracting LV. Only the file-name is mandatory, for a more detailed overview look into the API reference for the CardiacMeshDataset.

custom_theme = pyvista.themes.DefaultTheme()
custom_theme.background = 'white'
custom_theme.font.family = 'times'
custom_theme.font.color = 'black'
custom_theme.outline_color = "white"
custom_theme.edge_color = 'grey'
custom_theme.transparent_background = True

pyvista.close_all()
# pyvista.start_xvfb()  # Enable for Linux
plotter = pyvista.Plotter(off_screen=True, theme=custom_theme, window_size=(600, 600))
mesh_dataset.render_input_animation(filename="original_mesh", plotter=plotter, mesh_kwargs=dict(opacity=0.8, show_edges=True))

plotter = pyvista.Plotter(off_screen=True, theme=custom_theme, window_size=(600, 600))
mesh_dataset_refined.render_input_animation(filename="refined_mesh", plotter=plotter, mesh_kwargs=dict(opacity=0.8, show_edges=True))

import imageio

gif1 = imageio.v3.imread(f"original_mesh.gif")
gif2 = imageio.v3.imread(f"refined_mesh.gif")
combined_gif = imageio.get_writer(f'meshmesh.gif')

for f1, f2 in tqdm(zip(gif1, gif2), total=gif1.shape[0]):
    combined_gif.append_data(np.concatenate([f1, f2], axis=1))
combined_gif.close()

b64 = base64.b64encode(open("meshmesh.gif",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))

Select Slice

Usually only a subset of the entire mesh in necessary. The dataset offers a convenient way to select a slice from the refined mesh.

slice_dict = dict(slice_thickness=Quantity(10, "mm"),
                  slice_normal=np.array((0., 0., 1.)),
                  slice_position=Quantity([0, 0, -3], "cm"),
                  reference_time=Quantity(0, "ms"))

slice_mesh = mesh_dataset_refined.select_slice(**slice_dict)
slice_mod =  cmrsim.datasets.MeshDataset(slice_mesh, mesh_dataset_refined.timing)

plotter = pyvista.Plotter(off_screen=True, theme=custom_theme, window_size=(600, 600))
slice_mod.render_input_animation(filename="slice_mesh", plotter=plotter, mesh_kwargs=dict(opacity=0.8, show_edges=True))

b64 = base64.b64encode(open("slice_mesh.gif",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))

Retreiving information and field data

The MeshDataset implemetns a couple of convenient methods to get information as numpy arrays. For a full list refer to the API documentation. The following cells are meant to demonstrate only a few.

t, r = slice_mod.get_trajectories(start=slice_dict["reference_time"], end=Quantity(slice_mod.timing[-1], "ms"))
t.shape, r.shape

((10,), (4658, 10, 3))

r = slice_mod.get_positions_at_time(Quantity(0., "ms"))
r.shape

(4658, 3)

t, disps = slice_mod.get_field("displacements", range(5))
t, disps.shape

(array([ 0.  ,  4.04,  6.06,  8.08, 10.1 ]) <Unit('millisecond')>,
 (5, 4658, 3))