Embedded Quantitative MRI T_1ρ Mapping Using Non-Linear Primal-Dual Proximal Splitting

Abstract

Quantitative MRI (qMRI) methods allow reducing the subjectivity of clinical MRI by providing numerical values on which diagnostic assessment or predictions of tissue properties can be based. However, qMRI measurements typically take more time than anatomical imaging due to requiring multiple measurements with varying contrasts for, e.g., relaxation time mapping. To reduce the scanning time, undersampled data may be combined with compressed sensing (CS) reconstruction techniques. Typical CS reconstructions first reconstruct a complex-valued set of images corresponding to the varying contrasts, followed by a non-linear signal model fit to obtain the parameter maps. We propose a direct, embedded reconstruction method for T${}_{1\rho }$ mapping. The proposed method capitalizes on a known signal model to directly reconstruct the desired parameter map using a non-linear optimization model. The proposed reconstruction method also allows directly regularizing the parameter map of interest and greatly reduces the number of unknowns in the reconstruction, which are key factors in the performance of the reconstruction method. We test the proposed model using simulated radially sampled data from a 2D phantom and 2D cartesian ex vivo measurements of a mouse kidney specimen. We compare the embedded reconstruction model to two CS reconstruction models and in the cartesian test case also the direct inverse fast Fourier transform. The T${}_{1\rho }$ RMSE of the embedded reconstructions was reduced by 37–76% compared to the CS reconstructions when using undersampled simulated data with the reduction growing with larger acceleration factors. The proposed, embedded model outperformed the reference methods on the experimental test case as well, especially providing robustness with higher acceleration factors.

1. Introduction

Magnetic resonance imaging (MRI) is one of the most important tools for the clinical diagnosis of various diseases due to its excellent and versatile soft tissue contrast. Clinical MRI is based on expert interpretation of anatomical images of varying contrasts and thus tends to retain a level of subjectivity. Quantitative MRI (qMRI) methods, such as measurements of different relaxation times, allow reducing the subjectivity by providing numerical values on which diagnostic assessment or predictions of tissue properties can be based on.

However, such quantitative MRI measurements necessarily take more time than standard anatomical imaging. For example, in T${}_{1\rho }$ mapping [1,2], typically, 5–7 sets of measurements with varying spin lock times are collected to estimate the T${}_{1\rho }$ map. Such measurements will thus take 5–7 times longer than acquiring similar anatomical images, often approaching 10 min for a stack of quantitative 2D images.

T${}_{1\rho }$ imaging is based on tilting the magnetization into the $xy$-plane and then locking the magnetization with a spin-lock pulse of a certain amplitude and duration. Quantitative mapping, i.e., the measurement of the T${}_{1\rho }$ relaxation time constant, is realized by repeating the T${}_{1\rho }$ preparation with several different durations of the spin-lock pulse and collecting the full MR image for each of these preparations. The T${}_{1\rho }$ MRI contrast is particularly sensitive to molecular processes occurring at the frequency (${\omega }_1$) of the spin-lock pulse corresponding to the amplitude of the pulse: ${\omega }_1=\gamma B_1$, where $\gamma$ is the gyromagnetic ratio, which ties the magnetic field strength (of the radio frequency (RF) pulse) $B_1$ to its resonance frequency. Generally, spin-lock pulses operate at and are limited to frequencies that correspond to slow molecular processes that are often both biologically important and altered in disease-related changes. The T${}_{1\rho }$ relaxation time has been reported as a promising biomarker for numerous tissues and diseases, such as different disorders of the brain [3,4], cardiomyopathy [5], liver fibrosis [6], musculoskeletal disorders [2,7,8] and many others. For a broader overview of T${}_{1\rho }$ relaxation and its applications, the reader is referred to the reviews by Gilani and Sepponen [1], Wang and Regatte [7] and Borthakur et al. [2].

Staying still in the scanner for extended periods of time can prove to be challenging, for example, for pediatric patients. The excessively long data acquisition times are also operationally impractical because they lead to a small number of studies that can be performed daily with a single MRI device. Quantitative MRI and T${}_{1\rho }$ imaging in particular can thus greatly benefit from using undersampled measurements, which are a natural and efficient way to reduce the scanning time for a single qMRI experiment. When using undersampled data, conventional MR image reconstruction methods, such as regridding [9], may lead to insufficient reconstruction quality. The usage of compressed sensing (CS) [10,11] methods, where an iterative reconstruction method is used together with a sparsifying transform of the image, has proven highly successful with undersampled MRI data [12].

