Compare Fourier Model and T2* Model for EPI

Note

Go to the end to download the full example code.

Compare Fourier Model and T2* Model for EPI#

This examples walks through the elementary components of SNAKE.

Here we proceed step by step and use the Python interface. A more integrated alternative is to use the CLI snake-main

# Imports
import numpy as np
from snake.core.simulation import SimConfig, default_hardware, GreConfig
from snake.core.phantom import Phantom
from snake.core.smaps import get_smaps
from snake.core.sampling import EPI3dAcquisitionSampler

sim_conf = SimConfig(
    max_sim_time=6,
    seq=GreConfig(TR=100, TE=30, FA=3),
    hardware=default_hardware,
    fov_mm=(181, 217, 181),
    shape=(60, 72, 60),
)
sim_conf.hardware.n_coils = 8

phantom = Phantom.from_brainweb(sub_id=4, sim_conf=sim_conf, tissue_file="tissue_7T")

Setting up Acquisition Pattern and Initializing Result file.#

# The next piece of simulation is the acquisition trajectory.
# Here nothing fancy, we are using a EPI (fully sampled), that samples a 3D
# k-space (this akin to the 3D EPI sequence of XXXX)

sampler = EPI3dAcquisitionSampler(accelz=1, acsz=0.1, orderz="top-down")

smaps = None
if sim_conf.hardware.n_coils > 1:
    smaps = get_smaps(sim_conf.shape, n_coils=sim_conf.hardware.n_coils)

Acquisition with Cartesian Engine#

The generated file example_EPI.mrd does not contains any k-space data for now, only the sampling trajectory. let’s put some in. In order to do so, we need to setup the acquisition engine that models the MR physics, and get sampled at the specified k-space trajectory.

SNAKE comes with two models for the MR Physics:

  • model=”simple” :: Each k-space shot acquires a constant signal, which is the image contrast at TE.

  • model=”T2s” :: Each k-space shot is degraded by the T2* decay induced by each tissue.

# Here we will use the "simple" model, which is faster.
#
# SNAKE's Engine are capable of simulating the data in parallel, by distributing
# the shots to be acquired to a set of processes. To do so , we need to specify
# the number of jobs that will run in parallel, as well as the size of a job.
# Setting the job size and the number of jobs can have a great impact on total
# runtime and memory consumption.
#
# Here, we have a single frame to acquire with 60 frames (one EPI per slice), so
# a single worker will do.

from snake.core.engine import EPIAcquisitionEngine

engine = EPIAcquisitionEngine(model="simple")

engine("example_EPI.mrd", sampler, phantom, sim_conf, smaps=smaps)

engine(
    "example_EPI.mrd",
    sampler,
    phantom,
    sim_conf,
    smaps=smaps,
    worker_chunk_size=60,
    n_workers=1,
)
engine_t2s = EPIAcquisitionEngine(model="T2s")

engine_t2s(
    "example_EPI_t2s.mrd",
    sampler,
    phantom,
    sim_conf,
    smaps=smaps,
    worker_chunk_size=60,
    n_workers=1,
)

Existing example_EPI.mrd it will be overwritten

  0%|          | 0/60 [00:00<?, ?it/s]
  2%|▏         | 1/60 [00:00<00:17,  3.33it/s]
  5%|▌         | 3/60 [00:00<00:06,  8.19it/s]
  8%|▊         | 5/60 [00:00<00:04, 11.10it/s]
 12%|█▏        | 7/60 [00:00<00:04, 12.67it/s]
 15%|█▌        | 9/60 [00:00<00:03, 13.97it/s]
 18%|█▊        | 11/60 [00:00<00:03, 14.81it/s]
 22%|██▏       | 13/60 [00:01<00:03, 15.29it/s]
 25%|██▌       | 15/60 [00:01<00:02, 15.73it/s]
 28%|██▊       | 17/60 [00:01<00:02, 15.78it/s]
 32%|███▏      | 19/60 [00:01<00:02, 16.04it/s]
 35%|███▌      | 21/60 [00:01<00:02, 16.25it/s]
 38%|███▊      | 23/60 [00:01<00:02, 16.38it/s]
 42%|████▏     | 25/60 [00:01<00:02, 16.47it/s]
 45%|████▌     | 27/60 [00:01<00:01, 16.56it/s]
 48%|████▊     | 29/60 [00:01<00:01, 16.63it/s]
 52%|█████▏    | 31/60 [00:02<00:01, 16.67it/s]
 55%|█████▌    | 33/60 [00:02<00:01, 16.70it/s]
 58%|█████▊    | 35/60 [00:02<00:01, 16.31it/s]
 62%|██████▏   | 37/60 [00:02<00:01, 16.45it/s]
 65%|██████▌   | 39/60 [00:02<00:01, 16.55it/s]
 68%|██████▊   | 41/60 [00:02<00:01, 16.48it/s]
 72%|███████▏  | 43/60 [00:02<00:01, 16.50it/s]
 75%|███████▌  | 45/60 [00:02<00:00, 16.56it/s]
 78%|███████▊  | 47/60 [00:03<00:00, 16.56it/s]
 82%|████████▏ | 49/60 [00:03<00:00, 16.62it/s]
 85%|████████▌ | 51/60 [00:03<00:00, 16.19it/s]
 88%|████████▊ | 53/60 [00:03<00:00, 16.34it/s]
 92%|█████████▏| 55/60 [00:03<00:00, 16.45it/s]
 95%|█████████▌| 57/60 [00:03<00:00, 16.53it/s]
 98%|█████████▊| 59/60 [00:03<00:00, 16.58it/s]
