Proper orthogonal decomposition (POD) trajectories

The POD trajectory module computes the proper orthogonal decomposition for mesh data of shape (#Mesh nodes, #Time-steps, #Channels) and fit a taylor expansion to weight-functions. In this notebook the use-case of a contracting LV mesh including fiber and sheet orientation is demonstrated.

Import

from IPython.display import display, HTML, clear_output
import base64
import imageio

import tensorflow as tf
# gpu = tf.config.get_visible_devices("GPU")[0]
# tf.config.set_visible_devices(gpu, device_type="GPU")
# tf.config.experimental.set_memory_growth(gpu, True)

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

import os
import math
from typing import List, Union, Tuple, Sequence
from pint import Quantity
from scipy.optimize import curve_fit

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

sys.path.append("../")
import local_functions

2022-12-20 13:40:44.352780: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
2022-12-20 13:40:47.264574: E tensorflow/compiler/xla/stream_executor/cuda/cuda_driver.cc:267] failed call to cuInit: UNKNOWN ERROR (34)
2022-12-20 13:40:47.264629: I tensorflow/compiler/xla/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (0fc0ac3d5ebf): /proc/driver/nvidia/version does not exist
2022-12-20 13:40:47.264877: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.

Load cardiac mesh

files = [f'../../../cmr-random-diffmaps/notebooks/example_data/example_mesh/Displacements/Displ{i:03}.vtk' for i in range(0, 200, 1)]
timing = Quantity(np.loadtxt("../../../cmr-random-diffmaps/notebooks/example_data/example_mesh/PV_loop_reordered.txt")[:, 1], "ms")

original_mesh = cmrsim.datasets.CardiacMeshDataset.from_list_of_meshes(files, timing[:len(files)], mesh_numbers=(4, 40, 30), time_precision_ms=3)

Explanation/Exploration of the POD results

To give an insight into how the POD module works, the first section demonstrates the underlying computation and highlights intermediate results.

Caclulate POD

t, r = original_mesh.get_trajectories(start=Quantity(0, "ms"), end=Quantity(5000, "ms"))
n_modes=5
phi, weights = cmrsim.trajectory.PODTrajectoryModule.calculate_pod(t.m_as("ms"), data=r.m_as("m"), n_modes=n_modes)
print(phi.shape, weights.shape)

(14412, 5) (199, 5)

Visualize the modes of computed POD

pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(off_screen=True, window_size=(1600, 400), shape=(1, n_modes))
for n in range(n_modes):
    plotter.subplot(0, n)
    plotter.add_title(f"Mode {n}", font_size=12)
    plotter.add_mesh(phi.T[n].reshape(-1, 3), render_points_as_spheres=True)

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

Plot mode contributions over time and render the reconstructed mesh

Reconstructing a motion state \(u(t)\) of the mesh is done by calculated the weighted sum of the mesh modes for each time point:

\[u(t) = \Sigma_j^{N_{modes}} \phi_j w_j(t),\]

where \(\phi_j\) are the computed basis functions (modes) shown above and \(w_j(t)\) are the corresponding mode-weights as funtion of time shown below. The motion states reconstructed by 5 modes is shown as animation below.

# Plot mode weights as function of time
fig, ax = plt.subplots(1, 1, figsize=(12, 5), dpi=100)
ax.plot(t.m, weights);
ax.set_xlabel("Time (ms)"), ax.grid(alpha=0.5), ax.set_title("Mode contributions over time")
ax.legend(["$\phi_{i}$".format(i=n) for n in range(0, phi.shape[1])], ncol=5)
fig.tight_layout()
fig.savefig("modecontributions.png")
plt.close("all")
clear_output()

# Render reconstructed mesh
t_indices = np.linspace(0, t.shape[0]-1, 50).astype(int)
phi_dash = phi.T.reshape(n_modes, -1, 3)
trajectories = np.stack([np.sum((phi_dash * weights[t_idx].reshape(-1, 1, 1)), axis=0) for t_idx in tqdm(t_indices)])
plotter = pyvista.Plotter(off_screen=False, window_size=(500, 500), shape=(1, 1))
local_functions.animate_trajectories(plotter, trajectories,
                                     timing=[t.m_as("ms")[_] for _ in t_indices],
                                     filename="PODRecon.gif",
                                     mesh_kwargs=dict(render_points_as_spheres=True),
                                     visualize_trajectories=False,
                                     trajectory_kwargs=dict(opacity=0.5, color="red"))