Usage of CS methods for T${}_{1\rho }$ imaging have been previously studied, for example, in [13,14,15,16]. In [13], the authors used principal component analysis and dictionary learning in the first in vivo application of CS to T${}_{1\rho }$ reconstruction. In [14], the authors used spatial total variation (TV) together with Autocalibrating Reconstruction for Cartesian sampling (ARC) to accelerate the measurements. In [15], the authors compared 12 different sparsifying transforms in 3D T${}_{1\rho }$ mapping. The regularization model combining spatial TV with second-order contrast TV was found to perform the best, with satisfactory results with an acceleration factor (AF, i.e., the number of datapoints in full data divided by the number of data used in the reconstruction) up to 10 when using cartesian 3D sampling together with parallel imaging. In [16], both cartesian and radial data were reconstructed using various different regularization methods. The authors reached acceptable accuracy with AF up to 4 for the cartesian data, whereas with the radial data, the accuracy was acceptable with AF up to 10.

When using CS for T${}_{1\rho }$ mapping, the image series with varying spin-lock durations $T_{\mathrm{SL}}$ is first reconstructed, followed by a pixel-by-pixel non-linear least squares fit of a monoexponential (or a biexponential) signal model to the reconstructed image intensity data to obtain the desired T${}_{1\rho }$ relaxation time map. Since the exponential signal model combining the T${}_{1\rho }$ and varying $T_{\mathrm{SL}}$ is well known, a direct, embedded model may also be used to reconstruct the desired T${}_{1\rho }$ map directly from the k-space measurement data without the intermediate step of reconstructing the separate intensity maps corresponding to different $T_{\mathrm{SL}}$. Figure 1 shows a schematic of the CS T${}_{1\rho }$ mapping method as well as the direct, embedded model.

The direct one-step reconstruction utilizing the embedded model has clear advantages over the sequential two-step reconstruction model. First, it reduces the number of unknowns in the reconstruction problem significantly; for example, for measurements with seven spin-lock times, the number of unknowns may be reduced from $14N$ (one complex image for each contrast) to just $3N$ (T${}_{1\rho }$, S${}_0$ and a single phase map), where N is the number of pixels or voxels in a single image. Secondly, it allows the regularization of the parameter map of interest, i.e., the T${}_{1\rho }$ parameter map in the case of T${}_{1\rho }$ mapping instead of posing regularization on the complex-valued images corresponding to different contrasts in the intermediate step. Thirdly, since the signal model is embedded in the reconstruction, there is no need to decide what type of a contrast regularization model fits the data best.

A disadvantage of the embedded model is that it transforms the MRI inversion into a non-linear problem, which is not necessarily convex and thus requires proper initialization. The resulting non-linear and possibly non-convex optimization problem can, however, be solved conveniently with, for example, the non-linear primal-dual proximal splitting algorithm [17].

Alternatively, various deep learning approaches have also been proposed for different aspects of quantitative MRI. For example, in [18], the authors propose the use of deep learning neural networks to reduce the number of contrasts required for an accurate model fit in myocardial T${}_1$ mapping. Additionally, a model-guided self-supervised deep learning MRI reconstruction framework for direct T${}_1$ and T${}_2$ parameter mapping has been proposed [19]. For an overview of the usage of deep learning in MR relaxometry, see [20].

In this work, we propose an embedded parameterization model to directly reconstruct the T${}_{1\rho }$, S${}_0$, and phase maps from the k-space measurement data and use the non-linear primal-dual proximal splitting algorithm to solve the problem. The proposed model is tested with 2D simulated radial phantom data and 2D cartesian ex vivo mouse kidney data. The proposed embedded model is compared with two CS models: one with spatial TV and TV over the $T_{\mathrm{SL}}$ contrasts, which, we believe, is generally the most commonly used CS model in MRI, and a second CS model with spatial TV and second-order contrast TV, which in [15] was found to perform the best out of 12 different CS models for T${}_{1\rho }$ mapping. The first CS model is labeled “CS S1+C1”, and the second CS model is labeled “CS S1C2” throughout the paper. The models are named slightly different since in the first model, the spatial and contrast TV components are separate with two different regularization parameters, and in the second model, the spatial TV and the second-order contrast TV are under the same root with a single regularization parameter. In the cartesian test case, results from a direct inverse fast Fourier transform (iFFT) model are also shown as a reference. Reconstructions from both the CS models and the iFFT model are followed by the mono-exponential pixel-by-pixel T${}_{1\rho }$ fit.

2. Reconstruction Methods

2.1. Embedded T${}_{1\rho }$ Model

The signal model in T${}_{1\rho }$ mapping is

$$S_{c,k}=S_{0,k}exp(-T_{{\mathrm{SL}}_c}/{\mathrm{T}}_{1{\rho }_k}),$$

(1)

where $S_{c,k}$ is the signal intensity with spin-lock time $T_{{\mathrm{SL}}_c}$, where c denotes the contrast index and k denotes the pixel index, and $S_0$ is the proton density map, i.e., the signal intensity when $T_{\mathrm{SL}}=0$. For the recovery of the T${}_{1\rho }$ map, k-space measurement data are collected by scanning the target with multiple spin-lock times $T_{{\mathrm{SL}}_c}$. The measurement model mapping the $S_0$, T${}_{1\rho }$, and the phase map $\theta$ to the k-space measurements then reads

