Flow BSSFP#
Import modules and set TF-GPU configuration#
[1]:
!pip install -U cmrseq==0.23 --index-url https://gitlab.ethz.ch/api/v4/projects/30300/packages/pypi/simple
Defaulting to user installation because normal site-packages is not writeable
Looking in indexes: https://gitlab.ethz.ch/api/v4/projects/30300/packages/pypi/simple
Requirement already satisfied: cmrseq==0.23 in /home/jweine/.local/lib/python3.11/site-packages (0.23)
Requirement already satisfied: numpy in /usr/local/lib/python3.11/dist-packages (from cmrseq==0.23) (1.26.2)
Requirement already satisfied: pint in /usr/local/lib/python3.11/dist-packages (from cmrseq==0.23) (0.23)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.11/dist-packages (from cmrseq==0.23) (3.8.3)
Requirement already satisfied: tqdm in /usr/local/lib/python3.11/dist-packages (from cmrseq==0.23) (4.66.2)
Requirement already satisfied: scipy in /usr/local/lib/python3.11/dist-packages (from cmrseq==0.23) (1.12.0)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.11/dist-packages (from matplotlib->cmrseq==0.23) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.11/dist-packages (from matplotlib->cmrseq==0.23) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.11/dist-packages (from matplotlib->cmrseq==0.23) (4.50.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /usr/local/lib/python3.11/dist-packages (from matplotlib->cmrseq==0.23) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.11/dist-packages (from matplotlib->cmrseq==0.23) (23.2)
Requirement already satisfied: pillow>=8 in /usr/local/lib/python3.11/dist-packages (from matplotlib->cmrseq==0.23) (10.2.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/lib/python3/dist-packages (from matplotlib->cmrseq==0.23) (2.4.7)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.11/dist-packages (from matplotlib->cmrseq==0.23) (2.9.0.post0)
Requirement already satisfied: typing-extensions in /usr/local/lib/python3.11/dist-packages (from pint->cmrseq==0.23) (4.8.0)
Requirement already satisfied: six>=1.5 in /usr/lib/python3/dist-packages (from python-dateutil>=2.7->matplotlib->cmrseq==0.23) (1.16.0)
[notice] A new release of pip is available: 23.3.1 -> 24.0
[notice] To update, run: python3 -m pip install --upgrade pip
[2]:
import sys
sys.path.append("../")
sys.path.insert(0, "../../")
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)
[3]:
from IPython.display import clear_output
import pyvista
from pint import Quantity
import numpy as np
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
%matplotlib inline
import cmrsim
from cmrsim.utils import particle_properties as part_factory
import cmrseq
import local_functions
clear_output()
Define Sequence#
[4]:
# Define MR-system specifications
system_specs = cmrseq.SystemSpec(max_grad=Quantity(50, "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"))
# Define Fourier-Encoding parameters
fov = Quantity([194, 122], "mm")
matrix_size = (97, 61)
res = fov / matrix_size
print("Resolution:", res)
n_dummy = 25
center_index = np.array(matrix_size) // 2 + np.mod(matrix_size, 2)
print(center_index)
# Define imaging slice parameters
slice_thickness = Quantity(6, "mm")
slice_position = Quantity([0., 0., 0.], "m")
adc_duration = Quantity(4., "ms")
readout_direction = np.array([1., 0., 0.]) # readout along scanner Z, (FH)
slice_normal = np.array([0., 1., 0.])
slice_position_offset = Quantity(0.5, "cm")
# Define RF-exication parameters (coupled with seeding slice atm)
flip_angle = Quantity(np.pi / 4, "rad")
pulse_duration = Quantity(1., "ms")
RF_slice_normal = np.array([0., 1., 0.]) # Slice along scanner Y, (AP)
RF_slice_position = Quantity([0., 0., 0.], "m")
FA = np.pi/4
# Define Reseeding slice parameters
seeding_slice_position = Quantity([0.03, 0.005, -0.025], "m")
seeding_slice_normal = np.array([0.3, 0., 1.])
seeding_slice_thickness = Quantity(35, "mm")
seeding_pixel_spacing = np.array((0.5e-3, 0.5e-3, 0.5e-3)) # resolution of seeding mesh, mm
seeding_slice_bb = Quantity([6, 6, 4], "cm")
seeding_particle_density = 10 # per 1/mm*3
# Define initial filling
initial_fill_normal = np.array([0., 1., 0.])
initial_fill_position = Quantity([0., 0.005, 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
[49 31]
[5]:
seq_list = cmrseq.seqdefs.sequences.balanced_ssfp(system_specs,
matrix_size,
inplane_resolution=res,
slice_thickness=slice_thickness,
adc_duration=adc_duration,
flip_angle=flip_angle,
pulse_duration=pulse_duration,
slice_position_offset=slice_position_offset,
repetition_time=Quantity(0, "ms"), # defaults to shortest
dummy_shots=n_dummy)
rf_rotation_matrix = cmrseq.utils.get_rotation_matrix(RF_slice_normal, np.array([0., 0., 1.]))
ro_rotation_matrix = cmrseq.utils.get_rotation_matrix(slice_normal, readout_direction)
for rep_idx, seq in enumerate(seq_list):
[b.rotate_gradients(rf_rotation_matrix) for b in seq.get_block(partial_string_match="slice")]
[b.rotate_gradients(rf_rotation_matrix) for b in seq.get_block(partial_string_match=["slice", "rf", "adc"], invert_pattern=True)]
_, k_adc, t_adc = seq_list[n_dummy + center_index[1]].calculate_kspace()
plt.close("all")
f, (a1, a2) = plt.subplots(1, 2, figsize=(16, 6), gridspec_kw={"width_ratios": (2, 1)})
cmrseq.plotting.plot_sequence(seq_list[0], axes=a1, format_axes=True)
for s in seq_list[n_dummy:]:
cmrseq.plotting.plot_sequence(s, axes=[f.axes[-1], a1, a1, a1], format_axes=False)
for s in seq_list[n_dummy:]:
cmrseq.plotting.plot_kspace_2d(s, ax=a2, k_axes=[0, 2])
CMRSeq Warning : bSSFP Sequence: Repetition time too short to be feasible, set TR to 6.2 millisecond
Set up phantom#
[6]:
# Load Aorta Example.
mesh = pyvista.read("../example_resources/pv_computed_averaged_P14.vtk")
# 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["U1"], "m/s").m_as("m/ms") * 0.4
# To match the usual geometry (Head <-> Foot along z-axis) flip the axes of the mesh
mesh.points = mesh.points[:, [2, 1, 0]]
mesh.points[:, 1:] *= -1
mesh["velocity"] = mesh["velocity"][:, [2, 1, 0]]
mesh["velocity"][:, 1:] *= -1
global_mesh_offset = [0.02, 0.01, -0.02]
mesh.translate(global_mesh_offset, inplace=True)
# Rotate mesh and flow by 60 degrees
mesh.rotate_z(60, (0., 0., 0.), inplace=True)
vx = np.cos(60*np.pi/180)*mesh["velocity"][:, 0] - np.sin(60 * np.pi / 180) * mesh["velocity"][:, 1]
mesh["velocity"][:, 1] = np.sin(60 * np.pi / 180) * mesh["velocity"][:, 0] + np.cos(60 * np.pi / 180) * mesh["velocity"][:, 1]
mesh["velocity"][:, 0] = vx
mesh = mesh.cell_data_to_point_data()
# clip to restricted length to reduce mesh-size
mesh.clip(normal='z', origin=[0., 0., 0.09], inplace=True)
mesh.clip(normal='-z', origin=[0., 0., -0.09], inplace=True)
# setup initial fill dictionary
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')}
# Add artificial channel for off_resonance
TR = seq_list[-1].duration
mesh["off_res"] = (mesh.points[:, 0] < 0) / 2 / TR.m_as("ms") / system_specs.gamma.m_as("MHz/T")
# Setup flow dataset with seeding mesh information
# Here T1, T2 and T2' are set
particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
"T1": part_factory.uniform(default_value=1000.),
"T2": part_factory.uniform(default_value=200.),
"omega_t2s": part_factory.t2star_lorentzian(T2p=40, prctile_cutoff=0.002),
"M0": part_factory.uniform(default_value=1/seeding_particle_density)}
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"),
seeding_vol_spacing=seeding_pixel_spacing,
lookup_map_spacing=lookup_mesh_spacing,
field_list=[("velocity", 3), ("off_res", 1)])
mesh
Updated Mesh
Updated Slice Position
[6]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Inspect slice configuration
[7]:
pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(notebook=True, window_size=(800, 500))
# Plot input mesh
plotter.add_mesh(mesh, color="gray", opacity=0.6)
local_functions.add_custom_axes(plotter)
# Plot slice selective RF
_scale = 0.1
rf_mesh = pyvista.Box(bounds=(-_scale, _scale, -_scale, _scale, -slice_thickness.m_as("m") / 2, slice_thickness.m_as("m") / 2))
rf_mesh.points = cmrseq.utils.mps_to_xyz(rf_mesh.points[np.newaxis], slice_normal=RF_slice_normal, readout_direction=readout_direction)[0]
rf_mesh.translate(RF_slice_normal * slice_position_offset.m_as("m"), inplace=True)
plotter.add_mesh(rf_mesh, opacity=0.1, show_edges=True, color="red")
# # Plot projection plane of the cartesian image slice with M-S directions
image_slice_mesh = pyvista.Box(bounds=(-fov[0].m_as("m")/2, fov[0].m_as("m")/2, -fov[1].m_as("m")/2, fov[1].m_as("m")/2, 0, 0))
image_slice_mesh.points = cmrseq.utils.mps_to_xyz(image_slice_mesh.points[np.newaxis], slice_normal=slice_normal, readout_direction=readout_direction)[0]
plotter.add_mesh(image_slice_mesh, opacity=1, show_edges=True)
# Add Initial fill volume
plotter.add_mesh(dataset.gridded_seeding_volume, color="orange", opacity=1)
plotter.show(jupyter_backend="static")
Simulation#
Instantiate Bloch operators#
[8]:
from cmrsim.simulation_templates.flow import FlowSimulation
# 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 = [T2s, offres]
time_cat, rf_grid_cat, grad_grid_cat, _ = [np.stack(v) for v in cmrseq.utils.grid_sequence_list([seq_list[0]], force_uniform_grid=False)]
time_dummy, rf_grid_dummy, grad_grid_dummy, _ = [np.stack(v) for v in cmrseq.utils.grid_sequence_list(seq_list[1:(n_dummy+1)], force_uniform_grid=False)]
time, rf_grid, grad_grid, adc_on_grid = [np.stack(v) for v in cmrseq.utils.grid_sequence_list(seq_list[(n_dummy+1):], force_uniform_grid=False)]
print(time.shape, rf_grid.shape, grad_grid.shape, adc_on_grid.shape)
bloch_catalyst = cmrsim.bloch.GeneralBlochOperator(name="rf", gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
time_grid=time_cat[0], gradient_waveforms=grad_grid_cat,
rf_waveforms=rf_grid_cat, device="GPU:0",
submodules=submodules)
bloch_dummy = cmrsim.bloch.GeneralBlochOperator(name="rf", gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
time_grid=time_dummy[0], gradient_waveforms=grad_grid_dummy,
rf_waveforms=rf_grid_dummy, device="GPU:0",
submodules=submodules)
bloch_readout = cmrsim.bloch.GeneralBlochOperator(name="rf", gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
time_grid=time[0], gradient_waveforms=grad_grid,
rf_waveforms=rf_grid, adc_events=adc_on_grid,
device="GPU:0", submodules=submodules)
(61, 717) (61, 717) (61, 717, 3) (61, 717, 2)
Optional: Check Slice configuration#
Run Simulation#
[9]:
# Instantiate simulator:
simulator = FlowSimulation(flow_dataset=dataset, prep_modules=[bloch_catalyst, bloch_dummy], readout_module=bloch_readout)
n_steps_total = grad_grid.shape[1]
print(f"Total time-steps per TR: {n_steps_total}")
r, m = simulator(particle_density=seeding_particle_density,
dropping_distance=0.15, averages=1,
reseed_threshold=0.3, use_inital_filling=True,
initial_filling_dict=initial_fill_dict,
return_final_magnetization=True)
Total time-steps per TR: 717
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1712949523.638421 3187983 device_compiler.h:186] Compiled cluster using XLA! This line is logged at most once for the lifetime of the process.
Reconstruct#
First obtain the accumulated time-domain signal from the encoding module. This signal is discretized on the simulation raster. The encoding module yields the centered spatial projection corresponding to the time-domain signals fourier transform.
To obtain the projection with “removed oversampling” the projections are linearly interpolated to the positions corresponding to the instantaneous encoding times t_adc as defined above. By calling the encoding modules get_kspace() these cropped projections are again fourier transformed, yielding the 2D cartesian k-space.
[10]:
time_signal = tf.stack(bloch_readout.time_signal_acc, axis=0).numpy().reshape(matrix_size[::-1])
# Add noise based on max (DC) time signal at SNR of >100. This does not affect the magnitude image,
# but the phase image will not have structures in regions of zero signal
n = (np.random.normal(0, 1, time_signal.shape)) * np.max(np.abs(time_signal)) / 1000
n = n*np.exp(1j*np.random.normal(0, 1, time_signal.shape))
time_signal = time_signal + n
image = tf.signal.fftshift(tf.signal.ifft2d(tf.signal.ifftshift(time_signal, axes=(0, 1))), axes=(0, 1)).numpy()
time_per_seg = t_adc - t_adc[center_index[0]]
freq_per_seg = np.fft.fftfreq(n=time_per_seg.shape[0], d=time_per_seg[1]-time_per_seg[0])
print(time_signal.shape, image.shape)
(61, 97) (61, 97)
[11]:
# Show images
fig, axes = plt.subplots(1, 3, figsize=(12, 2.5))
axes[0].imshow(np.log(np.abs(np.squeeze(time_signal))), cmap="gray")
axes[1].imshow(np.abs(np.squeeze(image)), cmap="gray")
axes[2].imshow(np.angle(np.squeeze(image)), cmap="gray")
[_.axis("off") for _ in axes]
[_.set_title(t) for _, t in zip(axes, ["k-space", "Magnitude", "Phase"])]
fig.tight_layout()
Visualize final particle density#
This is a good way to ensure that there are no reseeding or particle trajectory errors, such as density voids or pileups
[12]:
fig, axes = plt.subplots(1, 1, figsize=(9, 5))
plot_axes = [0, 2] # Change these to change the axes of the projection
A, xe, ye, _ = plt.hist2d(r[:, plot_axes[0]], r[:, plot_axes[1]], 600)
fig.clear()
asp = (np.max(xe) - np.min(xe)) / (np.max(ye) - np.min(ye)) * 9
fig.set_size_inches(9, asp)
plt.imshow(A, aspect='auto', vmin=0, vmax=50)
plt.axis('off')
plt.show()
