Spoiled GRE Cartesian Readout

In this notebook the sequence is composed with cmrseq and the input is a shepp-logan phantom, constructed with the phantominator package link Reconstruction is performed with regular iFFt, as k-space sampling is uniform.

Imports

!pip install cmrseq --index-url https://gitlab.ethz.ch/api/v4/projects/30300/packages/pypi/simple

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)

from copy import deepcopy
from tqdm.notebook import tqdm
from pint import Quantity
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# Project library cmrsim
import cmrsim
import cmrseq
from cmrsim.utils import particle_properties as part_factory

Construct Sequence

In this example a spoiled gradient echo (FLASH / sGRE) is used. To obtain the sequence definition, corresponding function of the cmrseq package is used. The number of dummy shots (n_dummy) defines the amount of time used for reaching steady state, which has a significant impact on the simulation result.

# Define MR-system specifications
system_specs = cmrseq.SystemSpec(max_grad=Quantity(40, "mT/m"),
                                 max_slew=Quantity(200., "mT/m/ms"),
                                 grad_raster_time=Quantity(0.01, "ms"),
                                 rf_raster_time=Quantity(0.01, "ms"))

# Define most important sequence parameters
fov = Quantity([202, 202], "mm")
matrix_size = (51, 51)
res = fov / matrix_size
adc_duration = Quantity(2., "ms")
pulse_duration = Quantity(1.5, "ms")
slice_thickness = Quantity(2, "cm")
flip_angle = Quantity(12, "degree").to("rad")

print("Resolution:", res)

# Call sequence definitions
n_dummy = 5
sequence_list = cmrseq.seqdefs.sequences.flash(system_specs,
                                               slice_thickness=slice_thickness,
                                               flip_angle=flip_angle,
                                               pulse_duration=pulse_duration,
                                               time_bandwidth_product=6,
                                               matrix_size=matrix_size,
                                               inplane_resolution=res,
                                               adc_duration=adc_duration,
                                               echo_time=Quantity(2., "ms"),
                                               repetition_time=Quantity(4.5, "ms"),
                                               dummy_shots=n_dummy)

# Add crusher gradients
crusher_area = Quantity(np.pi * 6, "rad") / system_specs.gamma_rad / res[0]
crusher = cmrseq.bausteine.TrapezoidalGradient.from_area(system_specs, orientation=np.array([1., 0., 0.]),
                                                         area=crusher_area.to("mT/m*ms"), delay=Quantity(4.5, "ms"))
sequence_list = [seq + cmrseq.Sequence([deepcopy(crusher)], system_specs) for seq in sequence_list]