$$m=K({\mathrm{S}}_0,{\mathrm{T}}_{1\rho },\theta )+e,$$

(2)

where the vectors S${}_0$, T${}_{1\rho }$, and $\theta \in {\mathbb{R}}^N$ are the parameter maps to be reconstructed, and the complex measurement vector $m\in {\mathbb{C}}^{CM}$ is composed of k-space data with C spin-lock times, each consisting of M measurements. Further, we denote the complex-valued measurement noise by $e\in {\mathbb{C}}^{CM}$ and the non-linear forward model by $K:{\mathbb{R}}^{3N}\to {\mathbb{C}}^{CM}$.

The non-linear forward model can be further decomposed to

$$K({\mathrm{S}}_0,{\mathrm{T}}_{1\rho },\theta )=AB\left(D({\mathrm{S}}_0,{\mathrm{T}}_{1\rho },\theta )\right),$$

(3)

where A is the block-diagonal matrix containing the Fourier transform operations. In the case of cartesian measurements, the blocks of A read $A^c=U^c\mathcal{F}$, where $U^c$ is the undersampling pattern used with the measurements with contrast index c, and $\mathcal{F}$ is the Fourier transform. In the case of non-cartesian measurements, we approximate the forward model using the non-uniform fast Fourier transform (NUFFT [21]), i.e., $A^c=P^c\mathcal{F}{L}^c$, where $P^c$ is an interpolation and sampling matrix, and $L^c$ is a scaling matrix. Furthermore, D maps the S${}_0$ and T${}_{1\rho }$ parameter maps to magnitude images as

$$\begin{array}{c}\hfill D({\mathrm{S}}_0,{\mathrm{T}}_{1\rho },\theta )=\left[\begin{array}{c}{\mathrm{S}}_0\odot exp(-T_{{\mathrm{SL}}_1}/{\mathrm{T}}_{1\rho })\\ ⋮\\ {\mathrm{S}}_0\odot exp(-T_{{\mathrm{SL}}_C}/{\mathrm{T}}_{1\rho })\\ \theta \end{array}\right]:=\left[\begin{array}{c}{r}_1\\ ⋮\\ r_C\\ \theta \end{array}\right],\end{array}$$

(4)

where ⊙ is the Hadamard product, i.e., elementwise multiplication, and the exponentiation and the division of the scalars $T_{{\mathrm{SL}}_i}$ by the vector T${}_{1\rho }$ are to be interpreted as elementwise operations. Moreover, B maps the magnitude and phase components of the images to real and complex components and can be expressed as

$$B(r_1,...,r_C,\theta )=\left[\begin{array}{c}{r}_1\odot cos\theta \\ r_1\odot sin\theta \\ ⋮\\ r_C\odot cos\theta \\ r_C\odot sin\theta \end{array}\right].$$

(5)

Here too, the sin and cos are to be interpreted as elementwise operations. Note that if the phase maps vary between contrasts, the model can be easily modified to reconstruct separate phase maps for all contrasts instead of reconstructing only a single phase map. In a typical T${}_{1\rho }$ measurement, however, the contrast preparation is usually well-separated from the imaging segment, and thus, the phase can be expected to be the same between the otherwise identical image acquisition segments.

In the embedded reconstruction, we use total variation regularization for the S${}_0$ and T${}_{1\rho }$ maps and L2-norm regularization for the spatial gradient of the phase map. TV regularization has been shown to be one of the best performing approaches with CS in T${}_{1\rho }$ mapping [15], and thus, for a fair comparison, it was chosen as regularization for the S${}_0$ and T${}_{1\rho }$ maps. Additionally, gradient L2 regularization was used for the phase map since the phase maps are most often smooth. We also limit the S${}_0$ and T${}_{1\rho }$ parameter maps above a small positive value. With these, the minimization problem reads

$$\begin{array}{cc}\underset{{\mathrm{S}}_0,{\mathrm{T}}_{1\rho },\theta }{min}\left|\right|K({\mathrm{S}}_0,{\mathrm{T}}_{1\rho },\theta )-m{\left|\right|}_2^2+{\alpha }_1{\mathrm{TV}}_{\mathrm{S}}\left({\mathrm{S}}_0\right)+\hfill & {\alpha }_2{\mathrm{TV}}_{\mathrm{S}}\left({\mathrm{T}}_{1\rho }\right)\hfill \\ & \hspace{1em}\hspace{1em}+{\alpha }_3\left|\right|{\nabla }_{\mathrm{S}}\theta {\left|\right|}_2^2+{\delta }_{\mathrm{a}}1\left({\mathrm{S}}_0\right)+{\delta }_{\mathrm{a}}2\left({\mathrm{T}}_{1\rho }\right)\hfill \end{array}$$

(6)

