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#
[1]:
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#
[2]:
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#
[3]:
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)
[4]:
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
[4]:
Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Slice Volume to mimic slice-selective excitation#
[5]:
[6]:
/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#
[7]:
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()