Turbulent flow with velocity encoded spoiled GRE#
Imports#
[ ]:
!pip install -U cmrseq --index-url https://gitlab.ethz.ch/api/v4/projects/30300/packages/pypi/simple
[3]:
import sys
sys.path.append("../")
sys.path.insert(0, "../../")
# sys.path.insert(0, "../../../cmrseq")
import os
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
os.environ['CUDA_VISIBLE_DEVICES'] = "3"
import tensorflow as tf
gpu = tf.config.get_visible_devices("GPU")
if gpu:
tf.config.experimental.set_memory_growth(gpu[0], True)
[4]:
from IPython.display import display, clear_output, Image
import pyvista
from pint import Quantity
import numpy as np
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
%matplotlib inline
import cmrsim
import cmrseq
sys.path.append("../")
import local_functions
clear_output()
Sequence Definition#
[5]:
# Define MR-system specifications
system_specs = cmrseq.SystemSpec(max_grad=Quantity(80, "mT/m"),
max_slew=Quantity(200., "mT/m/ms"),
grad_raster_time=Quantity(0.01, "ms"),
rf_raster_time=Quantity(0.01, "ms"),
adc_raster_time=Quantity(0.01, "ms"))
# Sequence Parameters (deliberately set too short --> get shortest)
TR = Quantity(5, 'ms')
TE = Quantity(0, 'ms')
# Define Fourier-Encoding parameters
fov = Quantity([222, 142], "mm")
matrix_size = (111, 71)
res = fov / matrix_size
print("Resolution:", res)
# Define imaging slice parameters
slice_thickness = Quantity(5, "mm")
slice_position = Quantity([0., 0., 0.], "m")
adc_duration = Quantity(4., "ms")
readout_direction = np.array([0., 0., 1.]) # readout along scanner Z, (FH)
slice_normal = np.array([1., 0., 0.])
slice_position_offset = Quantity(0, "cm")
rotation_matrix_to_mps = cmrseq.utils.get_rotation_matrix(slice_normal, readout_direction, target_orientation='mps').T
rotation_matrix_to_xyz = cmrseq.utils.get_rotation_matrix(slice_normal, readout_direction, target_orientation='xyz').T
print("MPS-rotation matrix to mps: \n", *rotation_matrix_to_mps)
print("MPS-rotation matrix to xyz: \n", *rotation_matrix_to_xyz)
# Define RF-exication parameters (coupled with seeding slice atm)
flip_angle = Quantity(np.pi / 12, "rad")
pulse_duration = Quantity(1, "ms")
RF_slice_normal = np.array([1., 0., 0.]) # Slice along scanner Y, (AP)
RF_slice_position = Quantity([0., 0., 0.], "m")
time_bandwidth_product = 4
# Define crusher area
crusher_duration = Quantity(1, "ms")
crusher_area = Quantity(np.pi * 4, "rad") / system_specs.gamma_rad / res[0]
# Velocity encoding settings:
venc_list = [Quantity(0., 'm/s'), Quantity(5.5, 'm/s'), Quantity(0.5, 'm/s'),
Quantity(2.5, 'm/s'), Quantity(0.5, 'm/s'), Quantity(1, 'm/s'), Quantity(0.5, 'm/s')]
all_venc_directions = [np.array([0., 0., 1.]), *np.repeat(np.eye(3, 3)[::-1, :], 2, axis=0)]
print("Venc directions:\n", all_venc_directions)
venc_duration = Quantity(1, 'ms') # Duration = 0 will result in shortest possible venc
all_venc_directions_MPS = [rotation_matrix_to_mps@d for d in all_venc_directions]
print("Rotated vencs:\n", all_venc_directions_MPS)
Resolution: [2.0 2.0] millimeter
MPS-rotation matrix to mps:
[ 0. -0. 1.] [ 0. -1. 0.] [ 1. -0. 0.]
MPS-rotation matrix to xyz:
[0. 0. 1.] [ 0. -1. 0.] [1. 0. 0.]
Venc directions:
[array([0., 0., 1.]), array([0., 0., 1.]), array([0., 0., 1.]), array([0., 1., 0.]), array([0., 1., 0.]), array([1., 0., 0.]), array([1., 0., 0.])]
Rotated vencs:
[array([1., 0., 0.]), array([1., 0., 0.]), array([1., 0., 0.]), array([ 0., -1., 0.]), array([ 0., -1., 0.]), array([0., 0., 1.]), array([0., 0., 1.])]
Each list of sequence objects in the variable multi_sequence_list corresponds to all repetitions for one VENC. The sequence illustration tries to highlight this structure.
[6]:
multi_sequence_list: list[list[cmrseq.Sequence]] = cmrseq.contrib.pc_gre_multivenc(
system_specs, slice_thickness=slice_thickness,
flip_angle=flip_angle, pulse_duration=pulse_duration,
time_bandwidth_product=time_bandwidth_product,
matrix_size=matrix_size, inplane_resolution=res,
adc_duration=adc_duration, echo_time=TE,
repetition_time=TR, slice_position_offset=slice_position_offset,
venc_list=venc_list, venc_direction_list=all_venc_directions,
venc_duration=venc_duration, crusher_area=crusher_area,
crusher_duration=crusher_duration)
# Velocity encoding is already defined in XYZ, so only rotate imaging gradients to XYZ
for list_per_venc in multi_sequence_list:
for seq_per_tr in list_per_venc:
tmp_seq = seq_per_tr.partial_sequence(copy_blocks=False, partial_string_match="velocity", invert_pattern=True)
tmp_seq.rotate_gradients(rotation_matrix_to_xyz)
clear_output()
[8]:
# This animation is created in the appendix of this notebook
display(Image("mr_sequence_animation.gif"))
<IPython.core.display.Image object>
Phantom preparation#
Load U-bend mesh#
[9]:
mesh = pyvista.read("../example_resources/sUbend_219_9_rotated.vtu")
mesh
[9]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Create 3D Rendering of mean velocity and the trace of Reynolds stress tensor on the meshed domain#
[15]:
_tmpmesh = mesh.copy()
_tmpmesh["abs_vel"] = np.linalg.norm(mesh["velocity"], axis=1)
_tmpmesh["trace_rst"] = np.mean(mesh["RST"][:, 0:3], axis=1)
pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(off_screen=True, window_size=(800, 400), shape=(1, 2))
plotter.subplot(0, 0)
plotter.add_mesh(_tmpmesh, scalars="abs_vel", scalar_bar_args=dict(title=r"$\| \vec{v} \|$"), cmap="twilight", opacity=0.5)
plotter.add_mesh(_tmpmesh.slice_orthogonal(), scalars="abs_vel", show_scalar_bar=False, cmap="twilight")
plotter.subplot(0, 1)
plotter.add_mesh(_tmpmesh, scalars="trace_rst", scalar_bar_args=dict(title="Tr(RST)"), cmap="rainbow", opacity=0.5)
plotter.add_mesh(_tmpmesh.slice_orthogonal(), scalars="trace_rst", show_scalar_bar=False, cmap="rainbow")
plotter.screenshot("input_mesh.png");
display(Image("input_mesh.png"))
Define Seeding volumes and options#
[16]:
# Define Reseeding slice parameters
seeding_slice_position = Quantity([0, 0.048, -0.045], "m")
seeding_slice_normal = np.array([1., 0., 0.])
seeding_slice_thickness = Quantity(50, "mm")
seeding_pixel_spacing = np.array((0.25e-3, 0.25e-3, 0.25e-3)) # resolution of seeding mesh, mm
seeding_slice_bb = Quantity([5, 5.1, 5], "cm")
seeding_particle_density = 5 # per 1/mm*3
# Define initial filling
initial_fill_normal = np.array([1., 0., 0.])
initial_fill_position = Quantity([-0.0005, 0., 0.], "m")
initial_fill_thickness = Quantity(50, "mm")
initial_fill_bounding_box = Quantity([25, 25, 25], "cm")
lookup_mesh_spacing = np.array((0.5e-3, 0.5e-3, 0.5e-3)) # resolution of lookup mesh, mm
Create 3D rendering of Re-seeding Volume at the inlet of the meshed domain#
[19]:
pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(off_screen=True, window_size=(800, 400), shape=(1, 2))
plotter.add_mesh(mesh, scalars=np.linalg.norm(mesh["velocity"], axis=1),
scalar_bar_args=dict(title=r"$\| \vec{v} \|$"), cmap="twilight", opacity=0.9)
box = local_functions.transformed_box(np.eye(3,3), slice_normal=seeding_slice_normal, readout_direction=np.array([0, 0, 1]),
slice_position=seeding_slice_position, extend=seeding_slice_bb)
box2 = local_functions.transformed_box(np.eye(3,3), slice_normal=initial_fill_normal, readout_direction=np.array([0, 0, 1]),
slice_position=initial_fill_position, extend=initial_fill_bounding_box)
box3 = local_functions.transformed_box(np.eye(3,3), slice_normal=slice_normal, readout_direction=readout_direction,
slice_position=slice_position, extend=Quantity([*fov.m_as("cm"), slice_thickness.m_as("cm")], "cm"))
plotter.add_mesh(box, opacity=0.9, color="orange")
# plotter.add_mesh(box2, opacity=0.6, color="green")
plotter.subplot(0, 1)
plotter.add_mesh(mesh, scalars=np.linalg.norm(mesh["velocity"], axis=1),
scalar_bar_args=dict(title=r"$\| \vec{v} \|$"), cmap="twilight", opacity=0.9)
plotter.add_mesh(box3, opacity=0.6, color="gray")
local_functions.add_custom_axes(plotter)
plotter.screenshot("reseeding_box.png");
display(Image("reseeding_box.png"))
Instantiate RefillingFlowDataset#
[20]:
import cmrsim.utils.particle_properties as part_factory
# setup initial fill dict
initial_fill_dict = {'slice_normal': initial_fill_normal,
'slice_position': initial_fill_position.m_as('m'),
'slice_thickness': initial_fill_thickness.m_as('m'),
'slice_bounding_box': initial_fill_bounding_box.m_as('m')}
# Setup flow dataset with seeding mesh information
particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
"T1": part_factory.uniform(default_value=300.),
"T2": part_factory.uniform(default_value=100.),
"M0": part_factory.uniform(default_value=1),
"omega_t2s": part_factory.t2star_lorentzian(T2p=50, prctile_cutoff=0.002)
}
dataset = cmrsim.datasets.RefillingFlowDataset(mesh,
particle_creation_callables=particle_creation_fncs,
slice_position=seeding_slice_position.m_as("m"),
slice_normal=seeding_slice_normal,
slice_thickness=seeding_slice_thickness.m_as("m"),
slice_bounding_box=seeding_slice_bb.m_as("m"),
lookup_map_spacing=lookup_mesh_spacing,
seeding_vol_spacing=seeding_pixel_spacing,
field_list=[("velocity", 3), ("chol", 6), ("Frequency", 1), ("Reflectionmap", 1)])
[21]:
# Setup turbulent particle trajectory module
lookup_table_3d, map_dimensions = dataset.get_lookup_table()
trajectory_module = cmrsim.trajectory.TurbulentTrajectory(velocity_field=lookup_table_3d.astype(np.float32),
map_dimensions=map_dimensions.astype(np.float32),
additional_dimension_mapping=dataset.field_list[3:],
device="GPU:0")
Simulation#
Setup modules#
Define submodules to incorporate \(T_2^*\) effects
[22]:
T2s = cmrsim.bloch.submodules.T2starDistributionModel()
submodules = [T2s, ]
Grid sequence objects (per TR)
[23]:
time_list = []
rf_list = []
grad_list = []
adc_list = []
for sequence_list in tqdm(multi_sequence_list, desc="Grid sequences per VENC..."):
time, rf_grid, grad_grid, adc_on_grid = [np.stack(v) for v in cmrseq.utils.grid_sequence_list(sequence_list, force_uniform_grid=False)]
time_list.append(time)
rf_list.append(rf_grid)
grad_list.append(grad_grid)
adc_list.append(adc_on_grid)
time_ref = time_list[0][0]
adc_all = np.stack(adc_list, axis=1)
grad_all = np.stack(grad_list, axis=1)
rf_all = np.stack(rf_list, axis=1)
print("Bloch definition arrays: \n", time_ref.shape, "\n", adc_all.shape, "\n", grad_all.shape, "\n", rf_all.shape)
Instantiate a Bloch-operator per TR (parallel execution of VENCS)
[24]:
bloch_operator_list = []
for tr in range(grad_all.shape[0]):
bloch_operator_list.append(cmrsim.bloch.GeneralBlochOperator(name="rf", gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
time_grid=time_ref, gradient_waveforms=grad_all[tr],
rf_waveforms=rf_all[tr], adc_events=adc_all[tr], device="GPU:0",
submodules=submodules))
Instantiate Simulator#
[25]:
from cmrsim.simulation_templates.flow import FlowSimulationMultiVenc
simulator = FlowSimulationMultiVenc(dataset, prep_modules=[], readout_modules=bloch_operator_list, trajectory_module=trajectory_module)
Run#
[26]:
r, m = simulator(particle_density=seeding_particle_density, dropping_distance=0.2,
averages=1, reseed_threshold=1, use_inital_filling=True,
initial_filling_dict=initial_fill_dict, return_final_magnetization=True)
[27]:
"Total number of particles at the end of simulation:", r.shape
[27]:
('Total number of particles at the end of simulation:', (1537101, 3))
Reconstruction#
[29]:
signal_list = []
for bloch_mod in bloch_operator_list:
signal_list.append(tf.stack(bloch_mod.time_signal_acc, axis=0).numpy())
signals = np.stack(signal_list,axis=1)
image_list = []
kspace_list = []
for ind in range(np.shape(signals)[0]):
time_signal = np.squeeze(signals[ind, :, :])
n = (np.random.normal(0, 1, time_signal.shape)) * 0
n = n*np.exp(1j*np.random.normal(0, 1, time_signal.shape))
time_signal = time_signal + n
center_index = np.array(matrix_size) // 2 + np.mod(matrix_size, 2)
centered_projection = tf.signal.fft(time_signal).numpy()
centered_k_space = tf.signal.fft(centered_projection).numpy()
image_list.append(tf.signal.fftshift(tf.signal.ifft2d(tf.roll(tf.signal.ifftshift(centered_k_space, axes=(0, 1)), -1, axis=1)), axes=(0,1)).numpy())
kspace_list.append(centered_k_space)
[30]:
kspace_stacked = np.stack(kspace_list)
image_stacked = np.stack(image_list)
np.save("turb_images.npy", image_stacked)
np.save("turb_kspace.npy", kspace_stacked)
image_stacked.shape, kspace_stacked.shape
[30]:
((7, 71, 111), (7, 71, 111))
Infer velocity and turbulence#
Get mask from reference amplitude:
[31]:
from skimage.morphology import opening, closing, remove_small_holes, disk
kspace_stacked = np.load("turb_kspace.npy")
image_stacked = np.load("turb_images.npy")
reference_image = image_stacked[0]
mask = np.abs(reference_image) > 550
mask = closing(remove_small_holes(opening(mask, disk(2))))
mask.shape, reference_image.shape
[31]:
((71, 111), (71, 111))
Calculate velocity projections from phase differences
[32]:
phase_difference = np.mod(np.angle(image_stacked) - np.angle(reference_image) + np.pi, 2 * np.pi) - np.pi
masked_phase_difference = np.ma.masked_less(mask[np.newaxis], 0.5) * phase_difference
Turbulent kinetic energy
[33]:
turb_estimation = np.log(np.abs(reference_image) / np.abs(image_stacked))
masked_trub_estimation = mask[np.newaxis] * turb_estimation
Plot results
[37]:
plt.rcParams["font.family"] = "serif"
kspace_stacked = np.load("turb_kspace.npy")
image_stacked = np.load("turb_images.npy")
concatenated_mask = np.ma.masked_less(np.concatenate([mask ] * 7, axis=1), 0.5)
concatenated_mask[concatenated_mask < 1] = np.nan
# Show images
plt.close("all")
fig, axes = plt.subplots(4, 1, figsize=(13, 7.), constrained_layout=True)
# axes[0].imshow(np.log(np.abs(np.concatenate(kspace_stacked, 1))), cmap="gray")
axes[0].imshow(np.abs(np.concatenate(image_stacked, axis=1)), cmap="gray", vmin=0,
vmax=np.percentile(np.abs(np.squeeze(image_stacked[0])), 99.9))
axes[0].text(42, 10, "z", color="red", fontsize=12)
axes[0].arrow(x=10, y=14, dx=30, dy=0, color="red", width=1)
axes[0].text(12, 40, "y", color="green", fontsize=12)
axes[0].arrow(x=10, y=14, dx=0, dy=30, color="green", width=1)
phase_artist = axes[1].imshow(np.angle(np.concatenate(image_stacked, axis=1)), cmap="RdBu", vmin=-np.pi, vmax=np.pi)
fig.colorbar(phase_artist, location="bottom", orientation="horizontal", shrink=0.7, pad=0.04, fraction=0.8, aspect=100, label="rad")
axes[2].imshow(np.abs(np.concatenate(image_stacked, axis=1)) * concatenated_mask, cmap="gray", vmin=0, vmax=np.percentile(np.abs(np.squeeze(image_stacked[0])),99.9))
phase_arts = axes[2].imshow(np.concatenate(phase_difference, axis=1) * concatenated_mask, cmap="RdBu", vmin=-np.pi, vmax=np.pi, alpha=0.8)
fig.colorbar(phase_arts, location="bottom", orientation="horizontal", shrink=0.7, pad=0.04, fraction=0.8, aspect=100, label="rad")
axes[3].imshow(np.abs(np.concatenate(image_stacked, axis=1)) * concatenated_mask, cmap="gray", vmin=0, vmax=np.percentile(np.abs(np.squeeze(image_stacked[0])),99.9))
turb_arts = axes[3].imshow(np.concatenate(turb_estimation, axis=1) * concatenated_mask, cmap="gist_heat", vmin=0, vmax=1, alpha=0.8)
fig.colorbar(turb_arts, location="bottom", orientation="horizontal", shrink=0.7, pad=0.04, fraction=0.8, aspect=100, label="a.u.")
# axes[0].set_title((" " * 22).join([f"VENC #{i}" for i in range(image_stacked.shape[0])]))
title_string = (" " * 20).join(["No VENC"] + [f"{v.m_as('cm/s'):0.0f} cm/s" for v in venc_list[1:]])
title_string += "\n"
title_string += (" " * 20).join([" ------ "] + [f"{v}" for v in all_venc_directions[1:]])
axes[0].set_title(title_string)
[_.set_ylabel(t, fontsize=14) for _, t in zip(axes, ["Magnitude", "$\phi$", "$\phi - \phi_{Ref}$", "TKE"])]
[_.set_xticks([]) for _ in axes.flatten()];
[_.set_yticks([]) for _ in axes.flatten()];
[axes[0].text(i/7, 1.325, "\u2193"+f"{t})", horizontalalignment='left', verticalalignment='top',
transform=axes[0].transAxes, fontsize=16, fontweight="bold", color="k")
for i, t in enumerate("abcdefg")]
# fig.savefig("Figure7.tiff", dpi=500)abs
fig.savefig("bloch_flow_venc_sGRE_turbUBend.png", dpi=200)
display(Image("bloch_flow_venc_sGRE_turbUBend.png"))
Appendix - Rendering#
Illustrate sequences#
[ ]:
import imageio
plt.ioff()
for tr_idx in tqdm(range(len(multi_sequence_list[0])), leave=False):
plt.close("all");
figure, axes = plt.subplot_mosaic("AABBCCDD;.EEFFGG.", figsize=(12, 5), sharex=True, sharey=True)
axes = [axes[k] for k in "ABCDEFG"]
for venc_idx in range(len(multi_sequence_list)):
cmrseq.plotting.plot_sequence(multi_sequence_list[venc_idx][tr_idx], axes=axes[venc_idx], add_legend=False)
[ax.set_title(f"VENC #{i}") for i, ax in enumerate(axes)]
[(ax.set_ylabel(""), ax.set_yticklabels([])) for ax in axes[1:]]
[(ax.set_ylabel(""), ax.set_yticklabels([])) for ax in figure.axes[8:]]
figure.tight_layout()
figure.savefig(f"venc_seq{tr_idx:02}.png")
plt.close("all")
figure, axis = plt.subplots(1, 1, figsize=(5, 5))
cmrseq.plotting.plot_kspace_2d(multi_sequence_list[0][-1], plot_raster_trajectory=False, ax=axis, k_axes=(1, 2), map_sampling_times="global", add_colorbar=True)
cmrseq.plotting.plot_kspace_2d(multi_sequence_list[0][0], plot_raster_trajectory=False, ax=axis, k_axes=(1, 2), map_sampling_times="global", add_colorbar=False)
xlims = axis.get_xlim()
ylims = axis.get_ylim()
figure.tight_layout()
axis.clear()
for idx, seq in tqdm(enumerate(multi_sequence_list[0]), total=71, leave=False):
cmrseq.plotting.plot_kspace_2d(seq, plot_raster_trajectory=False, ax=axis, k_axes=(1, 2), map_sampling_times="global")
axis.set_ylim(ylims); axis.set_xlim(xlims)
figure.savefig(f"kspace{idx:02}")
_mr_seq_def_pngs = [imageio.v3.imread(f"venc_seq{idx:02}.png") for idx in tqdm(range(len(multi_sequence_list[0])), leave=False)]
_mr_kspace_pngs = [imageio.v3.imread(f"kspace{idx:02}.png") for idx in tqdm(range(len(multi_sequence_list[0])), leave=False)]
comb_pngs = [np.concatenate([i1, i2], axis=1) for i1, i2 in zip(_mr_seq_def_pngs, _mr_kspace_pngs)]
imageio.mimsave("mr_sequence_animation.gif", comb_pngs, duration=100, loop=0)
[os.remove(f"venc_seq{idx:02}.png") for idx in tqdm(range(len(multi_sequence_list[0])), leave=False)]
[os.remove(f"kspace{idx:02}.png") for idx in tqdm(range(len(multi_sequence_list[0])), leave=False)]