where ${\mathrm{TV}}_{\mathrm{S}}$ denotes spatial total variation, ${\nabla }_{\mathrm{S}}$ is the spatial discrete difference operator, and ${\delta }_{\mathrm{a}}\mathrm{i}$ are step functions with an infinite value below the parameter $a_i$ and 0 above or equal to the parameter. Further, ${\alpha }_1$, ${\alpha }_2$, and ${\alpha }_3$ are the regularization parameters for the S${}_0$, T${}_{1\rho }$, and phase maps, respectively, and $a_1$ and $a_2$ are the small positive constraints on the S${}_0$ and T${}_{1\rho }$ maps, respectively.

Solving the Embedded T${}_{1\rho }$ Reconstruction Problem

The non-linear, non-smooth optimization problem in Equation (6) is solved using the non-linear primal-dual proximal splitting algorithm proposed in [17], which is described in Algorithm 1 in its most general form. Here, the non-linear mapping $H:{\mathbb{R}}^{3N}\to {\mathbb{C}}^{CM+6N}$ contains the non-linear forward model K and the discrete difference matrices. The algorithm applied to the embedded T${}_{1\rho }$ reconstruction is described in more detail in the Appendix A.

Algorithm 1 Non-linear primal-dual proximal splitting presented in [17] (Algorithm 2.1)

$\mathrm{Choose}\omega \ge 0,\mathrm{and}\tau ,\sigma$ $\mathrm{s}.\mathrm{t}.\tau \sigma \left({sup}_{k=1,...,i}{\Vert \nabla H\left(x^k\right)\Vert }^2\right)<1.$ while Not reached stopping criterion do     $x^{i+1}:={(I+\tau \partial G)}^{-1}(x^i-\tau {[\nabla H\left(x^i\right)]}^{*}{y}^i)$     ${\overline{x}}^{i+1}:=x^{i+1}+\omega (x^{i+1}-x^i)$     $y^{i+1}:={(I+\sigma \partial F^{*})}^{-1}(y^i+\sigma H\left({\overline{x}}^{i+1}\right))$ end while

In our implementation, $x={({\mathrm{S}}_0^{\mathrm{T}},{\mathrm{T}}_{1\rho }^{\mathrm{T}},{\theta }^{\mathrm{T}})}^{\mathrm{T}}$, and we initialize the S${}_0$ and phase parts of $x^0$ using iFFT or adjoint of NUFFT of the $T_{\mathrm{SL}}=0$ measurements. T${}_{1\rho }$ was initialized to a constant value of 20, and the dual variable y was initialized to 0. Initializing the S${}_0$ map with a constant value instead of the iFFT or adjoint of NUFFT, the algorithm generally fails to converge to feasible solutions, whereas initializing the T${}_{1\rho }$ map with different feasible values, the algorithm converges to nearly the same solution with differences mainly in convergence speed.

In addition, we use varying primal step sizes for the different blocks of the embedded reconstruction, i.e., different ${\tau }_i$ parameters for the S${}_0$, T${}_{1\rho }$, and phase updates [22]. This essentially replaces the scalar step length parameter $\tau$ in Algorithm 1 with the diagonal matrix

$$\mathrm{T}=\left[\begin{array}{ccc}{\tau }_1{I}_N& 0& 0\\ 0& {\tau }_2{I}_N& 0\\ 0& 0& {\tau }_3{I}_N\end{array}\right].$$

(7)

The step parameters ${\tau }_1$, ${\tau }_2$, and ${\tau }_3$ are derived from the norm of the corresponding block of the matrix $\nabla H$. Here, however, we only use the non-linear part K of H to estimate the step lengths, as the linear part of H has only a minor impact on the norm of $\nabla H$. We set the parameter $\sigma$ to $\sigma =1/max\left({\tau }_i\right)$ and use $\omega =1$ for the relaxation parameter.

Since the block-diagonal matrix A is linear and can be normalized to 1, we have $\Vert \nabla K\Vert =\Vert J_B{J}_D\Vert$. Furthermore, the product of the Jacobians writes

