T1rho Mapping#

T1rho mapping using a spin echo sequence where a single readout line is acquired after a spin-lock preparation pulse. A long repetition time (TR) is used to ensure full signal recovery before the next k-space line is acquired.

Imports#

import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
import MRzeroCore as mr0
import numpy as np
import torch
from cmap import Colormap
from einops import rearrange
from mrpro.algorithms.reconstruction import DirectReconstruction
from mrpro.data import KData
from mrpro.data import SpatialDimension
from mrpro.data.traj_calculators import KTrajectoryCartesian
from mrpro.operators import DictionaryMatchOp
from mrpro.operators.models import MonoExponentialDecay

from mrseq.scripts.t1rho_se_single_line import main as create_seq
from mrseq.utils import sys_defaults

Settings#

We are going to use a numerical phantom with a matrix size of 128 x 128.

image_matrix_size = [128, 128]

tmp = tempfile.TemporaryDirectory()
fname_mrd = Path(tmp.name) / 't1rho.mrd'

Create the digital phantom#

We use the standard Brainweb phantom from MRzero, but we set the B0- and B1-field to be constant everywhere.

phantom = mr0.util.load_phantom(image_matrix_size)
phantom.B1[:] = 1.0
phantom.B0[:] = 0.0

Create the Spin-lock SE sequence#

To create the Spin-lock SE sequence, we use the previously imported t1rho_se_single_line script.

sequence, fname_seq = create_seq(
    system=sys_defaults,
    test_report=False,
    timing_check=False,
    fov_xy=float(phantom.size.numpy()[0]),
    n_readout=image_matrix_size[0],
    n_phase_encoding=image_matrix_size[1],
)

Minimum TE: 7.465 ms

Saving sequence file 't1rho_se_single_line_200fov_128nx_128ny_3TIs.seq' into folder '/home/runner/work/mrseq/mrseq/examples/output'.
_images/19db5f4174b377f6651e1f1ae2c3beb79ba079264db504ca31538f67c8b522ba.png _images/a5d38be7f1b27529f8cccf7c600352b841b5160b7508326c5b4d368c25d20303.png

Simulate the sequence#

Now, we pass the sequence and the phantom to the MRzero simulation and save the simulated signal as an (ISMR)MRD file.

mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))
signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-5)
mr0.sig_to_mrd(fname_mrd, signal, sequence)

>>>> Rust - compute_graph(...) >>>

Converting Python -> Rust: 0.008651544 s
Compute Graph
Computing Graph: 0.16550848 s
Analyze Graph
Analyzing Graph: 0.00685703 s
Converting Rust -> Python: 0.07146446 s
<<<< Rust <<<<

Reconstruct the images at different spin-lock times#

We use MRpro for the image reconstruction.

kdata = KData.from_file(fname_mrd, trajectory=KTrajectoryCartesian())
kdata.header.encoding_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0])
kdata.header.recon_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0])
recon = DirectReconstruction(kdata, csm=None)
idata = recon(kdata)

/opt/hostedtoolcache/Python/3.12.12/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)

We can now plot the images at different inversion times.

spin_lock_times = sequence.get_definition('TSL')
idat = idata.data.abs().numpy().squeeze()
fig, ax = plt.subplots(1, idat.shape[0], figsize=(5 * idata.shape[0], 5))
for i in range(idat.shape[0]):
    ax[i].imshow(idat[i, :, :], cmap='gray')
    ax[i].set_title(f'Spin-lock time = {int(spin_lock_times[i] * 1000)} ms')
    ax[i].set_xticks([])
    ax[i].set_yticks([])
_images/587c434a1b4935607933e0d9c81a4a3e1d00e9e944137156f4207528f718843d.png

Estimate the T1rho maps#

We use a dictionary matching approach to estimate the T1rho maps. This phantom does not have a specific T1rho model and therefore T1rho is assumed to be equal to T2. However, T2 for CSF is too long to be estimated properly, so we exclude it from the comparison.

dictionary = DictionaryMatchOp(
    MonoExponentialDecay(decay_time=sequence.get_definition('TSL')), index_of_scaling_parameter=0
)
dictionary.append(torch.tensor(1.0), torch.linspace(0.01, 0.8, 1000)[None, :])
m0_match, t1rho_match = dictionary(idata.data[:, 0, 0])

t2_input = np.roll(rearrange(phantom.T2.numpy().squeeze()[::-1, ::-1], 'x y -> y x'), shift=(1, 1), axis=(0, 1))
obj_mask = np.zeros_like(t2_input)
obj_mask[(t2_input > 0) & (t2_input < 0.1)] = 1
t2_input *= obj_mask
t1rho_measured = t1rho_match.numpy().squeeze() * obj_mask

fig, ax = plt.subplots(1, 3, figsize=(15, 3))
for cax in ax:
    cax.set_xticks([])
    cax.set_yticks([])

im = ax[0].imshow(t2_input, vmin=0, vmax=0.3, cmap=Colormap('navia').to_mpl())
fig.colorbar(im, ax=ax[0], label='Input T2 (s)')

im = ax[1].imshow(t1rho_measured, vmin=0, vmax=0.3, cmap=Colormap('navia').to_mpl())
fig.colorbar(im, ax=ax[1], label='Measured T2 (s)')

im = ax[2].imshow(t1rho_measured - t2_input, vmin=-0.03, vmax=0.03, cmap='bwr')
fig.colorbar(im, ax=ax[2], label='Difference T2 (s)')

relative_error = np.sum(np.abs(t2_input - t1rho_measured)) / np.sum(np.abs(t2_input))
print(f'Relative error {relative_error}')
assert relative_error < 0.2

Relative error 0.16339872777462006
_images/a038d9915d604896661f9815a13cecc69cdc2fc571d6ea8a0b91b9790b5cf80d.png
Show Source