Note
Go to the end to download the full example code.
Single anatomical volume with 2D EPI with SNAKE-fMRI#
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
Setting up the base simulation Config. This configuration holds all key parameters for the simulation, describing the scanner parameters.
sim_conf = SimConfig(
max_sim_time=3,
seq=GreConfig(TR=50, TE=30, FA=3),
hardware=default_hardware,
)
sim_conf.hardware.n_coils = 8
sim_conf.fov.res_mm = (3,3,3)
sim_conf
Creating the base Phantom#
The simulation acquires the data describe in a phantom. A phantom consists of fuzzy segmentation of head tissue, and their MR intrinsic parameters (density, T1, T2, T2*, magnetic susceptibilities)
Here we use Brainweb reference mask and values for convenience.
phantom = Phantom.from_brainweb(sub_id=4, sim_conf=sim_conf, output_res=1)
# Here are the tissue availables and their parameters
phantom
No matrix size found in the header. The header is probably missing.
Phantom[brainweb-{sub_id:02d}]: (12, 5)
tissue name PropTissueEnum.T1PropTissueEnum.T2PropTissueEnum.T2sPropTissueEnum.rhoPropTissueEnum.chi
vessels 1600.080.070.0 1.0-9.0
wm 500.070.061.00.7699999809265137-9.050000190734863
muscles 900.050.040.0 1.0-9.050000190734863
skull 300.0 1.0 1.00.10000000149011612-8.859999656677246
dura 1000.0100.050.00.6000000238418579-9.050000190734863
marrow 1000.050.040.0 1.0-9.050000190734863
csf 4000.02000.02000.0 1.0-9.050000190734863
muscles_skin 2569.0329.058.0 1.0-9.050000190734863
gm 900.090.069.00.8600000143051147-9.050000190734863
fat 250.070.060.0 1.0-7.789999961853027
bck 800.0800.0800.0 0.00.36000001430511475
fat2 250.070.060.0 1.0-7.789999961853027
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.
Here, we have a single frame to acquire with 60 slice (one EPI per slice), so a single worker will be enough.
from snake.core.engine import EPIAcquisitionEngine
engine = EPIAcquisitionEngine(model="T2s", slice_2d=True) # use the 2D epi.
engine(
"example_EPI2D.mrd",
sampler=sampler,
phantom=phantom,
sim_conf=sim_conf,
worker_chunk_size=20,
n_workers=2,
)
Using 2D acquisition, the TR_eff is updated to TR_vol
Existing example_EPI2D.mrd it will be overwritten
[Errno 2] No such file or directory: 'example_EPI2D.mrd'
0%| | 0/60 [00:00<?, ?it/s]
2%|▏ | 1/60 [00:00<00:25, 2.35it/s]
3%|▎ | 2/60 [00:00<00:14, 4.06it/s]
5%|▌ | 3/60 [00:00<00:10, 5.28it/s]
7%|▋ | 4/60 [00:00<00:09, 6.15it/s]
8%|▊ | 5/60 [00:00<00:08, 6.77it/s]
10%|█ | 6/60 [00:01<00:07, 7.21it/s]
12%|█▏ | 7/60 [00:01<00:07, 7.52it/s]
13%|█▎ | 8/60 [00:01<00:06, 7.70it/s]
15%|█▌ | 9/60 [00:01<00:06, 7.87it/s]
17%|█▋ | 10/60 [00:01<00:06, 7.99it/s]
18%|█▊ | 11/60 [00:01<00:06, 8.05it/s]
20%|██ | 12/60 [00:01<00:05, 8.11it/s]
22%|██▏ | 13/60 [00:01<00:05, 8.16it/s]
23%|██▎ | 14/60 [00:02<00:05, 8.20it/s]
25%|██▌ | 15/60 [00:02<00:05, 8.12it/s]
27%|██▋ | 16/60 [00:02<00:05, 8.17it/s]
28%|██▊ | 17/60 [00:02<00:05, 8.19it/s]
30%|███ | 18/60 [00:02<00:05, 8.22it/s]
32%|███▏ | 19/60 [00:02<00:04, 8.24it/s]
33%|███▎ | 20/60 [00:02<00:04, 8.26it/s]
35%|███▌ | 21/60 [00:02<00:04, 8.25it/s]
37%|███▋ | 22/60 [00:02<00:04, 8.26it/s]
38%|███▊ | 23/60 [00:03<00:04, 8.27it/s]
40%|████ | 24/60 [00:03<00:04, 8.27it/s]
42%|████▏ | 25/60 [00:03<00:04, 8.26it/s]
43%|████▎ | 26/60 [00:03<00:04, 8.27it/s]
45%|████▌ | 27/60 [00:03<00:03, 8.26it/s]
47%|████▋ | 28/60 [00:03<00:03, 8.27it/s]
48%|████▊ | 29/60 [00:03<00:03, 8.25it/s]
50%|█████ | 30/60 [00:03<00:03, 8.05it/s]
52%|█████▏ | 31/60 [00:04<00:03, 8.06it/s]
53%|█████▎ | 32/60 [00:04<00:03, 8.06it/s]
55%|█████▌ | 33/60 [00:04<00:03, 8.05it/s]
57%|█████▋ | 34/60 [00:04<00:03, 8.11it/s]
58%|█████▊ | 35/60 [00:04<00:03, 8.16it/s]
60%|██████ | 36/60 [00:04<00:02, 8.09it/s]
62%|██████▏ | 37/60 [00:04<00:02, 8.14it/s]
63%|██████▎ | 38/60 [00:04<00:02, 8.16it/s]
65%|██████▌ | 39/60 [00:05<00:02, 8.13it/s]
67%|██████▋ | 40/60 [00:05<00:02, 8.15it/s]
68%|██████▊ | 41/60 [00:05<00:02, 8.13it/s]
70%|███████ | 42/60 [00:05<00:02, 8.17it/s]
72%|███████▏ | 43/60 [00:05<00:02, 8.15it/s]
73%|███████▎ | 44/60 [00:05<00:01, 8.17it/s]
75%|███████▌ | 45/60 [00:05<00:01, 8.17it/s]
77%|███████▋ | 46/60 [00:05<00:01, 8.19it/s]
78%|███████▊ | 47/60 [00:06<00:01, 8.14it/s]
80%|████████ | 48/60 [00:06<00:01, 8.18it/s]
82%|████████▏ | 49/60 [00:06<00:01, 8.15it/s]
83%|████████▎ | 50/60 [00:06<00:01, 8.16it/s]
85%|████████▌ | 51/60 [00:06<00:01, 8.20it/s]
87%|████████▋ | 52/60 [00:06<00:00, 8.04it/s]
88%|████████▊ | 53/60 [00:06<00:00, 8.11it/s]
90%|█████████ | 54/60 [00:06<00:00, 8.15it/s]
92%|█████████▏| 55/60 [00:07<00:00, 8.18it/s]
93%|█████████▎| 56/60 [00:07<00:00, 8.21it/s]
95%|█████████▌| 57/60 [00:07<00:00, 8.24it/s]
97%|█████████▋| 58/60 [00:07<00:00, 8.25it/s]
98%|█████████▊| 59/60 [00:07<00:00, 8.27it/s]
100%|██████████| 60/60 [00:07<00:00, 8.26it/s]
100%|██████████| 60/60 [00:07<00:00, 7.87it/s]
No coil_cov found in the dataset.
0%| | 0/60 [00:00<?, ?it/s]
33%|███▎ | 20/60 [00:19<00:39, 1.02it/s]
67%|██████▋ | 40/60 [00:20<00:08, 2.37it/s]
100%|██████████| 60/60 [00:32<00:00, 1.93it/s]
100%|██████████| 60/60 [00:33<00:00, 1.80it/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
with CartesianFrameDataLoader("example_EPI2D.mrd") as data_loader:
mask, kspace_data = data_loader.get_kspace_frame(0)
Reconstructing a Single Frame of fully sampled EPI boils down to performing a 3D IFFT:
from scipy.fft import ifftn, ifftshift, fftshift
axes = (-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))
import matplotlib.pyplot as plt
from snake.toolkit.plotting import axis3dcut
fig, ax = plt.subplots()
axis3dcut(image_data.squeeze().T, None, None, cbar=False, cuts=(0.5,0.5,0.5), ax=ax)
plt.show()
Total running time of the script: (0 minutes 44.632 seconds)