$$\begin{array}{c}\hfill J_B{J}_D=\left[\begin{array}{ccc}\mathrm{diag}(cos\theta )E_1& \mathrm{diag}((T_{\mathrm{SL}}1/({\mathrm{T}}_{1\rho }\odot {\mathrm{T}}_{1\rho }))\odot cos\theta \odot r_1)& -\mathrm{diag}(sin\theta \odot r_1)\\ \mathrm{diag}(sin\theta )E_1& \mathrm{diag}((T_{\mathrm{SL}}1/({\mathrm{T}}_{1\rho }\odot {\mathrm{T}}_{1\rho }))\odot sin\theta \odot r_1)& \mathrm{diag}(cos\theta \odot r_1)\\ \mathrm{diag}(cos\theta )E_2& \mathrm{diag}((T_{\mathrm{SL}}2/({\mathrm{T}}_{1\rho }\odot {\mathrm{T}}_{1\rho }))\odot cos\theta \odot r_2)& -\mathrm{diag}(sin\theta \odot r_2)\\ \mathrm{diag}(sin\theta )E_2& \mathrm{diag}((T_{\mathrm{SL}}2/({\mathrm{T}}_{1\rho }\odot {\mathrm{T}}_{1\rho }))\odot sin\theta \odot r_2)& \mathrm{diag}(cos\theta \odot r_2)\\ ⋮& ⋮& ⋮\\ \mathrm{diag}(cos\theta )E_C& \mathrm{diag}((T_{\mathrm{SL}}\mathrm{C}/({\mathrm{T}}_{1\rho }\odot {\mathrm{T}}_{1\rho }))\odot cos\theta \odot r_C)& -\mathrm{diag}(sin\theta \odot r_C)\\ \mathrm{diag}(sin\theta )E_C& \mathrm{diag}((T_{\mathrm{SL}}\mathrm{C}/({\mathrm{T}}_{1\rho }\odot {\mathrm{T}}_{1\rho }))\odot sin\theta \odot r_C)& \mathrm{diag}(cos\theta \odot r_C)\end{array}\right],\end{array}$$

(8)

where $E_i=\mathrm{diag}(exp(-T_{\mathrm{SL}}\mathrm{i}/{\mathrm{T}}_{1\rho }))$ and $r_i={\mathrm{S}}_0\odot exp(-T_{\mathrm{SL}}\mathrm{i}/{\mathrm{T}}_{1\rho })$. Now, since the matrix $J_B{J}_D$ consists of only diagonal blocks, and the index of the maximum value is the same for all $E_i$, it is straightforward to estimate the ${\tau }_i$ from the norms of the maximum values of the column-blocks of Equation (8) yielding

$$\begin{array}{cc}\hfill {\tau }_1=& \sqrt{\sum _{i=1}^C\Vert {exp(-T_{\mathrm{SL}}\mathrm{i}/{\mathrm{T}}_{1\rho })\Vert }_{\infty }^2}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\tau }_2=& \sqrt{\sum _{i=1}^C\Vert {r_i\odot (T_{\mathrm{SL}}\mathrm{i}/({\mathrm{T}}_{1\rho }\odot {\mathrm{T}}_{1\rho }))\Vert }_{\infty }^2}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\tau }_3=& \sqrt{\sum _{i=1}^C\Vert {{\mathrm{S}}_0\odot exp(-T_{\mathrm{SL}}\mathrm{i}/{\mathrm{T}}_{1\rho })\Vert }_{\infty }^2}.\hfill \end{array}$$

(11)

In addition, we calculate the norms in every iteration and update the used ${\tau }_i$ and $\sigma$ if the step is smaller than the previously used step.

In our experience, these step lengths may, however, prove to be too small, and in some cases, larger step lengths, especially for the T${}_{1\rho }$ update step, may be used to obtain faster convergence. In this work, we used a multiplier of 50 for the T${}_{1\rho }$ update step ${\tau }_2$ in the radial simulation. Note that the step length criterion of Algorithm 1 still holds with the multiplier since ${\tau }_2·\sigma$ remains small due to the selection of $\sigma$.

2.2. Compressed Sensing Reference Methods

We compare the embedded model to two CS models, which include a complex valued reconstruction of the images with different spin-lock times, followed by a pixel-by-pixel non-linear least squares fit of the monoexponential signal model to obtain the T${}_{1\rho }$ and S${}_0$ parameter maps. The first CS reconstruction model uses spatial total variation together with first-order total variation over the varying $T_{SL}$ contrasts (labeled CS S1+C1), and the second one uses spatial total variation together with second-order total variation over the varying $T_{SL}$ contrasts (labeled CS S1C2).

The measurement model for a single contrast image is

$$m^c=A^c{u}^c+e^c,$$

(12)

where the superscript c denotes the contrast index, $m^c\in {\mathbb{C}}^M$ is the k-space data vector for contrast index c, $u^c\in {\mathbb{C}}^N$ is the image vector, $e^c\in {\mathbb{C}}^M$ is the complex valued noise vector, and $A^c$ is the forward model, which depends on the measurement sequence and undersampling pattern and is described in more detail in Section 2.1.

With the measurement model of Equation (12), spatial total variation, and total variation over the contrasts, the CS minimization problem reads

$$u^{*}=\underset{u}{arg\; min}\Vert {Au-m\Vert }_2^2+\alpha {\mathrm{TV}}_{\mathrm{S}}\left(u\right)+\beta {\mathrm{TV}}_{\mathrm{C}}\left(u\right),$$

(13)

where A is a block-diagonal matrix containing the forward transforms $A^c$ corresponding to each image, $u\in {\mathbb{C}}^{NC}$ is all the images vectorized, such that C is the number of contrasts, and $m\in {\mathbb{C}}^{MC}$ is all the k-space measurements vectorized. Further, ${\mathrm{TV}}_{\mathrm{S}}$ denotes spatial total variation, ${\mathrm{TV}}_{\mathrm{C}}$ denotes total variation over contrasts, and $\alpha$ and $\beta$ are the regularization parameters of spatial and contrast TV, respectively.

