Static Phantom with off-resonance

In this notebook, a simulation of a balanced steady state free precession with a spherical static phantom with air inclusions is performed. To do so the following steps are performed:

1. Define the Sequence (with cmrseq)

2. Define a 3D spherical phantom

3. Run simulation

4. Reconstruct images

Import modules and set TF-GPU configuration

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

2024-04-11 11:57:09.786775: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX512F FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2024-04-11 11:57:12.741889: I external/local_xla/xla/stream_executor/cuda/cuda_executor.cc:887] could not open file to read NUMA node: /sys/bus/pci/devices/0000:b2:00.0/numa_node
Your kernel may have been built without NUMA support.
2024-04-11 11:57:13.271896: I external/local_xla/xla/stream_executor/cuda/cuda_executor.cc:887] could not open file to read NUMA node: /sys/bus/pci/devices/0000:b2:00.0/numa_node
Your kernel may have been built without NUMA support.
2024-04-11 11:57:13.272238: I external/local_xla/xla/stream_executor/cuda/cuda_executor.cc:887] could not open file to read NUMA node: /sys/bus/pci/devices/0000:b2:00.0/numa_node
Your kernel may have been built without NUMA support.

import sys
sys.path.insert(0, "../../")
from copy import deepcopy
import base64
from IPython.display import display, HTML, 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
import cmrseq

sys.path.append("../")
import local_functions
clear_output()

1. Define Sequence

Configure Simulation Parameters

# 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"),
                                 b0=Quantity(1.5, "T")
                                )


# Define Fourier-Encoding parameters
fov = Quantity([21, 21], "cm")
matrix_size = (151, 151)
res = fov / matrix_size
print("Resolution:", res)
n_dummy = 51
center_index = np.array(matrix_size) // 2 + np.mod(matrix_size, 2)

# Define imaging slice parameters
slice_thickness = Quantity(5, "mm")
slice_position = Quantity([0., 0., 0.], "m")
adc_duration = Quantity(2., "ms")
readout_direction = np.array([0., 1., 0.])
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/4, "rad")
pulse_duration = Quantity(0.6, "ms")

Resolution: [0.1390728476821192 0.1390728476821192] centimeter

Instantiate Sequence

Note: As a \(\alpha/2-TR/2\) preparation module is used, the first acquisition happens in the n_dummy+1 entry of the sequence list

seq_list = cmrseq.seqdefs.sequences.balanced_ssfp(system_specs,
                                                  matrix_size,
                                                  inplane_resolution=res,
                                                  repetition_time=Quantity(0, "ms"),
                                                  slice_thickness=slice_thickness,
                                                  adc_duration=adc_duration,
                                                  flip_angle=flip_angle,
                                                  pulse_duration=pulse_duration,
                                                  slice_position_offset=slice_position_offset,
                                                  dummy_shots=n_dummy)

CMRSeq Warning : bSSFP Sequence: Repetition time too short to be feasible, set TR to 4.2 millisecond

Plot sequence

_, k_adc, t_adc = seq_list[n_dummy + center_index[1]].calculate_kspace()
full_seq = deepcopy(seq_list[0])
full_seq.extend(deepcopy(seq_list[1:]))
clear_output()

plt.close("all")
f = plt.figure(constrained_layout=True, figsize=(15, 8))
axes = f.subplot_mosaic("AAB;CCC")

cmrseq.plotting.plot_sequence(seq_list[0], axes=axes["A"], format_axes=True, add_legend=False)
for s in seq_list[n_dummy:]:
    cmrseq.plotting.plot_sequence(s, axes=[f.axes[-1], axes["A"], axes["A"], axes["A"]], format_axes=False)
for s in seq_list[n_dummy:]:
    cmrseq.plotting.plot_kspace_2d(s, ax=axes["B"], k_axes=[0, 1])
cmrseq.plotting.plot_sequence(full_seq, axes=axes["C"])

2. Define Phantom

Create 3D object

phantom = local_functions.create_spherical_phantom(dimensions=(61, 61, 61), spacing=(0.004, 0.004, 0.004))
clear_output()

# Create field to signal if regular mesh-points are within original unstructured grid
phantom["in_mesh"] = np.ones(phantom.points.shape[0])

