Echo Planar Imaging#
Qualitative imaging with an EPI readout.
Imports#
import importlib
if not importlib.util.find_spec('mrseq'):
%pip install mrseq[examples]
import tempfile
from pathlib import Path
import matplotlib.pyplot as plt
import MRzeroCore as mr0
import numpy as np
from mrpro.algorithms.reconstruction import DirectReconstruction
from mrpro.data import KData
from mrpro.data.traj_calculators import KTrajectoryIsmrmrd
from raw2ismrmrd.utils import combine_ismrmrd_files
from mrseq.sequences.epi2d_fid import main as create_seq_fid
from mrseq.sequences.epi2d_se import main as create_seq_se
from mrseq.utils import sys_defaults
from mrseq.utils.EpiReadout import EpiReadout
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
Settings#
We are going to use a numerical phantom with a matrix size of 128 x 128.
image_matrix_size = [100, 100]
flip_angle_degree = 12
n_dummy_excitations = 200
tmp = tempfile.TemporaryDirectory()
Create the digital phantom#
We use the standard Brainweb phantom from MRzero, but we choose the B1-field to be constant everywhere.
phantom = mr0.util.load_phantom(image_matrix_size)
phantom.B1[:] = 1.0
Downloaded simulation data
EPI readouts#
Single-shot EPI readouts are a very efficient way to obtain a full image after a single excitation. There are different options of how they can be carried out.
We can have a symmetric readout where we alternate between going from left to right and right to left through k-space.
epi_readout = EpiReadout(n_readout=12, n_phase_encoding=12, ramp_sampling=False, readout_type='symmetric')
fig = epi_readout.plot_trajectory()
Here we are sampling only during the flat time of the readout gradient. We can make the sequence faster by also sampling during the ramp-up and ramp-down part of the readout gradient. Nevertheless, this means that we don’t have a strictly Cartesian trajectory anymore and gridding/nufft is required during image reconstruction.
epi_readout = EpiReadout(n_readout=12, n_phase_encoding=12, ramp_sampling=True, readout_type='symmetric')
fig = epi_readout.plot_trajectory()
We can shorten the echo time by applying partial Fourier along the phase encoding direction.
epi_readout = EpiReadout(
n_readout=12, n_phase_encoding=12, ramp_sampling=True, partial_fourier_factor=0.75, readout_type='symmetric'
)
fig = epi_readout.plot_trajectory()
Another approach is to always go back to one side of k-space after each readout such that the readout direction is always from left to right. This is less efficient but minimized any artifacts due to asymmetry between positive and negative readout gradients.
epi_readout = EpiReadout(
n_readout=12, n_phase_encoding=12, partial_fourier_factor=0.75, ramp_sampling=True, readout_type='flyback'
)
fig = epi_readout.plot_trajectory()
EPI sequences#
To create the EPI FID and SE sequences, we use the previously imported epi2d_fid script and epi2d_se script.
We are going to simulate different types of EPI sequences and compare the obtained images.
EPI FID with symmetric readout without ramp sampling#
sequence, fname_seq = create_seq_fid(
system=sys_defaults,
test_report=False,
timing_check=False,
show_plots=False,
te=0.073,
fov_xy=float(phantom.size.numpy()[0]),
n_readout=image_matrix_size[0],
n_phase_encoding=image_matrix_size[0],
partial_fourier_factor=0.75,
ramp_sampling=False,
readout_type='symmetric',
)
mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))
signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)
fname_mrd = Path(tmp.name) / 'epi_fid_sym_no_ramp.mrd'
mr0.sig_to_mrd(fname_mrd, signal, sequence)
combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))
kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())
recon = DirectReconstruction(kdata, csm=None)
idata_fid_sym_no_ramp = recon(kdata)
Saving sequence file 'epi2d_fid_200fov_100px_sym_fid_2ro_nors_0p75pf_withnoise_withnav.seq' into folder '/home/runner/work/mrseq/mrseq/examples/output'.
>>>> Rust - compute_graph(...) >>>
Converting Python -> Rust: 0.000059608 s
Compute Graph
Computing Graph: 0.001833782 s
Analyze Graph
Analyzing Graph: 0.000001241 s
Converting Rust -> Python: 0.000008743 s
<<<< Rust <<<<
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/data/KData.py:166: UserWarning: No vendor information found. Assuming Siemens time stamp format.
warnings.warn('No vendor information found. Assuming Siemens time stamp format.', stacklevel=1)
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/operators/FourierOp.py:124: UserWarning: If both FFT and NUFFT dims are present, Cartesian FFT dims need to be aligned with the k-space dimension, i.e. kx along k0, ky along k1 and kz along k2. We are going to use NUFFT for all dimensions. Creating your own FourierOp with the desired combination of FFT and NUFFT dimensions might be more efficient.
return cls(
EPI FID with symmetric readout and ramp sampling#
sequence, fname_seq = create_seq_fid(
system=sys_defaults,
test_report=False,
timing_check=False,
show_plots=False,
te=0.073,
fov_xy=float(phantom.size.numpy()[0]),
n_readout=image_matrix_size[0],
n_phase_encoding=image_matrix_size[0],
partial_fourier_factor=0.75,
readout_type='symmetric',
)
mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))
signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)
fname_mrd = Path(tmp.name) / 'epi_fid_sym_ramp.mrd'
mr0.sig_to_mrd(fname_mrd, signal, sequence)
combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))
kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())
recon = DirectReconstruction(kdata, csm=None)
idata_fid_sym_ramp = recon(kdata)
Saving sequence file 'epi2d_fid_200fov_100px_sym_fid_2ro_rs_0p75pf_withnoise_withnav.seq' into folder '/home/runner/work/mrseq/mrseq/examples/output'.
>>>> Rust - compute_graph(...) >>>
Converting Python -> Rust: 0.000135429 s
Compute Graph
Computing Graph: 0.002329183 s
Analyze Graph
Analyzing Graph: 0.000001513 s
Converting Rust -> Python: 0.000008282 s
<<<< Rust <<<<
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/data/KData.py:166: UserWarning: No vendor information found. Assuming Siemens time stamp format.
warnings.warn('No vendor information found. Assuming Siemens time stamp format.', stacklevel=1)
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/operators/FourierOp.py:124: UserWarning: If both FFT and NUFFT dims are present, Cartesian FFT dims need to be aligned with the k-space dimension, i.e. kx along k0, ky along k1 and kz along k2. We are going to use NUFFT for all dimensions. Creating your own FourierOp with the desired combination of FFT and NUFFT dimensions might be more efficient.
return cls(
EPI FID with flyback readout and ramp sampling#
sequence, fname_seq = create_seq_fid(
system=sys_defaults,
test_report=False,
timing_check=False,
show_plots=False,
te=0.073,
fov_xy=float(phantom.size.numpy()[0]),
n_readout=image_matrix_size[0],
n_phase_encoding=image_matrix_size[0],
partial_fourier_factor=0.75,
readout_type='flyback',
)
mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))
signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)
fname_mrd = Path(tmp.name) / 'epi_fid_flyback_ramp.mrd'
mr0.sig_to_mrd(fname_mrd, signal, sequence)
combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))
kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())
recon = DirectReconstruction(kdata, csm=None)
idata_fid_flyback_ramp = recon(kdata)
Saving sequence file 'epi2d_fid_200fov_100px_flyb_fid_2ro_rs_0p75pf_withnoise_withnav.seq' into folder '/home/runner/work/mrseq/mrseq/examples/output'.
>>>> Rust - compute_graph(...) >>>
Converting Python -> Rust: 0.000160777 s
Compute Graph
Computing Graph: 0.002320018 s
Analyze Graph
Analyzing Graph: 0.000001302 s
Converting Rust -> Python: 0.000009785 s
<<<< Rust <<<<
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/data/KData.py:166: UserWarning: No vendor information found. Assuming Siemens time stamp format.
warnings.warn('No vendor information found. Assuming Siemens time stamp format.', stacklevel=1)
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/operators/FourierOp.py:104: UserWarning: K-space points lie outside of the encoding_matrix and will be ignored. Increase the encoding_matrix to include these points.
self._cart_sampling_op: CartesianSamplingOp | None = CartesianSamplingOp(
EPI SE with symmetric readout and ramp sampling#
sequence, fname_seq = create_seq_se(
system=sys_defaults,
test_report=False,
timing_check=False,
show_plots=False,
n_slices=1,
fov_xy=float(phantom.size.numpy()[0]),
n_readout=image_matrix_size[0],
n_phase_encoding=image_matrix_size[0],
partial_fourier_factor=0.75,
readout_type='symmetric',
)
mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))
signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-6)
fname_mrd = Path(tmp.name) / 'epi_se_sym_ramp.mrd'
mr0.sig_to_mrd(fname_mrd, signal, sequence)
combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))
kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())
recon = DirectReconstruction(kdata, csm=None)
idata_se_sym_ramp = recon(kdata)
Saving sequence file 'epi2d_se_200fov_100px_sym_se_2ro_rs_0p75pf_withnoise_withnav.seq' into folder '/home/runner/work/mrseq/mrseq/examples/output'.
>>>> Rust - compute_graph(...) >>>
Converting Python -> Rust: 0.000090694 s
Compute Graph
Computing Graph: 0.002220831 s
Analyze Graph
Analyzing Graph: 0.000001453 s
Converting Rust -> Python: 0.000009705 s
<<<< Rust <<<<
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/data/KData.py:166: UserWarning: No vendor information found. Assuming Siemens time stamp format.
warnings.warn('No vendor information found. Assuming Siemens time stamp format.', stacklevel=1)
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/mrpro/operators/FourierOp.py:124: UserWarning: If both FFT and NUFFT dims are present, Cartesian FFT dims need to be aligned with the k-space dimension, i.e. kx along k0, ky along k1 and kz along k2. We are going to use NUFFT for all dimensions. Creating your own FourierOp with the desired combination of FFT and NUFFT dimensions might be more efficient.
return cls(
fig, ax = plt.subplots(2, 4, figsize=(16, 8))
for cax in ax.flatten():
cax.set_xticks([])
cax.set_yticks([])
for idx, (idata, label) in enumerate(
zip(
[idata_fid_sym_no_ramp, idata_fid_sym_ramp, idata_fid_flyback_ramp, idata_se_sym_ramp],
['symmetric no ramp', 'symmetric ramp', 'flyback ramp', 'se symmetric ramp'],
strict=True,
)
):
idat = idata.data.abs().squeeze().numpy()
idat /= np.sort(idat.flatten())[int(idat.size * 0.99)]
if idx == 0:
idat_ref = idat
relative_error = 0.0
relative_error += np.sum(np.abs(idat - idat_ref)) / np.sum(np.abs(idat_ref))
ax[0, idx].imshow(idat, vmin=0, vmax=1, cmap='gray')
ax[0, idx].set_title(label)
ax[1, idx].imshow(idat - idat_ref, vmin=-1, vmax=1, cmap='bwr')
relative_error /= 3
print(f'Relative error {relative_error}')
assert relative_error < 0.20
Relative error 0.18769095838069916