The second CS minimization problem, which uses the single regularization parameter version of combined spatial TV and second-order contrast TV, reads

$$u^{*}=\underset{u}{arg\; min}\Vert {Au-m\Vert }_2^2+\alpha {\mathrm{TV}}_{\mathrm{SC}}\left(\mathrm{u}\right),$$

(14)

where

$${\mathrm{TV}}_{\mathrm{SC}}\left(u\right)=\sum _k\sqrt{{\left({\nabla }_{\mathrm{x}}u\right)}_k^2+{\left({\nabla }_{\mathrm{y}}u\right)}_k^2+{\left({\nabla }_{\mathrm{c}}^{2}u\right)}_k^2},$$

where ${\nabla }_{\mathrm{x}}$ and ${\nabla }_{\mathrm{y}}$ are the horizontal and vertical direction spatial discrete forward difference operators, respectively, ${\nabla }_{\mathrm{c}}^2$ is the second order contrast direction discrete difference operator, and k is an index that goes through all the pixels in the set of images.

Both of the minimization problems (Equations (13) and (14)) are solved using the popular primal-dual proximal splitting algorithm of Chambolle and Pock [23].

Finally, in the CS models (and the iFFT model), we fit the mono-exponential T${}_{1\rho }$ signal equation

$$[{\mathrm{T}}_{1\rho ,k}^{*},{\mathrm{S}}_{0,k}^{*}]=\underset{{\mathrm{T}}_{1\rho },{\mathrm{S}}_0}{arg\; min}\Vert {|u_k|-{\mathrm{S}}_{0}exp(-T_{\mathrm{SL}}/{\mathrm{T}}_{1\rho })\Vert }^2$$

(15)

pixel by pixel to the reconstructed intensity images obtained by solving either Equation (13) or Equation (14). Here, $|u_k|=|u_k^1|,...,|u_k^C|$ is the vector of reconstruction intensity values at pixel location k with $T_{\mathrm{SL}}$ contrasts 1 to C, and similarly, $T_{\mathrm{SL}}$ is the vector of $T_{\mathrm{SL}}$ values of contrasts 1 to C. Note that the final S${}_0$ estimate is obtained from the mono-exponential model fit instead of taking the intensity values from the reconstructions with $T_{\mathrm{SL}}=0$.

3. Materials and Methods

3.1. Simulated Golden Angle Radial Data

The simulation of the radial measurement data was based on the Shepp–Logan phantom in dimensions $128\times 128$, which was zero-filled to dimensions $192\times 192$. The T${}_{1\rho }$ values of the target were set to between 20 and 120. The intensity with $T_{\mathrm{SL}}=0$ was set to a maximum of 1, and the phase of the target was set $2\pi x/192$, where x is the horizontal coordinate of the pixel. The images of the simulated T${}_{1\rho }$, S${}_0$, and phase maps are shown in Figure 2. To generate the varying $T_{\mathrm{SL}}$ measurements, spin lock times of 0, 4, 8, 16, 32, 64, and 128 ms were used. For each $T_{\mathrm{SL}}$, 302 (i.e., $\sim 192·\pi /2$) golden angle [24] spokes were generated. This corresponds to full sampling for equispaced radial spokes with image dimensions $192\times 192$ in the sense that the distance between spokes at their outermost points satisfies the Nyquist criterion [25]. Finally, complex Gaussian noise at 5 % of the mean of the absolute values of the full noiseless simulation was added to the simulated measurements.

3.2. Cartesian Data from Ex Vivo Mouse Kidney

Experimental ex vivo data from a mouse kidney was acquired from a separate study. The data were collected in compliance with ethical permits (ESAVI/270/04.10.07/2017) at 9.4 T using a 19 mm quadrature RF volume transceiver (RAPID Biomedical GmbH, Rimpar, Germany) and VnmrJ3.1 Varian/Agilent DirectDrive console. T${}_{1\rho }$ relaxation data were collected using a refocused T${}_{1\rho }$ preparation scheme [26] with a spin-lock frequency of 500 Hz and $T_{\mathrm{SL}}$ = 0, 8, 16, 32, 64, and 128 ms. The T${}_{1\rho }$-prepared data, i.e., T${}_{1\rho }$-weighted images, were collected using a fast spin echo sequence with a repetition time of 5 s, effective echo time of 5.5 ms, echo train length of 8, slice thickness of 1mm, field-of-view of 17 × 17 mm and acquisition matrix of 192 × 192. Eventually, only spin-lock times up to 64 ms were used in the reconstruction as the signal intensity of the longest spin-lock time was close to the noise level and had minimal or no effect on the reconstruction.

3.3. Reconstruction Specifics