100%|██████████| 60/60 [00:03<00:00, 15.50it/s]
No coil_cov found in the dataset.

  0%|          | 0/60 [00:00<?, ?it/s]
100%|██████████| 60/60 [00:03<00:00, 19.39it/s]
100%|██████████| 60/60 [00:03<00:00, 17.26it/s]
Existing example_EPI.mrd it will be overwritten

  0%|          | 0/60 [00:00<?, ?it/s]
  3%|▎         | 2/60 [00:00<00:03, 14.54it/s]
  7%|▋         | 4/60 [00:00<00:03, 15.80it/s]
 10%|█         | 6/60 [00:00<00:03, 16.16it/s]
 13%|█▎        | 8/60 [00:00<00:03, 16.28it/s]
 17%|█▋        | 10/60 [00:00<00:03, 16.34it/s]
 20%|██        | 12/60 [00:00<00:02, 16.35it/s]
 23%|██▎       | 14/60 [00:00<00:02, 16.39it/s]
 27%|██▋       | 16/60 [00:00<00:02, 16.40it/s]
 30%|███       | 18/60 [00:01<00:02, 16.14it/s]
 33%|███▎      | 20/60 [00:01<00:02, 16.22it/s]
 37%|███▋      | 22/60 [00:01<00:02, 16.29it/s]
 40%|████      | 24/60 [00:01<00:02, 16.33it/s]
 43%|████▎     | 26/60 [00:01<00:02, 16.36it/s]
 47%|████▋     | 28/60 [00:01<00:01, 16.38it/s]
 50%|█████     | 30/60 [00:01<00:01, 16.41it/s]
 53%|█████▎    | 32/60 [00:01<00:01, 16.43it/s]
 57%|█████▋    | 34/60 [00:02<00:01, 16.10it/s]
 60%|██████    | 36/60 [00:02<00:01, 15.99it/s]
 63%|██████▎   | 38/60 [00:02<00:01, 16.16it/s]
 67%|██████▋   | 40/60 [00:02<00:01, 16.28it/s]
 70%|███████   | 42/60 [00:02<00:01, 16.38it/s]
 73%|███████▎  | 44/60 [00:02<00:00, 16.48it/s]
 77%|███████▋  | 46/60 [00:02<00:00, 16.64it/s]
 80%|████████  | 48/60 [00:02<00:00, 16.74it/s]
 83%|████████▎ | 50/60 [00:03<00:00, 16.77it/s]
 87%|████████▋ | 52/60 [00:03<00:00, 16.52it/s]
 90%|█████████ | 54/60 [00:03<00:00, 16.68it/s]
 93%|█████████▎| 56/60 [00:03<00:00, 16.75it/s]
 97%|█████████▋| 58/60 [00:03<00:00, 16.83it/s]
100%|██████████| 60/60 [00:03<00:00, 16.76it/s]
100%|██████████| 60/60 [00:03<00:00, 16.41it/s]
No coil_cov found in the dataset.

  0%|          | 0/60 [00:00<?, ?it/s]"Unable to synchronously open object (object 'waveforms' doesn't exist)"

