Turbulent flow with velocity encoded spoiled GRE#
Imports#
[ ]:
!pip install cmrseq
[1]:
# Load TF and check for GPUs
import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
import tensorflow as tf
gpus = tf.config.list_physical_devices('GPU')
for gpu in gpus:
tf.config.experimental.set_memory_growth(gpu, True)
logical_gpus = tf.config.list_logical_devices('GPU')
print(len(gpus), "Physical GPUs,", len(logical_gpus), "Logical GPUs")
# 3rd Party dependencies
from pint import Quantity
import pyvista
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
%matplotlib inline
from IPython.display import display, Image, clear_output
import imageio
from tqdm.notebook import tqdm
# Project library cmrsim
import sys
sys.path.append("../../")
sys.path.append("../")
sys.path.insert(0, "/scratch/jweine/cmrseq")
import cmrseq
import cmrsim
import cmrsim.utils.particle_properties as part_factory
import local_functions
1 Physical GPUs, 1 Logical GPUs
Sequence Definition#
[3]:
# 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") # 1
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") # 0.5
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.
[4]:
multi_sequence_list: list[list[cmrseq.Sequence]] = cmrseq.contrib.quantitative_flow.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)
/scratch/jweine/cmrseq/cmrseq/parametric_definitions/velocity.py:134: UserWarning: Bipolar gradient duration set too short
warn("Bipolar gradient duration set too short")
/opt/conda/lib/python3.10/site-packages/pint/facets/plain/quantity.py:986: RuntimeWarning: invalid value encountered in double_scalars
magnitude = magnitude_op(new_self._magnitude, other._magnitude)
/scratch/jweine/cmrseq/cmrseq/contrib/quantitative_flow.py:107: UserWarning: Repetition time too short to be feasible, set TR to 8.68 millisecond
warn(f"Repetition time too short to be feasible, set TR to {minimal_tr}")
/scratch/jweine/cmrseq/cmrseq/contrib/quantitative_flow.py:115: UserWarning: Echo time too short to be feasible, set TE to 5.0600000000000005 millisecond
warn(f"Echo time too short to be feasible, set TE to {minimal_te}")
<IPython.core.display.Image object>
Phantom preparation#
Load U-bend mesh#
The
[5]:
mesh = pyvista.read("../example_resources/sUbend_219_9_rotated.vtu")
mesh
[5]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
[8]:
_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=(600, 300), shape=(1, 2))
plotter.subplot(0, 0)
plotter.add_mesh(_tmpmesh, scalars="abs_vel", scalar_bar_args=dict(title="$\| \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#
[10]:
# 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
[190]:
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="$\| \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="$\| \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#
[12]:
# 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)])
[13]:
# 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
[14]:
T2s = cmrsim.bloch.submodules.T2starDistributionModel()
submodules = [T2s, ]
Grid sequence objects (per TR)
[15]:
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)
[16]:
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#
[17]:
from cmrsim.simulation_templates.flow import FlowSimulationMultiVenc
simulator = FlowSimulationMultiVenc(dataset, prep_modules=[], readout_modules=bloch_operator_list, trajectory_module=trajectory_module)
Run#
[18]:
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)
[25]:
"Total number of particles at the end of simulation:", r.shape
[25]:
('Total number of particles at the end of simulation:', (1553888, 3))
Reconstruction#
[20]:
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)
[37]:
[37]:
((7, 71, 111), (7, 71, 111))
Infer velocity and turbulence#
Get mask from reference amplitude:
[124]:
reference_image = image_list[0]
mask = np.abs(reference_image) > 450
mask.shape, reference_image.shape
[124]:
((71, 111), (71, 111))
Calculate velocity projections from phase differences
[125]:
?
[126]:
turb_estimation = np.log(np.abs(reference_image) / np.abs(image_stacked))
masked_trub_estimation = mask[np.newaxis] * turb_estimation
Plot results
[189]:
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)
# Show images
plt.close("all")
fig, axes = plt.subplots(5, 1, figsize=(16, 9), constrained_layout=True)
axes[0].imshow(np.log(np.abs(np.concatenate(kspace_stacked, 1))), cmap="gray")
axes[1].imshow(np.abs(np.concatenate(image_stacked, axis=1)), cmap="gray", vmin=0, vmax=np.percentile(np.abs(np.squeeze(image)),99.9))
axes[1].text(42, 10, "z", color="red", fontsize=12)
axes[1].arrow(x=10, y=14, dx=30, dy=0, color="red", width=1)
axes[1].text(12, 40, "y", color="green", fontsize=12)
axes[1].arrow(x=10, y=14, dx=0, dy=30, color="green", width=1)
phase_artist = axes[2].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[3].imshow(np.abs(np.concatenate(image_stacked, axis=1)), cmap="gray", vmin=0, vmax=np.percentile(np.abs(np.squeeze(image)),99.9))
phase_arts = axes[3].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[4].imshow(np.abs(np.concatenate(image_stacked, axis=1)), cmap="gray", vmin=0, vmax=np.percentile(np.abs(np.squeeze(image)),99.9))
turb_arts = axes[4].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 = (" " * 25).join([f"{v.m_as('cm/s'):0.0f} cm/s" for v in venc_list])
title_string += "\n"
title_string += (" " * 25).join([f"{v}" for v in all_venc_directions])
axes[0].set_title(title_string)
[_.set_ylabel(t) for _, t in zip(axes, ["k-space", "Magnitude", "$\phi$", "$\phi - \phi_{Ref}$", "??"])]
[_.set_xticks([]) for _ in axes.flatten()];
[_.set_yticks([]) for _ in axes.flatten()];
fig.savefig("simulation_results.svg")
Appendix - Rendering#
Illustrate sequences#
[ ]:
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)]
Backup#
[4]:
# Load UBend Example.
mesh = pyvista.read("../example_resources/sUbend_219_9.vtu")
# The velocity [m/s] is stored in the mesh as U1. In cmrsim all definitions are done in [ms] therefore
# We add the velocity field as scaled data to the mesh
mesh["velocity"] = Quantity(mesh["U"], "m/s").m_as("m/ms")
mesh.points = mesh.points[:,[2,1,0]]
mesh["velocity"] = mesh["velocity"][:,[2,1,0]]
mesh["RST"] = mesh["RST"][:,[5,4,2,3,1,0]]
mesh["Reflectionmap"] = mesh["U"][:,0]*0 + 1 # in ms
tfreq = np.sqrt(np.sum(mesh["RST"][:,[0,3,5]],1));
tfreq = tfreq / np.percentile(tfreq,90) * 1000
mesh["Frequency"] = tfreq # in Hz
mesh["RST"] = mesh["RST"] # m^2/s^2
# RST form xx xy xz yy yz zz
a11=np.sqrt(np.clip(mesh["RST"][:,0],0,None) + np.min(np.abs(mesh["RST"][:,0])))
a21=mesh["RST"][:,1]/a11
a22=np.sqrt(np.clip(mesh["RST"][:,3]-a21**2,0,None) + np.min(np.abs(mesh["RST"][:,3]-a21**2)))
a31=mesh["RST"][:,2]/a11
a32=(mesh["RST"][:,4]-a21*a31)/a22
a33=np.sqrt(np.clip(mesh["RST"][:,5]-a31**2-a32**2,0,None))
# Chol form
mesh["chol"] = np.stack([a11,a21,a22,a31,a32,a33],axis=1)/1000/5 # m/ms
global_mesh_offset = [0, 0.04, -0.05]
mesh.translate(global_mesh_offset, inplace=True)
mesh = mesh.cell_data_to_point_data()
# clip to restricted length to reduce mesh-size
mesh.clip(normal='-z', origin=[0., 0., -0.1], inplace=True)
mesh.save("../example_resources/sUbend_219_9_rotated.vtu")
[4]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|