The radial data from the 2D phantom were reconstructed with the embedded model and the two CS models with acceleration factors of 1, 5, 10, 20, 30, 50, and 101 (rounded to the nearest integer). In T${}_{1\rho }$ imaging, the images measured with varying spin-lock times are expected to have high redundancy in the sense that the images are expected to be structurally similar with decreasing intensity as $T_{\mathrm{SL}}$ increases, making complementary k-space sampling warranted. In complementary k-space sampling, the subsampling with any measured contrast is different from the others, meaning that each sampling adds to spatial information gained at other contrasts. The golden angle radial sampling is especially well suited for this as the measurements are inherently complementary (i.e., each new spoke has a different path in the k-space compared to the previous ones), and each measured spoke traverses through the central (low-frequency) part of the k-space, which contains significant information on the intensity level of the images. Thus, we sampled the golden angle data such that, for example, with an acceleration factor of 0, the first contrast used the first 30 spokes out of the 302 total, the second contrast used spokes 31 through 60 and so on to achieve complementary k-space sampling. Examples of the radial sampling pattern for an acceleration factor of 20 and the cartesian sampling pattern for acceleration factor of 5 are shown in Figure 3.

In the embedded model, the phase regularization parameter was set to a constant value at 0.01, and the other regularization parameters were varied over a wide range. In the CS models, the regularization parameters were also varied over a wide range to find the best parameters. The reconstructions shown use the regularization parameters that yielded the smallest T${}_{1\rho }$ RMSE with respect to the ground truth phantom.

The NUFFT operator used in the radial data reconstructions was implemented using the Michigan Image Reconstruction Toolbox (MIRT) [27]. The interpolator used was the minmax:kb interpolator with a neighbourhood size of 4 and scaling factor of 2.

The cartesian ex vivo mouse kidney data were reconstructed with the embedded, the iFFT, and the two CS methods with acceleration factors of 2, 3, 4, and 5 (rounded to the nearest integer). Undersampling was conducted by taking a number of full k-space rows corresponding to the desired acceleration factor since cartesian data collection in MRI scanners is carried out line by line. For the undersampled reconstructions, 1/4 of the total sampled k-space rows were taken from around the center to include zero frequency and enough low-frequency data in all contrasts. Half of the rest 3/4 were taken from the top part and the other half from the bottom part. To achieve complementary sampling, the rows from the top and bottom parts were selected such that all rows were first selected once in random order before continuing to sample from the full set of rows again (Figure 3).

In the ex vivo test case, too, the phase regularization parameter of the embedded model was set to a constant level, which was 0.0001, and the other parameters of the embedded and both CS models were varied over a wide range to find the optimal T${}_{1\rho }$ estimate. The embedded model reconstructions were compared to the embedded reconstruction with full data, and likewise, the CS and iFFT model reconstructions were compared to the corresponding reconstructions with full data as the true T${}_{1\rho }$ map is not available. Thus, the RMSEs reflect each model’s relative tolerance for undersampling compared to the situation where fully sampled data are available for the particular reconstruction model.

4. Results

4.1. Simulated Golden Angle Radial Data

With the radial simulated phantom data, all the methods produce reconstructions with similar RMSEs when using full data (acceleration factor 1). With undersampled data, the embedded model outperforms both the CS models as measured by RMSE of both the T${}_{1\rho }$ (Figure 4) and S${}_0$ (Figure 5) maps with all acceleration factors and the improvement increases with larger acceleration factors.

The T${}_{1\rho }$ maps computed using the CS models are also visibly noisier as the model does not allow direct regularization of the T${}_{1\rho }$ map (Figure 4). With an acceleration factor of 101, reconstructions of both CS models start to break down, whereas the embedded model reconstruction still reconstructs the target reasonably well, with RMSE values below those of the the CS models at an acceleration factor of 20–30 (Figure 4, Figure 5 and Figure 6).

4.2. Cartesian Data from Ex Vivo Mouse Kidney

In the cartesian ex vivo test case, the performance of the embedded and CS models in their relative tolerance for undersampling is similar with an acceleration factor of 2, and both CS models perform slightly worse than the embedded model with an acceleration factor of 3 (Figure 7, Figure 8 and Figure 9). With an acceleration factor of 4, the performance of the CS models is already clearly worse than the performance of the embedded model, and while both of the CS models fail in the reconstruction with an acceleration factor of 5, the embedded model still produces similar tolerance for undersampling as with the smaller acceleration factors. The undersampled iFFT reconstructions shown for reference perform worse than the CS or the embedded model reconstructions with all the acceleration factors.

5. Discussion

In this work, we proposed a non-linear, embedded T${}_{1\rho }$ model for direct quantitative T${}_{1\rho }$ reconstruction. The model is solved using the non-linear primal-dual proximal splitting algorithm [17]. We compared the embedded model reconstructions to two compressed sensing reconstructions followed by a mono-exponential T${}_{1\rho }$ fit in a radial simulated test case and a cartesian ex vivo test case. In the cartesian test case, we also show results from iFFT reconstructions followed by the T${}_{1\rho }$ fit.

