Cardiac CINE

In this notebook, a precomputed volume of the XCAT Phantom is used to select slices and evaluate a BSSFP signal model including off-resonance effects.

Imports

import sys
from typing import Tuple
sys.path.append('../../')
sys.path.append("../../tests/")

import tensorflow as tf
physical_devices = tf.config.list_physical_devices('GPU')
if physical_devices:
    print("Available GPUS: \n\t", "\n\t".join([str(_) for _ in physical_devices]))
    tf.config.experimental.set_visible_devices(physical_devices[1], 'GPU')
    tf.config.experimental.set_memory_growth(physical_devices[1], True)

from collections import OrderedDict
import numpy as np
import pyvista
from pint import Quantity
import math

import matplotlib.pyplot as plt
%matplotlib inline

import cmrseq
import cmrsim
import phantoms

Available GPUS:
         PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')
        PhysicalDevice(name='/physical_device:GPU:1', device_type='GPU')
        PhysicalDevice(name='/physical_device:GPU:2', device_type='GPU')

2022-08-15 09:59:15.218653: 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 AVX512F FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
2022-08-15 09:59:15.799503: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1532] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 22684 MB memory:  -> device: 1, name: NVIDIA TITAN RTX, pci bus id: 0000:3d:00.0, compute capability: 7.5

Define Simulation Class

class BSSFPSimulator(cmrsim.analytic.simulation.AnalyticSimulation):
    def build(self, fov: Tuple[float, float], matrix_size: Tuple[int, int], flip_angle: float, TE: float, TR: float, n_coils: int = 8):
        n_kx, n_ky = matrix_size

        # Encoding definition
        inplane_resolution = Quantity(fov, "m") / matrix_size
        sequence_list = cmrseq.seqdefs.sequences.balanced_ssfp(cmrseq.SystemSpec(), matrix_size, inplane_resolution,
                                                               slice_thickness=Quantity(10, "mm"),
                                                               adc_duration=Quantity(2, "ms"), flip_angle=Quantity(flip_angle, "degree"),
                                                               pulse_duration=Quantity(1, "ms"))
        encoding_module = cmrsim.analytic.encoding.GenericEncoding(name="bssfp_readout", sequence=sequence_list[1:], absolute_noise_std=0., device="GPU:0")

        # Signal module construction
        sequence_module = cmrsim.analytic.contrast.BSSFP(flip_angle=np.array([flip_angle, ]), echo_time=[TE, ], repetition_time=[TR, ], device="GPU:0")
        mdims = np.array([(-fov[0]/2, fov[0]/2), (-fov[1]/2, fov[1]/2), (0., 0.1)])
        coil_sens_module = cmrsim.analytic.contrast.BirdcageCoilSensitivities(map_shape=np.array((n_kx*4, n_ky*4, 1, 8)).astype(int),
                                                                              map_dimensions=mdims, device="GPU:0")

        # Forward model composition
        forward_model = cmrsim.analytic.CompositeSignalModel(sequence_module, coil_sens_module)
        return forward_model, encoding_module, None

Load input and compute the off-resonance map

from scipy.spatial.transform import Rotation

def select_slice(mesh: pyvista.UnstructuredGrid, slice_position: Quantity,
                 slice_normal: np.ndarray, slice_thickness: Quantity,
                 inplane_extend: Quantity, spacing: Quantity) -> pyvista.UnstructuredGrid:
    spatial_extend = np.array([*inplane_extend.m_as("m"), slice_thickness.m_as("m")])
    origin = -spatial_extend / 2
    dims = np.around(spatial_extend / spacing.m_as("m")).astype(int).tolist()
    slice_mesh = pyvista.UniformGrid(dims=dims,
                                     spacing=np.array(spacing.m_as("m")).tolist(),
                                     origin=np.array(origin).tolist())
    slice_mesh = slice_mesh.cast_to_structured_grid()
    # Rotate slice
    slice_normal = slice_normal / np.linalg.norm(slice_normal, keepdims=True)
    rot_vec = Rotation.align_vectors(np.array([[0., 0., 1.]]), slice_normal.reshape(1, 3))[0].as_rotvec().reshape(3)
    rot_angle = np.linalg.norm(rot_vec)
    print(rot_angle, rot_vec)
    prerot_r = slice_mesh.points
    slice_mesh = slice_mesh.rotate_vector(vector=rot_vec/rot_angle, angle=np.rad2deg(rot_angle), inplace=False)
    print(np.mean(slice_mesh.points - prerot_r))
    slice_mesh = slice_mesh.translate(slice_position.m_as("m"), inplace=False)
    return mesh.probe(slice_mesh)