k_grad, k_adc, t_adc = sequence_list[n_dummy+1].calculate_kspace()
t_center = Quantity(t_adc[matrix_size[0]//2], "ms")

combined_sequence_dummyshots = deepcopy(sequence_list[0])
combined_sequence_dummyshots.extend(sequence_list[1:n_dummy])

combined_sequence = deepcopy(sequence_list[n_dummy])
combined_sequence.extend(deepcopy(sequence_list[n_dummy+1:]))

# Plot the sequence
plt.close("all")

fig = plt.figure(constrained_layout=True, figsize=(15, 8))
gs = fig.add_gridspec(2, 2)
f3_ax1 = fig.add_subplot(gs[0, 0])
f3_ax1.set_title('Sequence Per TR')
f3_ax2 = fig.add_subplot(gs[0, 1])
f3_ax2.set_title('k-space-vectors')
f3_ax3 = fig.add_subplot(gs[1, 0:2])
f3_ax3.set_title('Full Sequence')

axes = fig.axes

cmrseq.plotting.plot_sequence(sequence_list[0],  axes[0], legend_position="lower left")
for tr_idx, s in enumerate(sequence_list[1:]):
    cmrseq.plotting.plot_sequence(s,  [fig.axes[-1], axes[0], axes[0], axes[0]], format_axes=False, add_legend=False)
    if tr_idx > n_dummy:
        cmrseq.plotting.plot_kspace_2d(s, ax=axes[1], k_axes=[0, 1])
_ = cmrseq.plotting.plot_sequence(combined_sequence, axes=axes[2], format_axes=True)

Resolution: [3.9607843137254903 3.9607843137254903] millimeter

CMRSeq Warning : FLASH Sequence: Echo time too short to be feasible, set TE to 2.14 millisecond
Extending Sequence: 100%|██████████| 4/4 [00:00<00:00, 340.63it/s]
Extending Sequence: 100%|██████████| 50/50 [00:00<00:00, 317.22it/s]

Set up Phantom

import phantominator
x, y, z = 600, 600, 5
m0, t1, t2 = phantominator.shepp_logan(N=(x, y, z), MR=True, B0=1)
r = np.stack(np.meshgrid(np.linspace(-0.1, 0.1, x), np.linspace(-0.1, 0.1, y),
                         np.linspace(-0.1, 0.1, z)), axis=-1).astype(np.float32)

choose = np.where(m0 > 0)
r, m0, t1, t2 = [_[choose] for _ in (r, m0, t1, t2)]
n_new = r.shape[0]
random_offsets = np.random.uniform(-1, 1, size=(n_new, 3)) * 1e-4
r += random_offsets

properties = {"M0": m0.flatten().astype(np.float32),
              "T1": t1.flatten().astype(np.float32) * 1000,
              "T2": t2.flatten().astype(np.float32) * 1000,
              "magnetization": part_factory.norm_magnetization()(n_new)}
m_init = properties.pop("magnetization")
print(f"Phantom with {n_new} particles")

Phantom with 466820 particles

Grid sequence and configure Bloch Module

The dummy shot simulation and the TRs containing an acquisition bloch are captured by one Bloch module each, where the TRs are indexed in the simulation loop using the rep_idx

t_comb_dummy, rf_comb_dummy, grad_comb_dummy, adc_on_comb_dummy = [np.stack(v) for v in cmrseq.utils.grid_sequence_list([combined_sequence_dummyshots, ])]
t_comb, rf_comb, grad_comb, adc_on_comb = [np.stack(v) for v in cmrseq.utils.grid_sequence_list(sequence_list[n_dummy:])]

print([v.shape for v in [t_comb, rf_comb, grad_comb, adc_on_comb]])

module_dummy = cmrsim.bloch.GeneralBlochOperator(name="dummy_shot_module", gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
                                                 time_grid=t_comb_dummy[0], gradient_waveforms=grad_comb_dummy,
                                                 rf_waveforms=rf_comb_dummy, device="GPU:0")

module_acq = cmrsim.bloch.GeneralBlochOperator(name="testmodule", gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
                                               time_grid=t_comb[0], gradient_waveforms=grad_comb,
                                               rf_waveforms=rf_comb, adc_events=adc_on_comb, device="GPU:0")

[(51, 566), (51, 566), (51, 566, 3), (51, 566, 2)]

Run Simulation

print(f"Simulating {t_comb_dummy.shape[0] + t_comb.shape[0]*t_comb.shape[1]} steps for {m_init.shape[0]} isochromates")

for rbatch, mbatch, dict_batch in tqdm(tf.data.Dataset.from_tensor_slices((r, m_init, properties)).batch(int(1.5e6)), "Phantom Batches"):
    mbatch, rbatch = module_dummy(initial_position=rbatch, magnetization=mbatch, repetition_index=0, **dict_batch)

    for rep_idx in tqdm(range(t_comb.shape[0]), desc="TR index", leave=False):
        _ = module_acq(initial_position=rbatch, magnetization=mbatch, repetition_index=rep_idx, run_parallel=False, **dict_batch)

Simulating 28867 steps for 466820 isochromates

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1712826959.366766 2053534 device_compiler.h:186] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.

Reconstruct Image

# Non uniform combined
time_signal = tf.stack(module_acq.time_signal_acc, axis=0).numpy().reshape(matrix_size[::-1])
centered_projection = tf.signal.fft(time_signal).numpy()
centered_k_space = tf.signal.fft(centered_projection).numpy()
image = tf.signal.fftshift(tf.signal.ifft2d(tf.roll(tf.signal.ifftshift(centered_k_space, axes=(0, 1)),  -1, axis=1)), axes=(0, 1)).numpy()
time_per_seg = t_adc - t_center.m_as("ms")
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, centered_projection.shape, time_per_seg.shape, freq_per_seg.shape, image.shape)

(51, 51) (51, 51) (51,) (51,) (51, 51)

Plot 2D representation of specific single line acquisitions

If the configured adc only acquires the exact k-space locations the acquired time signal equals the k-space line

# Plot 2D Lines
fig, axes = plt.subplots(1, 4, figsize=(14, 2.5))
[a.grid(True) for a in axes]
[a.set_title(t) for a, t in zip(axes, ["temporal signal", "specturm", "k-space line", "spatial projection"])]
[a.legend(["real", "Imaginary"]) for a in axes]
for ky in range(0, time_signal.shape[0], 10):
    axes[0].plot(time_per_seg, np.real(time_signal[ky]), color="C0", linestyle="-", alpha=0.7)
    axes[0].plot(time_per_seg, np.imag(time_signal[ky]), color="C1", linestyle="-", alpha=0.7)
    axes[1].plot(freq_per_seg, np.real(centered_projection[ky]), color="C0", linestyle="-", alpha=0.7)
    axes[1].plot(freq_per_seg, np.imag(centered_projection[ky]), color="C1", linestyle="-", alpha=0.7)
    axes[2].plot(np.arange(freq_per_seg.shape[0]), np.real(centered_k_space[ky]), color="C0", linestyle="-", alpha=0.7)
    axes[2].plot(np.arange(freq_per_seg.shape[0]), np.imag(centered_k_space[ky]), color="C1", linestyle="-", alpha=0.7)
    axes[3].plot(np.arange(freq_per_seg.shape[0]), np.real(image[ky]), color="C0", linestyle="-", alpha=0.7)
    axes[3].plot(np.arange(freq_per_seg.shape[0]), np.imag(image[ky]), color="C1", linestyle="-", alpha=0.7)