Instantiate a RegularGridDataset

Also assign MRI relevant material properties and comput the off-resonance distribution

# Create regular grids
dataset = cmrsim.datasets.RegularGridDataset.from_unstructured_grid(phantom, pixel_spacing=Quantity([0.5, 0.5, 0.5], "mm"), padding=Quantity([5, 5, 5], "cm"))

# Assign properties (everything not 'in-mesh' is assumed to be air)
dataset.mesh["chi"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]) * (-9.05), np.ones_like(dataset.mesh["in_mesh"]) * 0.36) * 1e-6   # susceptibilty in ppm
dataset.mesh["M0"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]), np.zeros_like(dataset.mesh["in_mesh"]))          # density in percent
dataset.mesh["T1"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]) * 1000, np.zeros_like(dataset.mesh["in_mesh"]))   # time in ms
dataset.mesh["T2"] = np.where(dataset.mesh["in_mesh"], np.ones_like(dataset.mesh["in_mesh"]) * 100, np.zeros_like(dataset.mesh["in_mesh"]))    # time in ms

# Calculate b0-offresonance field
dataset.compute_offresonance(b0=system_specs.b0, susceptibility_key="chi");

CMRSeq RuntimeWarning: /scratch/jweine/cmrsim/notebooks/bloch_simulation/../../cmrsim/datasets/_regular_grid.py:192 : invalid value encountered in divide
CMRSeq UnitStrippedWarning: /usr/local/lib/python3.11/dist-packages/pyvista/core/dataset.py:1975 : The unit of the quantity is stripped when downcasting to ndarray.

Select a slice

The slice coordinates are transformed into MPS coordinates as visualized below. The selected slice-mesh is wrapped again into a RegularGridDataset to subsequently extract the non-trivial material points as dictionary of arrays.

slice_mesh = dataset.select_slice(slice_normal=slice_normal, spacing=Quantity((0.25, 0.25, 1.5), "mm"),
                                  slice_position=slice_normal*slice_position_offset,
                                  field_of_view=Quantity([22, 22, slice_thickness.m_as("cm")*1.5], "cm"),
                                  readout_direction=readout_direction, in_mps=True)
slice_dataset = cmrsim.datasets.RegularGridDataset(slice_mesh)
phantom_dict = slice_dataset.get_phantom_def(filter_by="in_mesh", keys=["M0", "T1", "T2", "offres"],
                                             perturb_positions_std=Quantity(0.01, "mm").m_as("m"))

CMRSeq Warning : Optimal rotation is not uniquely or poorly defined for the given sets of vectors.

Plot 3D phantom with offresonance map

pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(window_size=(1200, 400), shape=(1, 2))
plotter.add_mesh(dataset.mesh.threshold(0.1, "M0"), scalars="offres", opacity=0.65,
                 cmap="twilight", scalar_bar_args=dict(title="Offresonance (T)"))
local_functions.add_custom_axes(plotter, scale=0.15)
box = local_functions.transformed_box(np.eye(3, 3), slice_normal, readout_direction, slice_position_offset*slice_normal,
                                      extend=Quantity([*fov.m_as("m"), slice_thickness.m_as("m")], "m"))
plotter.add_mesh(box)
plotter.subplot(0, 1)
plotter.add_mesh(pyvista.PolyData(phantom_dict["r_vectors"]),
                 scalars=phantom_dict["offres"],  opacity=0.65,
                 cmap="twilight", scalar_bar_args=dict(title="Offresonance (T)"))
local_functions.add_custom_axes(plotter, scale=0.15)
plotter.show(jupyter_backend="static")

CMRSeq Warning : Optimal rotation is not uniquely or poorly defined for the given sets of vectors.

3. Perfom Simulation

Instantiate simulation modules

From the preparation phase as well as the acquisition TRs grid the sequence and instantiate the modules including the off-resonance submodule

dummy_sequence = deepcopy(seq_list[0])
dummy_sequence.extend(seq_list[1:n_dummy+1])