In the simulated test case, where the RMSE metric with respect to the true target image is available, the embedded model outperformed both of the CS models with improvement increasing towards the higher acceleration factors. In the experimental test case with Cartesian ex vivo mouse kidney data, the RMSEs reflect the relative tolerance of the method with respect to the case where the fully sampled data were available for that particular method. In this case, the embedded model and the CS models had similar RMSEs for an acceleration factor of 2, and for higher acceleration factors, the embedded model clearly exhibited better tolerance for undersampling, indicating that the embedded model would allow the usage of higher acceleration factors than the CS models.

We believe the main factor for the better performance of the embedded model, especially with higher acceleration factors, is the reduction in the reconstructed parameters. In the simulation, in the standard CS approach, there are $14N$ unknowns, where N is the number of pixels in a single image, whereas in the embedded model, there are $3N$ unknowns. This is a reduction of 79% in the number of reconstructed unknowns. The same also holds true for the experimental case, where the reduction is 70%. Thus, when utilizing the embedded model, the problem is less undersampled—in the sense of the number of unknowns compared to the number of measurement points—than when using the CS models.

The two CS models perform quite similarly with the second-order contrast TV model CS S1C2 performing slightly better overall than the CS S1+C1 model in the simulated test case. The same observation can be made in the cartesian test case up to an acceleration factor of 4. In the Cartesian test case, the CS S1+C1 model has a smaller RMSE than CS S1C2 with an acceleration factor of 5, but in this case, both of the CS models failed to produce useful T${}_{1\rho }$ or S${}_0$ maps. From the practical point of view, the second-order contrast TV model with the implementation described in [15] is also more convenient than the CS S1+C1 model as it requires selecting only a single regularization parameter.

The embedded model is, however, slower to compute than the CS models. For example, our code implementation running on MATLAB (R2017b, The MathWorks, Inc., Natick, MA, USA) using an Intel Xeon E5-2630 CPU took 104 min for the embedded model and 26 min for the CS S1+C1 model with the radial simulation data with $\mathrm{AF}=5$. For the experimental cartesian data, the difference was bigger: for example, for $\mathrm{AF}=2$, the embedded model took 75 min to compute, while the CS S1+C1 model converged to stopping criterion in under a minute. The computation times could, however, be shortened, for example, by optimizing the code, running the code on a GPU, and also loosening the stopping criteria since we ran the iterations with rather strict criteria.

In the radial simulated test case, the embedded model reconstructs the target quite well even with an acceleration factor of 101, using only three spokes per $T_{SL}$ contrast, and 21 spokes in the whole reconstruction. In the cartesian test case, the acceleration factors that can be reached are much smaller. Even though the target used in the radial simulation is rather simple, it is evident that the radial sampling pattern, particularly with the golden angle sampling where k-space spokes are complementary and go through the center part of the k-space, allows much higher acceleration factors than a cartesian line-by-line sampling pattern. This is due to the undersampling artefacts in radial sampling (i.e., streaking) being more noise-like in the transform domain than the undersampling artefacts that arise in cartesian sampling [28,29]. This finding is aligned with the findings of [16].

Testing the proposed embedded model with radial experimental data, in vivo data, 3D data, and parallel imaging data are interesting future works, and our hypothesis is that similar results, where the embedded model outperforms the CS models, are to be expected. In addition, the embedded T${}_{1\rho }$ model could be tested with other regularizers, such as total generalized variation [30], which balances between minimizing the first- and second-order differences of the signal, making the results less piecewise constant, an issue for TV regularization, which is visible in the embedded reconstructions in, e.g., Figure 7. Other regularizers, which could alleviate the over-smoothness, include, for example, non-local means [31] or dictionary learning [32].

As the contrast manipulation scheme of the signal acquisition and the quantitative signal equation are the only major aspects that change between different qMRI contrasts, the proposed method can easily be adapted to fit other qMRI cases as well. Besides other qMRI methods, other aspects where embedded modelling could offer further benefits are T${}_{1\rho }$ dispersion imaging [33,34], where the data are acquired at multiple spin-locking amplitudes, and reducing RF energy deposition by asymmetric data reduction for the different spin-lock times (i.e., less data for long spin-lock pulses). More generally, shorter scan times may allow for higher spin-lock durations and/or higher amplitude pulses, as the specific absorption rate of RF energy can be minimized via acquiring less data for the most demanding pulses. Alternatively, multi-contrast embedded modelling could offer further avenues for data reduction.

6. Conclusions

In this work, we proposed an embedded T${}_{1\rho }$ reconstruction method, which directly reconstructs the T${}_{1\rho }$, S${}_0$, and phase maps from the measurement data. The reconstruction method also allows direct regularization of these parameter maps, and thus, a priori information about the parameter maps may be incorporated into the reconstruction. We also showed that the proposed method outperforms two compressed sensing models in two test cases, especially when using higher acceleration factors.