data = pyvista.read("../../../itet-stor/mritrans/Stefano/Synthetic_image_input/Image_01.vti")
suszeptibilities = dict({0: 2.9e-7, 1: -6e-7, 2: -6e-7, 3: -6e-7, 4: -6e-7, 5: -4e-7, 6: -7e-7, 7: -4e-7, 8: -7e-7, 13: -4e-7, 15: 0., 16: 0, 36: -4e-7, 37: -7e-7})
labels = data["labels"].reshape(data.dimensions, order="F")
inh = cmrsim.utils.offresonance.generate_inhomogeneity_map(labels, Quantity(1, "T"), sus_maps=suszeptibilities, padsize=10)
data["inhomogeneity"] = inh.flatten(order = "F")
data

/scratch/jweine/cmrsim/notebooks/analytic_simulation/../../cmrsim/utils/offresonance.py:29: RuntimeWarning: invalid value encountered in true_divide
  corterm = 3 * z ** 2 / (x ** 2 + y ** 2 + z ** 2) - 1

Header

Data Arrays

UniformGrid Information N Cells 44390570 N Points 44826624 X Bounds -2.597e-01, 2.352e-01 Y Bounds -2.894e-01, 2.055e-01 Z Bounds -9.976e-02, 6.489e-02 Dimensions 512, 512, 171 Spacing 9.685e-04, 9.685e-04, 9.685e-04 N Arrays 11

Name Field Type N Comp Min Max labels Points float64 1 0.000e+00 6.600e+01 LV_mask Points int8 1 0.000e+00 1.000e+00 PD Points float64 1 3.250e+01 9.754e+01 T1 Points float64 1 2.685e+01 1.400e+03 T2 Points float64 1 4.152e+00 2.848e+02 T2s Points float64 1 1.000e+00 6.600e+01 e_circ Points int8 1 0.000e+00 0.000e+00 e_radial Points int8 1 0.000e+00 0.000e+00 e_long Points int8 1 0.000e+00 0.000e+00 magnitude Points float64 1 1.417e-02 1.214e+00 inhomogeneity Points float64 1 -3.613e+00 3.781e+00

Slice Volume to mimic slice-selective excitation

slice_position = Quantity([0., 0., -3.5], "cm")
slice_thickness = Quantity(10, "mm")
slice_normal = np.array([0., -0.3, 1.])
readout_direction = np.array([1., 0., 0.])
inplane_extend = Quantity([40, 40], "cm")
spacing = Quantity([0.5, 0.5, 2], "mm")

data_slice = select_slice(data, slice_position, slice_normal, slice_thickness,
                          inplane_extend=inplane_extend, spacing=spacing)
data_slice["labels"] = np.around(data_slice["labels"]).astype(np.int64)
data_slice["PD"][np.where(data_slice["labels"] == 0)] = 0.

/tmp/ipykernel_22762/794569812.py:15: UserWarning: Optimal rotation is not uniquely or poorly defined for the given sets of vectors.
  rot_vec = Rotation.align_vectors(np.array([[0., 0., 1.]]), slice_normal.reshape(1, 3))[0].as_rotvec().reshape(3)

0.29145679447786704 [-0.29145679  0.          0.        ]
-5.42645902337312e-05

Visualize input volume

pyvista.close_all()
pyvista.set_jupyter_backend("panel")
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(500, 750))
origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025,
                               tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1)
                 for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]]

# plotter.add_axes_at_origin()
plotter.add_mesh(data, scalars="T2", opacity=0.2, cmap="gray")
plotter.add_mesh(data_slice, scalars="inhomogeneity", opacity=0.5, cmap="twilight")

plotter.add_mesh(origin_arrows[0], color="blue", label="x")
plotter.add_mesh(origin_arrows[1], color="green", label="y")
plotter.add_mesh(origin_arrows[2], color="red", label="z")
plotter.show()

Set up dataset

Note: The coordinates of the input must be in mps/(RO-PE-SS) coordinates

r_mps = cmrseq.utils.xyz_to_mps(data_slice.points - slice_position.m_as("m"), slice_normal=slice_normal, readout_direction=readout_direction)
r_mps -= np.mean(r_mps, axis=0, keepdims=True)
cell_volume = Quantity(np.product(spacing.m_as("m")), "m**3")
data_slice["M0"] = data_slice["PD"] * cell_volume