pyvista.close_all()

# Combine figures and display
gif1 = imageio.v3.imread("PODRecon.gif")
img = imageio.v3.imread("modecontributions.png")
combined_gif = imageio.get_writer("meshrecon.gif")
for f1 in tqdm(gif1, total=gif1.shape[0]):
    combined_gif.append_data(np.concatenate([img[..., :3], f1], axis=1))
combined_gif.close()
clear_output()
b64 = base64.b64encode(open("meshrecon.gif",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))

Calling the POD module

In the this second section the POD module is instantiated and called, where the mechanism above is used internally on instantiating the module. The mode-weights are represented by a Taylor-series of specified order, which allows to reconstruct the motion state arbitrary times within the time-interval of given mesh-snapshots.

tf.config.run_functions_eagerly(True)
# Get mesh-snap shots
t, r = original_mesh.get_trajectories(start=Quantity(0, "ms"), end=Quantity(5000, "ms"))

# Define number of modes and degree of Taylor series to represent the mode-weights
n_modes=5
poly_order=10

# Instantiate module
module = cmrsim.trajectory.PODTrajectoryModule(t.m_as("ms"), data=r.m_as("m"), n_modes=n_modes, poly_order=poly_order)
# Or do batched calls:
# module = cmrsim.trajectory.PODTrajectoryModule(t.m_as("ms"), data=r.m_as("m"), n_modes=n_modes, poly_order=poly_order, batch_size=1000)

# Define new time-points
new_time_grid = np.linspace(t[0].m, t[-1].m, 100).astype(np.float32)

# Get reconstructed mesh at new time-points
reconstructed_mesh_points, _ = module(new_time_grid)

# (batched call continued)
reconstructed_mesh_points = tf.concat([module(new_time_grid, batch_index=i)[0] for i in range(10)], axis=0)

clear_output()
reconstructed_mesh_points.shape

TensorShape([4804, 100, 3])

By accessing the super-class _evaluate_trajectory method (this is only for demonstration purpose) we can also get the mode-weights at the new timings

mode_weights = tf.transpose(module._taylor_module._evaluate_trajectory(new_time_grid)[:, :, 0], [1, 0])

Plot the results of reconstruction

# Plot mode weights as function of time
fig, ax = plt.subplots(1, 1, figsize=(12, 5), dpi=100)
ax.plot(new_time_grid, mode_weights);
ax.set_xlabel("Time (ms)"), ax.grid(alpha=0.5), ax.set_title("Mode contributions over time (module call)")
ax.legend(["$\phi_{i}$".format(i=n) for n in range(0, phi.shape[1])], ncol=5)
fig.tight_layout()
fig.savefig("modecontributions.png")
plt.close("all")
clear_output()

# Render reconstructed mesh
plotter = pyvista.Plotter(off_screen=False, window_size=(500, 500), shape=(1, 1))
local_functions.animate_trajectories(plotter, np.swapaxes(reconstructed_mesh_points.numpy(), 0, 1),
                                     timing=new_time_grid,
                                     filename="PODRecon.gif",
                                     mesh_kwargs=dict(render_points_as_spheres=True),
                                     visualize_trajectories=False,
                                     trajectory_kwargs=dict(opacity=0.5, color="red"))
pyvista.close_all()

# Combine figures and display
gif1 = imageio.v3.imread("PODRecon.gif")
img = imageio.v3.imread("modecontributions.png")
combined_gif = imageio.get_writer("meshrecon.gif")
for f1 in tqdm(gif1, total=gif1.shape[0]):
    combined_gif.append_data(np.concatenate([img[..., :3], f1], axis=1))
combined_gif.close()
clear_output()
b64 = base64.b64encode(open("meshrecon.gif",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))