100%|██████████| 60/60 [00:01<00:00, 32.35it/s]
100%|██████████| 60/60 [00:02<00:00, 27.22it/s]
Existing example_EPI_t2s.mrd it will be overwritten
[Errno 2] No such file or directory: 'example_EPI_t2s.mrd'

  0%|          | 0/60 [00:00<?, ?it/s]
  3%|▎         | 2/60 [00:00<00:03, 15.73it/s]
  7%|▋         | 4/60 [00:00<00:03, 16.29it/s]
 10%|█         | 6/60 [00:00<00:03, 16.50it/s]
 13%|█▎        | 8/60 [00:00<00:03, 16.63it/s]
 17%|█▋        | 10/60 [00:00<00:02, 16.67it/s]
 20%|██        | 12/60 [00:00<00:02, 16.70it/s]
 23%|██▎       | 14/60 [00:00<00:02, 16.71it/s]
 27%|██▋       | 16/60 [00:00<00:02, 16.73it/s]
 30%|███       | 18/60 [00:01<00:02, 16.45it/s]
 33%|███▎      | 20/60 [00:01<00:02, 16.54it/s]
 37%|███▋      | 22/60 [00:01<00:02, 16.63it/s]
 40%|████      | 24/60 [00:01<00:02, 16.68it/s]
 43%|████▎     | 26/60 [00:01<00:02, 16.72it/s]
 47%|████▋     | 28/60 [00:01<00:01, 16.73it/s]
 50%|█████     | 30/60 [00:01<00:01, 16.75it/s]
 53%|█████▎    | 32/60 [00:01<00:01, 16.77it/s]
 57%|█████▋    | 34/60 [00:02<00:01, 16.41it/s]
 60%|██████    | 36/60 [00:02<00:01, 16.50it/s]
 63%|██████▎   | 38/60 [00:02<00:01, 16.59it/s]
 67%|██████▋   | 40/60 [00:02<00:01, 16.58it/s]
 70%|███████   | 42/60 [00:02<00:01, 16.65it/s]
 73%|███████▎  | 44/60 [00:02<00:00, 16.71it/s]
 77%|███████▋  | 46/60 [00:02<00:00, 16.72it/s]
 80%|████████  | 48/60 [00:02<00:00, 16.75it/s]
 83%|████████▎ | 50/60 [00:03<00:00, 16.78it/s]
 87%|████████▋ | 52/60 [00:03<00:00, 16.35it/s]
 90%|█████████ | 54/60 [00:03<00:00, 16.49it/s]
 93%|█████████▎| 56/60 [00:03<00:00, 16.57it/s]
 97%|█████████▋| 58/60 [00:03<00:00, 16.61it/s]
100%|██████████| 60/60 [00:03<00:00, 16.62it/s]
100%|██████████| 60/60 [00:03<00:00, 16.61it/s]
No coil_cov found in the dataset.

  0%|          | 0/60 [00:00<?, ?it/s]"Unable to synchronously open object (object 'waveforms' doesn't exist)"

100%|██████████| 60/60 [00:06<00:00,  9.95it/s]
100%|██████████| 60/60 [00:06<00:00,  9.37it/s]

Simple reconstruction#

Getting k-space data is nice, but SNAKE also provides rudimentary reconstruction tools to get images (and check that we didn’t mess up the acquisition process). This is available in the companion package snake.toolkit.

Loading the .mrd file to retrieve all information can be done using the ismrmd python package, but SNAKE provides convenient dataloaders, which are more efficient, and take cares of managing underlying files access. As we are showcasing the API, we will do things manually here, and use only core SNAKE.

from snake.mrd_utils import CartesianFrameDataLoader
from scipy.fft import ifftn, ifftshift, fftshift


def reconstruct_frame(filename):
    with CartesianFrameDataLoader(filename) as data_loader:
        mask, kspace_data = data_loader.get_kspace_frame(0)
    axes = (-3, -2, -1)
    image_data = ifftshift(
        ifftn(fftshift(kspace_data, axes=axes), axes=axes, norm="ortho"), axes=axes
    )

    # Take the square root sum of squares to get the magnitude image (SSOS)
    image_data = np.sqrt(np.sum(np.abs(image_data) ** 2, axis=0))

    return image_data.squeeze().T


image_simple = reconstruct_frame("example_EPI.mrd")
image_T2s = reconstruct_frame("example_EPI_t2s.mrd")

Plotting the result#

import matplotlib.pyplot as plt
from snake.toolkit.plotting import axis3dcut

fig, axs = plt.subplots(1, 3, figsize=(30, 10))

for ax, img, title in zip(
    axs,
    (image_simple, image_T2s, abs(image_simple - image_T2s)),
    ("simple", "T2s", "diff"),
):
    axis3dcut(fig, ax, img, None, None, cbar=True, cuts=(40, 40, 40), width_inches=4)
    ax.set_title(title)

plt.show()
simple, T2s, diff

Total running time of the script: (0 minutes 31.371 seconds)

Gallery generated by Sphinx-Gallery

Contents