phantom_dict = OrderedDict(r_vectors=r_mps, M0=np.array(data_slice["M0"], dtype=np.complex64),
                           T1=np.array(data_slice["T1"]), T2=np.array(data_slice["T2"]),
                           T2star=data_slice["T2s"], theta=data_slice["inhomogeneity"])

phantom_dict = {k: v.astype(np.float32) if k!="M0" else v for k, v in phantom_dict.items()}
dataset = cmrsim.datasets.AnalyticDataset(phantom_dict, filter_inputs=True, expand_dimension=True)

/scratch/jweine/conda/envs/tf29_pyvista/lib/python3.9/site-packages/pyvista/core/dataset.py:1877: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  scalars = np.array(scalars)

Instantiate simulator

matrix_size = np.array([128, 128])
fov = Quantity([40, 40], "cm")
simulator = BSSFPSimulator(build_kwargs=dict(fov=fov.m_as("m").astype(np.float32),
                                             matrix_size=matrix_size, flip_angle=35,
                                             TE=6, TR=12))
_ = simulator.encoding_module.k_space_segments.assign(32)

/scratch/jweine/conda/envs/tf29_pyvista/lib/python3.9/site-packages/pint/quantity.py:1313: RuntimeWarning: invalid value encountered in double_scalars
  magnitude = magnitude_op(new_self._magnitude, other._magnitude)

From lookup base (1, 512, 512, 1, 16)

f, a = plt.subplots(1, 1)
kvec, sampling_times = simulator.encoding_module._calculate_trajectory()
sc = a.scatter(kvec.numpy()[:, 0], kvec.numpy()[:, 1], c=sampling_times.numpy())
plt.colorbar(sc, label="Time from pulse [ms]")
f.suptitle("K-space")
a.set_ylabel("$k_y$"), a.set_xlabel("$k_x$")

(Text(0, 0.5, '$k_y$'), Text(0.5, 0, '$k_x$'))

Run Simulation

%time kspace = simulator(dataset, voxel_batch_size=100000, execute_eager=False, unstack_repetitions=True)

Run Simulation:         |---------------|0/2339073      |---------------|1/128

2022-08-15 10:00:42.551211: I tensorflow/compiler/xla/service/service.cc:170] XLA service 0x7efb44b2f1e0 initialized for platform CUDA (this does not guarantee that XLA will be used). Devices:
2022-08-15 10:00:42.551261: I tensorflow/compiler/xla/service/service.cc:178]   StreamExecutor device (0): NVIDIA TITAN RTX, Compute Capability 7.5
2022-08-15 10:00:42.592145: W tensorflow/compiler/tf2xla/kernels/assert_op.cc:38] Ignoring Assert operator linear_interpolation_lookup/assert_less_equal/Assert/Assert
2022-08-15 10:00:42.595002: I tensorflow/compiler/mlir/tensorflow/utils/dump_mlir_util.cc:263] disabling MLIR crash reproducer, set env var `MLIR_CRASH_REPRODUCER_DIRECTORY` to enable.
2022-08-15 10:00:43.488638: I tensorflow/core/platform/default/subprocess.cc:304] Start cannot spawn child process: No such file or directory
2022-08-15 10:00:43.581242: I tensorflow/compiler/jit/xla_compilation_cache.cc:478] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.

Run Simulation:         |XXXXXXXXXXXXXX-|2300000/2339073        |---------------|1/128

2022-08-15 10:01:54.130409: W tensorflow/compiler/tf2xla/kernels/assert_op.cc:38] Ignoring Assert operator linear_interpolation_lookup/assert_less_equal/Assert/Assert

Run Simulation:         |XXXXXXXXXXXXXXX|2339073/2339073        |XXXXXXXXXXXXXXX|128/128 CPU times: user 41.2 s, sys: 8.42 s, total: 49.7 s
Wall time: 1min 17s

Reconstruct Single Coil images

simulator.encoding_module.absolute_noise_std.assign(tf.constant([1e-8]))
simulator._update()
kspace_n = simulator.encoding_module.add_noise(kspace[..., 0, :])
channel_wise_recon = cmrsim.reconstruction.Cartesian2D(sample_matrix_size=matrix_size)(kspace_n)
print(channel_wise_recon.shape)

(8, 1, 128, 128)

plt.close("all")
f, a = plt.subplots(2, 4, figsize=(16, 8))
a = a.flatten()
[a[i].imshow(np.abs(ch[-1]), cmap="gray") for i, ch in enumerate(channel_wise_recon)]
[_.axis("off") for _ in a.flatten()]
f.tight_layout()