time_dummy, rf_grid_dummy, grad_grid_dummy, _ = [np.stack(v) for v in cmrseq.utils.grid_sequence_list([dummy_sequence, ])]
time, rf_grid, grad_grid, adc_on_grid = [np.stack(v) for v in cmrseq.utils.grid_sequence_list(seq_list[n_dummy+1:])]
print(time_dummy.shape, rf_grid_dummy.shape, grad_grid_dummy.shape)
print(time.shape, rf_grid.shape, grad_grid.shape, adc_on_grid.shape)

# Define off-resonsance submodule
offres = cmrsim.bloch.submodules.OffResonance(gamma=system_specs.gamma.m_as("MHz/T"), device="GPU:0")


# Construct BlochOperators and pass offresonance module:
module_dummyshots = cmrsim.bloch.GeneralBlochOperator(name="dummy_shots", 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=[offres]
                                                      )

module_acquisition = cmrsim.bloch.GeneralBlochOperator(name="acquisition", 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=[offres]
                                                     )

Extending Sequence: 100%|██████████| 51/51 [00:00<00:00, 210.32it/s]

(1, 21631) (1, 21631) (1, 21631, 3)
(151, 571) (151, 571) (151, 571, 3) (151, 571, 2)

Create Bloch-dataset for batched data stream

properties = {"M0":  phantom_dict["M0"].astype(np.float32),
              "T1":  phantom_dict["T1"].astype(np.float32),
              "T2":  phantom_dict["T2"].astype(np.float32),
              "off_res": (phantom_dict["offres"]).astype(np.float32).reshape(-1, 1),
              "magnetization": cmrsim.utils.particle_properties.norm_magnetization()(phantom_dict["M0"].shape[0]),
              "initial_position": phantom_dict["r_vectors"].astype(np.float32)
              }
input_dataset = cmrsim.datasets.BlochDataset(properties, filter_inputs=True)
print(properties["M0"].shape)

(2368744,)

Call modules in loop-structure

# Instantiate simulator:
print(f"Total time-steps per TR: {time.shape[1]}")

for batch in tqdm(input_dataset(batchsize=int(3e6)).take(1)):
    print({k:v.shape for k, v in batch.items()})
    m_init = batch.pop("magnetization")
    initial_position = batch.pop("initial_position")
    m, r = module_dummyshots(initial_position=initial_position, magnetization=m_init, **batch)
    # m, r = m_init, initial_position
    for tr_index in tqdm(range(time.shape[0]), desc="Iterating TRs", leave=False):
        m, r = module_acquisition(initial_position=r, magnetization=m, repetition_index=tr_index, **batch)

Total time-steps per TR: 571

{'M0': TensorShape([2368744]), 'T1': TensorShape([2368744]), 'T2': TensorShape([2368744]), 'off_res': TensorShape([2368744, 1]), 'magnetization': TensorShape([2368744, 3]), 'initial_position': TensorShape([2368744, 3])}

2024-04-11 12:07:00.557728: I external/local_xla/xla/service/service.cc:168] XLA service 0xd71df30 initialized for platform CUDA (this does not guarantee that XLA will be used). Devices:
2024-04-11 12:07:00.557798: I external/local_xla/xla/service/service.cc:176]   StreamExecutor device (0): NVIDIA TITAN RTX, Compute Capability 7.5
2024-04-11 12:07:00.641387: I tensorflow/compiler/mlir/tensorflow/utils/dump_mlir_util.cc:269] disabling MLIR crash reproducer, set env var `MLIR_CRASH_REPRODUCER_DIRECTORY` to enable.
2024-04-11 12:07:00.809058: I external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:454] Loaded cuDNN version 8904
2024-04-11 12:07:00.907477: I external/local_tsl/tsl/platform/default/subprocess.cc:304] Start cannot spawn child process: No such file or directory
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1712830021.791145 2063062 device_compiler.h:186] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.

4. Reconstruct images

time_signal = tf.stack(module_acquisition.time_signal_acc, axis=0).numpy().reshape(matrix_size[::-1], order="F")
time_signal += tf.complex(*[tf.random.normal(shape=time_signal.shape, stddev=50) for i in range(2)])
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()

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

# kspace_plot = axes[0].imshow(np.log(np.abs(np.squeeze(centered_k_space))), cmap="gray")
kspace_plot = axes[0].imshow(np.angle(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()