Spin Echo with Single Shot EPI

In this notebook the sequence is composed with cmrseq and the input is a shepp-logan phantom, constructed with the phantominator package link

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
import pyvista
%matplotlib inline

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

Construct Sequence

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

sequence_kwargs = dict(
    field_of_view=Quantity([200, 200], "mm"),
    matrix_size=(51, 51),
    echo_time=Quantity(40, "ms"),
    slice_thickness=Quantity(2, "cm"),
    slice_orientation=np.array([0, 0, 1.]),
    pulse_duration=Quantity(1500, "us"),
    epi_slope_sampling=False,
    tbw_product=4.,
    max_epi_duration=None,
    epi_water_fat_shift='minimum',
    partial_fourier_lines=0,
    blip_direction='up',
)

seq = cmrseq.seqdefs.sequences.se_ssepi(system_specs, **sequence_kwargs)

res = sequence_kwargs["field_of_view"] / sequence_kwargs["matrix_size"]

print("Resolution:", res)

Resolution: [3.9215686274509802 3.9215686274509802] millimeter

# Plot Sequence
plt.close("all")
f, axes = plt.subplots(1, 2, figsize=(15, 5), gridspec_kw={"width_ratios": (2, 1)})
cmrseq.plotting.plot_sequence(seq, axes=axes[0])
cmrseq.plotting.plot_kspace_2d(seq, ax=axes[1])

<Axes: xlabel='$k_x$ $[1/m]$', ylabel='$k_y$ $[1/m]$'>

Set up Phantom

import phantominator
x, y, z = 750, 750, 15
m0, t1, t2 = phantominator.shepp_logan(N=(x, y, z), MR=True, B0=1)

image_box = np.array((0.2, 0.2, 0.15))
dims = np.array([x, y, z])
resolution = image_box / dims
origin = - image_box / 2
uniform_grid = pyvista.ImageData(dimensions=dims, spacing=resolution, origin=origin)
uniform_grid["M0"] = m0.flatten(order="F")
uniform_grid["T1"] = t1.flatten(order="F") * 1000
uniform_grid["T2"] = t2.flatten(order="F") * 1000

uniform_grid

CMRSeq PyVistaDeprecationWarning: /usr/local/lib/python3.11/dist-packages/pyvista/core/grid.py:549 : `dims` argument is deprecated. Please use `dimensions`.

Header	Data Arrays
ImageData Information N Cells 7854014 N Points 8437500 X Bounds -1.000e-01, 9.973e-02 Y Bounds -1.000e-01, 9.973e-02 Z Bounds -7.500e-02, 6.500e-02 Dimensions 750, 750, 15 Spacing 2.667e-04, 2.667e-04, 1.000e-02 N Arrays 3	Name Field Type N Comp Min Max M0 Points float64 1 0.000e+00 9.800e-01 T1 Points float64 1 0.000e+00 4.200e+03 T2 Points float64 1 0.000e+00 1.990e+03

properties = {"M0":  uniform_grid["M0"].astype(np.float32),
              "T1":  uniform_grid["T1"].astype(np.float32),
              "T2":  uniform_grid["T2"].astype(np.float32),
              "magnetization": part_factory.norm_magnetization()(x*y*z),
              "initial_position": uniform_grid.points.astype(np.float32)
              }
input_dataset = cmrsim.datasets.BlochDataset(properties, filter_inputs=True)
surface_mesh = uniform_grid.threshold(0.1, scalars="M0")

Grid sequence and configure Bloch Module

t_comb, rf_comb, grad_comb, adc_on_comb = [np.stack(v) for v in cmrseq.utils.grid_sequence_list([seq, ])]

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

module_comb = 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")

[(1, 7943), (1, 7943), (1, 7943, 3), (1, 7943, 2)]

Run Simulation

for batch in tqdm(input_dataset(batchsize=int(1e6))):
    _ = module_comb(repetition_index=0, **batch)

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

Reconstruct Image

# Flip uneven k-space lines as they are acquired in "reversed" time
time_signal = np.stack(module_comb.time_signal_acc, axis=0).reshape(sequence_kwargs["matrix_size"][::-1])
time_signal[1::2, :] = time_signal[1::2, ::-1]
centered_k_space = time_signal
image = tf.signal.fftshift(tf.signal.ifft2d(tf.signal.ifftshift(centered_k_space))).numpy() * (-1)

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()