fig.tight_layout()

Plot 3D representation of simulation results

If the configured adc only acquires the exact k-space locations the acquired time signal equals the k-space line

# Plot all PE-lines in 3D
plt.close("all")
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig, axes = plt.subplots(1, 4, subplot_kw={'projection': '3d'}, figsize=(14, 4))
[a.set_title(t) for a, t in zip(axes, ["temporal signal", "specturm", "k-space line", "spatial projection"])]

for ky in range(time_signal.shape[0]):
    axes[0].plot(time_per_seg, np.ones_like(time_signal[ky].real) * ky, np.abs(time_signal[ky]), color="C0", linestyle="-", alpha=0.7)
    axes[1].plot(freq_per_seg, np.ones_like(centered_projection[ky].real) * ky, np.abs(centered_projection[ky]), color="C0", linestyle="-", alpha=0.7)
    axes[2].plot(np.arange(freq_per_seg.shape[0]), np.ones_like(centered_k_space[ky].real) * ky, np.abs(centered_k_space[ky]), color="C0", linestyle="-", alpha=0.7)
    axes[3].plot(np.arange(freq_per_seg.shape[0]), np.ones_like(image[ky].real) * ky, np.abs(image[ky]), color="C0", linestyle="-", alpha=0.7)
fig.tight_layout()

Plot Image

# Show images
fig, axes = plt.subplots(1, 3, figsize=(12, 5))

kspace_plot = axes[0].imshow(np.log(np.abs(np.squeeze(centered_k_space))), cmap="gray")
fig.colorbar(kspace_plot, ax=axes[0], fraction=0.045, pad=0.04)

abs_plot = axes[1].imshow(np.abs(np.squeeze(image)), cmap="gray")
fig.colorbar(abs_plot, ax=axes[1], fraction=0.045, pad=0.04)
phase_plot = axes[2].imshow(np.angle(np.squeeze(image)), cmap="twilight", vmin=-np.pi, vmax=np.pi)
fig.colorbar(phase_plot, ax=axes[2], fraction=0.045, pad=0.04)
[_.axis("off") for _ in axes]
[_.set_title(t) for _, t in zip(axes, ["k-space", "Magnitude", "Phase"])]
fig.tight_layout()

Create Thumbnail

plt.close("all")

fig = plt.figure(constrained_layout=True, figsize=(15, 12))
gs = fig.add_gridspec(4, 6)
fig.add_subplot(gs[0, 0:3])
fig.add_subplot(gs[0, 3:])
fig.add_subplot(gs[1, 0:])
fig.add_subplot(gs[2:, 0:2])
fig.add_subplot(gs[2:, 2:4])
fig.add_subplot(gs[2:, 4:6])

axes = fig.axes
[a.set_title(t) for a, t in zip(axes, ['Sequence Per TR', 'k-space-vectors', 'Full Sequence', "k-space", "Magnitude", "Phase"])]

cmrseq.plotting.plot_sequence(sequence_list[0],  axes[0])
for tr_idx, s in enumerate(sequence_list[1:]):
    cmrseq.plotting.plot_sequence(s,  [fig.axes[-1], axes[0], axes[0], axes[0]], format_axes=False)
    if tr_idx > n_dummy:
        cmrseq.plotting.plot_kspace_2d(s, ax=axes[1], k_axes=[0, 1])
_ = cmrseq.plotting.plot_sequence(combined_sequence, axes=axes[2], format_axes=True)

kspace_plot = axes[3].imshow(np.log(np.abs(np.squeeze(centered_k_space))), cmap="gray")
fig.colorbar(kspace_plot, ax=axes[3], fraction=0.045, pad=0.04, location="bottom")
abs_plot = axes[4].imshow(np.abs(np.squeeze(image)), cmap="gray")
fig.colorbar(abs_plot, ax=axes[4], fraction=0.045, pad=0.04, location="bottom")
phase_plot = axes[5].imshow(np.angle(np.squeeze(image)), cmap="twilight", vmin=-np.pi, vmax=np.pi)
fig.colorbar(phase_plot, ax=axes[5], fraction=0.045, pad=0.04, location="bottom")
[_.axis("off") for _ in axes[3:]]
fig.savefig("static_bloc_spoiled_gre.svg")