Multi-Velocity encoded spoiled GRE laminar flow (parallel execution)#
Imports#
[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 widget
# Project library cmrsim
import sys
sys.path.append("../../")
sys.path.append("../../../cmrseq")
import cmrseq
import cmrsim
import cmrsim.utils.particle_properties as part_factory
1 Physical GPUs, 1 Logical GPUs
Sequence Definition#
[2]:
# 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
TR = Quantity(5, 'ms')
TE = Quantity(0, 'ms')
# Velocity encoding settings:
venc_list = [Quantity(0.,'m/s'), Quantity(5.5,'m/s'), Quantity(2.,'m/s'), Quantity(1.,'m/s')]
venc_direction_list = [np.array([0,0,1]), np.array([0,0,1]), np.array([0,1,0]), np.array([1,0,0])]
venc_duration = Quantity(1, 'ms') # Duration = 0 will result in shortest possible venc
crusher_duration = Quantity(0.3, "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(1., "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")
# Define RF-exication parameters (coupled with seeding slice atm)
flip_angle=Quantity(np.pi/12, "rad")
pulse_duration = Quantity(0.5, "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 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, 2], "cm")
seeding_particle_density = 10 # 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
Resolution: [2.0 2.0] millimeter
[3]:
rotation_matrix = cmrseq.utils.get_rotation_matrix(slice_normal, readout_direction, target_orientation='mps').T
venc_direction_list_MPS = [rotation_matrix@d for d in venc_direction_list]
crusher_area = Quantity(np.pi * 4, "rad") / system_specs.gamma_rad / res[0]
multi_list = cmrseq.contrib.sequences.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=venc_direction_list_MPS,
venc_duration=venc_duration,
crusher_area=crusher_area,
crusher_duration=crusher_duration)
for sequence_list in multi_list:
for seq in sequence_list:
seq.rotate_gradients(cmrseq.utils.get_rotation_matrix(slice_normal, readout_direction, target_orientation='xyz').T)
/scratch/jweine/conda/envs/tf29_pyvista/lib/python3.9/site-packages/cmrseq/parametric_definitions/velocity.py:134: UserWarning: Bipolar gradient duration set too short
warn("Bipolar gradient duration set too short")
/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)
/scratch/jweine/conda/envs/tf29_pyvista/lib/python3.9/site-packages/cmrseq/contrib/sequences.py:115: UserWarning: Echo time too short to be feasible, set TE to 3.1799999999999997 millisecond
warn(f"Echo time too short to be feasible, set TE to {minimal_te}")
[4]:
_, k_adc, t_adc = multi_list[2][0].calculate_kspace()
r = len(multi_list)
plt.close("all")
f, ax = plt.subplots(r, 2, figsize=(10, 3*r), gridspec_kw={'width_ratios': [5, 2]})
for seq,idx in zip(multi_list,range(r)):
cmrseq.plotting.plot_sequence(seq[0], axes=ax[idx,0], add_legend=False)
for s in seq[1:]:
cmrseq.plotting.plot_sequence(s, axes=[ax[idx,0], ax[idx,0], ax[idx,0], ax[idx,0]], format_axes=False, add_legend=False)
for s in seq:
cmrseq.plotting.plot_kspace_2d(s, ax=ax[idx,1], k_axes=[1, 2])
f.tight_layout()
[11]:
from cmrsim.simulation_templates.flow import FlowSimulation
from cmrsim.simulation_templates.flow import FlowSimulationMultiVenc
# Define off-resonsance and T2* modules
offres = cmrsim.bloch.submodules.OffResonance(gamma=system_specs.gamma.m_as("MHz/T"), device="GPU:0")
T2s = cmrsim.bloch.submodules.T2starDistributionModel()
submodules = [offres, T2s]
[12]:
time_list = []
rf_list = []
grad_list = []
adc_list = []
for sequence_list in multi_list:
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)
bloch_list = []
for tr in range(grad_all.shape[0]):
bloch_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))
[13]:
# 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")
# To match the usual geometry (Head <-> Foot along z-axis) flip the axes of the mesh
mesh.points = mesh.points[:,[2,1,0]]
mesh["velocity"] = mesh["velocity"][:,[2,1,0]]
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)
# 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')}
# Detailed artificial off resonance
#TR = seq_list[0].duration
#bandOR = 1/TR.m_as("ms")/ system_specs.gamma.m_as("MHz/T");
mesh["off_res"] = mesh.points[:, 1] * 0
# Add artificial channel for off_resonance
#banding_dist = Quantity(10, "cm")
#banding_shift = Quantity(3, "cm")
#TR = seq_list[0].duration
#mesh["off_res"] = (mesh.points[:, 2] / banding_dist.m_as("m") + banding_shift/banding_dist) / TR.m_as("ms") / system_specs.gamma.m_as("MHz/T")
# 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), ("off_res", 1)])
mesh
Updated Mesh
Updated Slice Position
[13]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
[14]:
simulator = FlowSimulationMultiVenc(dataset, prep_modules=[], readout_modules=bloch_list)
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)
WARNING:tensorflow:5 out of the last 5 calls to <function GeneralBlochOperator._call_core at 0x7fb252505e50> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for more details.
WARNING:tensorflow:6 out of the last 6 calls to <function GeneralBlochOperator._call_core at 0x7fb2520e2ca0> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for more details.
[15]:
signal_list = []
for bloch_mod in bloch_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))*100
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)
[16]:
# Show images
fig, axes = plt.subplots(len(venc_list), 3, figsize=(12, 2.5*len(venc_list)))
for image,k_space,idx in zip(image_list,kspace_list,range(len(venc_list))):
axes[idx,0].imshow(np.log(np.abs(np.squeeze(k_space))), cmap="gray")
axes[idx,1].imshow(np.abs(np.squeeze(image)), cmap="gray", vmin = 0, vmax = np.percentile(np.abs(np.squeeze(image)),99.9))
axes[idx,2].imshow(np.angle(np.squeeze(image)), cmap="RdBu",vmin=-np.pi,vmax = np.pi)
[_.axis("off") for _ in axes.flatten()]
[_.set_title(t) for _, t in zip(axes[0,:], ["k-space", "Magnitude", "Phase"])]
fig.tight_layout()
[17]:
fig, axes = plt.subplots(len(venc_list)-1, 2, figsize=(8, 2.5*(len(venc_list)-1)))
ref_im = image_list[0]
mask = np.abs(ref_im)>0
axes = axes.flatten()
for image,idx in zip(image_list[1:],range(len(venc_list)-1)):
axes[2*idx].imshow(mask*np.mod(np.angle(np.squeeze(image))-np.angle(np.squeeze(ref_im)) + np.pi,2*np.pi)-np.pi, cmap="RdBu", vmin = -np.pi,vmax = np.pi)
axes[2*idx+1].imshow(mask*np.log(np.abs(np.squeeze(ref_im))/np.abs(np.squeeze(image))), cmap="gist_heat",vmin=0,vmax=1)
[_.axis("off") for _ in axes.flatten()]
fig.tight_layout()
[18]:
fig, axes = plt.subplots(len(venc_list), 1, figsize=(5, 4*(len(venc_list)-1)))
ref_im = image_list[0]
mask = np.abs(ref_im)>0
axes = axes.flatten()
vlist = []
for image,idx in zip(image_list[1:],range(len(venc_list)-1)):
venc_cur = venc_list[idx+1].m_as("cm/s")
v = (np.mod(np.angle(np.squeeze(image))-np.angle(np.squeeze(ref_im)) + np.pi,2*np.pi)-np.pi)/np.pi*venc_cur
vlist.append(v)
axes[idx].imshow(v, cmap="RdBu", vmin = -venc_cur,vmax = venc_cur)
axes[-1].imshow(np.sqrt(np.sum(np.asarray(vlist)**2,axis=0)),cmap="viridis",vmin=0,vmax=600)
[_.axis("off") for _ in axes.flatten()]
fig.tight_layout()
[19]:
fig, axes = plt.subplots(1, 1, figsize=(12, 5))
A,xe,ye,_ = plt.hist2d(r[:,1],r[:,2],600)
fig.clear()
asp = (np.max(xe) - np.min(xe)) / (np.max(ye) - np.min(ye)) * 12
fig.set_size_inches(12,asp)
plt.imshow(A,aspect='auto',vmax=50)
plt.show()