A practical guideline for T 1 reconstruction from various flip angles in MRI

Abstract

A fast and robust method for T₁ estimation in MRI is the so-called variable flip angle technique. We introduce a novel family of T₁ reconstruction methods from data acquired with various flip angles and propose a family member which combines the robustness of a nonlinear- with the computational advantages of a linear reconstruction. The constructed family contains the most common approaches for T₁ estimation, namely a linear and a nonlinear approach. A general sensitivity analysis for arbitrary members of the family is established. Advantages of the optimized reconstruction are demonstrated on phantom- as well as real data, showing improvements of up to 24% as compared with the linear method. As a further means to stabilize T₁ estimation, spatial stabilization methods are compared. We demonstrate on phantom and on real data that improved results can be obtained if not only T₁ but also a second unknown M₀ in the reconstruction is stabilized.

Keywords

T1 estimation variable flip angles sensitivity spatial stabilization

Introduction

Dynamic Contrast Enhanced MRI (DCE-MRI) is an imaging technique, which can be used to quantify physiological function, as e.g. the glomerular filtration rate (GFR), a main indicator of kidney performance.^1,2 A first, crucial step for such an analysis is the estimation of the baseline longitudinal relaxation times T₁ of the tissue prior to the administration of contrast agent.^3–5 Knowledge of the baseline T₁ map makes it possible to determine the amount of contrast agent at a given location during the DCE-MRI sequence.^1,2,6 After the amount of contrast agent in the tissue has been determined, methods such as pharmacokinetic modelling can be used to extract the functional tissue parameters.^1,2,7,8

It is well known that it is crucial to determine the baseline T₁ map accurately, as errors in T₁ can lead to significant misestimation of the reconstructed functional tissue parameters.⁹ However, there are several factors which make T₁ estimation delicate: first, T₁ cannot be measured directly and needs to be reconstructed from data.^10–12 Second, T₁ maps are ideally determined for each individual separately as they depend on tissue properties and can vary between individuals.¹³ In the special case of abdominal MRI of, e.g., kidney or prostate, an additional requirement is measurement speed, as motion by breathing or heart-beat complicates the acquisition process.^9,14

A T₁ reconstruction method which fulfils these requirements is the variable flip angle method.^3,4,10–12 The variable flip angle method has been shown to yield both fast and accurate results and is based on so-called spoiled gradient echo (SPGR) sequences, which are common in clinical MRI.¹⁵ For the variable flip angle method, a number of two or more SPGR MR images is acquired with different flip angles.^10–12 T₁ is subsequently reconstructed by fitting a model to the data.

However, there are several issues which complicate T₁ measurements with the variable flip angle method: the use of fast imaging techniques and weak signal can lead to low-quality data with low signal to noise ratio.¹⁶ Additionally, we show in Lemma 4 that T₁ reconstruction is unstable in areas with low steady-state magnetization or high T₁. Stable reconstruction methods are hence crucial to avoid outliers. This work will focus on two reconstruction techniques: in the nonlinear method,¹⁰ a nonlinear optimization problem has to be solved. This method not only yields the most stable results^4,5 but also comes with some computational overhead, see Wang et al.¹⁷ for a GPU implementation. Another approach is to reformulate the nonlinear problem as a linear problem and to solve the reconstruction as a linear regression.^11,12 For the regression approach, there is a closed-form solution formula,^4,11,12 but as noise influences both the abscissa and the ordinate, the stability of the reconstruction suffers.⁴

In this article, we discuss different ways to stabilize and to simplify T₁ reconstruction. First, we present a novel framework, which unifies the nonlinear and the linear reconstruction method. We then determine a linear reconstruction method, which approximates the nonlinear one more closely than the standard linear reconstruction. Specifically, we introduce a family of weighted linear reconstruction methods and show that the nonlinear reconstruction method corresponds to a choice of nonlinear weights. We then present a general approach to evaluate arbitrary weighted reconstruction techniques, which is based on the implicit function theorem and the error propagation theorem. An evaluation is performed for the sensitivity of the reconstructed T₁ with respect to perturbation of the data. For the special case of the nonlinear and of the linear weights, our findings are confirmed by the results obtained in Wang et al.⁴ and Deoni et al.⁵ Based on the computed sensitivities, we then introduce linear weights which approximate the nonlinear ones. Our experiments on phantom and on real data confirm that the linear approximated weights allow to reconstruct T₁ using closed-form solutions with errors similar to the nonlinear approach. The resulting reconstruction method thus combines the speed and simplicity of the linear reconstruction with the robustness of the nonlinear reconstruction.

As an additional means to stabilize T₁ reconstruction, we also evaluate a class of approaches which add spatial stabilization to the reconstruction. We compare joint approaches, where reconstruction and spatial stabilization are performed simultaneously^13,18 with denoising approaches, and where spatial stabilization is performed subsequently.¹⁶ As spatial stabilization might be performed only with respect to T₁,²³ we also investigate the impact of additional stabilization of the second unknown of the reconstruction, the constant M₀, cf. ‘Data fit for a single location and weighting strategies’ section. On phantom experiments, we find that the joint approaches are only superior to the denoising approaches if stabilization is performed with respect to both variables. Also, we compare stabilization with respect to two different parametrizations used for T₁ reconstruction, cf. ‘Data fit for a single location and weighting strategies’ section, Wang et al.⁴ and Christensen et al.¹⁰ Here, we find that stabilization in the $z = (T 1, M 0)$ coordinate system yields more stable results as compared with stabilization in an alternative coordinate system $z' = (E 1, N)$ , which is often used for T₁ reconstruction.^4,11,12 We conclude by showing improved T₁ reconstruction results on real data by using a total variation (TV) approach, which is stabilized both with respect to T₁ and M₀.

The Matlab code used to determine the optimal weights is made available under http://www.mic.uni-luebeck.de/people/constantin-sandmann.html.

Stabilization methods for T1 estimation

In this section, we describe different ways to compute tissue properties from MR measurements. This section is organized as follows: in ‘Data fit for a single location and weighting strategies’ section, we introduce a family of T₁ reconstruction methods. This family includes the two most commonly approaches for T₁ reconstruction, namely a linear and a nonlinear approach. In ‘Sensitivity of the weighted reconstruction’ section, we introduce a general framework under which we can analyse the sensitivity of the reconstruction methods with respect to errors in the data. ‘A linear approximation of the nonlinear weights’ section is based on this analysis. Here, we introduce a novel method to construct family members which yield robust reconstruction but are still simple to apply.

In ‘Spatial stabilization’ section, we describe spatially stabilized reconstruction methods. Namely joint methods, where T₁ reconstruction and spatial stabilization are performed simultaneously^13,17 and denoising methods, where spatial stabilization is performed subsequently to the reconstruction.¹⁶

Data fit for a single location and weighting strategies

In this section, we describe different ways to compute tissue properties from MR measurements. Particularly, we are interested in computing the unknown longitudinal MR relaxation time T₁ and a constant M₀ from MR data $d = (d 1, \dots, d K)$ acquired with a known repetition time $TR$ and various flip angles $α = (α 1, \dots, α K)$ ; the latter being considered as control parameters. Note that data, control parameters and unknowns also depend on a spatial location x. However, we start our outline by freezing the location and address the spatial dependency explicitly in ‘Data fit for multiple locations’ section.

Our forward model is the so-called signal equation S(z) = d.¹⁹ The signal equation describes the magnitude of a measured MR signal in dependency of the local tissue parameters and is hence valid in every voxel. It is phrased in the variable z, which can either be $z = (T 1, M 0)$ or $z = (E 1, N)$ . These are related by the following change of coordinates:

E 1 : = e - TR / T 1, N : = M 0 (1 - E 1) or T 1 = - TR / \log E 1, M 0 = N / (1 - E 1)

For given α and $TR$ , the signal equation for spoiled gradient echo sequences reads¹⁹

S (z) = M 0 \sin α \frac{1 - e - TR / T 1}{1 - \cos α e - TR / T 1} = \frac{\sin α N}{1 - \cos α E 1}

(1)

Introducing the E₁ (or T₁) dependent factors $ω_{k}^{NL} : = \sin α k / (1 - \cos α k E 1)$ , in each component $k \in {1, \dots, K}$ , the residual $S (z) - d$ can be decomposed into a linear part and a weighting:

\frac{\sin α k N}{1 - \cos α k E 1} - d k = ω_{k}^{NL} (d k / \tan α k E 1 + N - d k / \sin α k) = ω_{k}^{NL} (A k z - b k)

where the k-th row of

A \in ℝ K, 2

and

b \in ℝ K

are defined by

A k : = (d k / \tan α k, 1)

and

b k : = d k / \sin α k

A straightforward approach to compute the unknowns z is therefore to solve for each spatial location a weighted K-by-2 least squares problem:²⁰

Minimize J (z) : = ‖ Az - b ‖ W 2 = \sum_{k = 1}^{K} ω_{k}^{2} (A k z - b k) 2

(2)

Gupta¹¹ proposed to use $ω k = 1$ , which results in a standard linear least squares problem with an explicit solution formula $z = φ (d) = A † b, A †$ the Moore–Penrose pseudo-inverse.²⁰ Note, however, that lower errors in reconstruction can be achieved using nonlinear weighting, as noise influences A_k as well as b_k.^3,4,17 In ‘A linear approximation of the nonlinear weights’ section, we introduce linear weights which approximate the nonlinear weights. As the linear weighting method comes with a closed form solution formula, it is preferable from a practical point of view.

We summarize our model assumptions and conclude with a Lemma stating existence of an optimal parameter pair for the linear weights:

Assumption 1 We assume that there exist at least two measurements with different flip angles, i.e. K > 1 and $α 1 \neq α 2$ . For the weights, we assume that $ω \in ℝ K$ with $ω i > 0$ . For the parameters, $z \in ℝ 2$ assume that $z i > 0$ .

Lemma 1 (Linear Reconstruction). Assume that ω is independent of z. Then there exists a function $φ W : ℝ K \to ℝ 2, φ W (d) : = A † b$ , returning the unique least squares solution of the linearly weighted data fit.

Proof. Follows from the properties of the pseudo inverse, cf. Lawson and Hanson.²⁰□

Sensitivity of the weighted reconstruction

In this section, we compute the sensitivities of the introduced weighted reconstruction techniques in dependency of the weights ω. More precisely, given some parameters $z *$ , we will calculate the variance of T₁, if it is reconstructed from noisy measurements $d \propto N (S (z *), σ d)$ . Let us hence assume that $d * = S (z *)$ and that we have a set of noisy observations $d = d * + n$ with white Gaussian noise n. Then the function $φ (d) : = argmin z | | Az - b | | W 2$ maps an observation d to the corresponding reconstructed parameter set z. The sensitivity of the reconstruction method φ is then given by $\nabla φ (d *)$ . Even more, the error propagation theorem states that if $d \propto N (d *, σ d)$ , the reconstructed parameters $z = (z 1, z 2)$ are distributed as $z i \propto N (z_{i}^{*}, σ z i)$ for $σ z i = σ d | \nabla φ i |$ .²¹ In this section, we present a unified framework to calculate $| \nabla φ i |$ for the complete family of approaches introduced in ‘Data fit for a single location and weighting strategies’ section. Note that all these approaches can be written as $φ (d) : = argmin z J (z, d)$ for some suitable J. We will obtain $\nabla φ$ by means of the following observation: by construction, the function $f (z, d) : = \nabla z J (z, d)$ is zero for all pairs $φ (d) = z$ in a neighbourhood of $d *$ since $φ (d)$ minimizes $J (\cdot, d)$ . It follows that $\nabla d [f (φ (d), d)] = 0$ . Application of the chain-rule yields $\nabla z f \cdot \nabla φ + \nabla d f = 0$ . Assuming that $\nabla z f$ is invertible, we obtain by another straight-forward calculation $\nabla φ = - [\nabla z f] - 1 \nabla d f$ . Note that the inverse function theorem²² gives the exact assumptions under which the above calculus is justified.

To begin with, we introduce the functions $f W, f NL : ℝ 2 \times ℝ K \to ℝ 2$ with

f W (z, d) : = \nabla z ‖ Az - d ‖ W 2 and f NL (z, d) : = \nabla z ‖ S (z) - d ‖ 2

(3)

Lemma 2 Let z be given with $S (z) = : d$ . Under Assumptions 1, the matrices $\nabla z f W (z, d), \nabla z f NL (z, d) \in ℝ 2, 2$ have full rank; cf. (3).

Proof. We start by noting that for the linear reconstruction, the weights do not depend on z. Let $W : = diag (ω 1, \dots, ω k)$ . A direct calculation shows

\frac{1}{2} \nabla z f W = A ⊤ WA

Since $c \tan α i = d i$ for all α_i if c = 0 and d = 0, it follows by assumption that A has full rank. The statement thus follows from W being a positive definite diagonal matrix.

For the nonlinear case,

\frac{1}{2} \nabla z f NL = [\nabla z S] ⊤ [\nabla z S] - (S (z) - d) ⊤ \nabla_{z}^{2} S = [\nabla z S] ⊤ [\nabla z S];

the simplification due to S(z) = d; the 2-by-2-by-K tensor $\nabla_{z}^{2} S$ is not explicitly outlined. The K-by-2 matrix $\nabla z S$ reads

\nabla z S = [d i s (α i), d i / M 0] i = 1 K with s (α i) = \frac{(\cos α i - 1) (\partial E 1 / \partial z 1)}{(1 - \cos α i E 1) (1 - E 1)}

Since there are at least two measurements with different flip-angles and $d \neq 0$ , it follows that $\nabla z S$ and hence $\nabla z f$ has full rank.□

Lemma 3 (Nonlinear Reconstruction). Let $z *$ be given and $S (z *) = : d *$ . Under Assumption 1, there exist neighbourhoods $U \subseteq ℝ 2$ of $z *$ and $V \subseteq ℝ K$ of $d *$ and a function $φ NL : V \to U$ returning a unique minimizer of the nonlinear data fit and

\nabla φ NL (d *) = - [\nabla z f NL (z *, d *)] - 1 \nabla d f NL (z *, d *)

(4)

Proof. As it follows from Lemma 2 that $f NL (z *, d *) = 0$ and that $\nabla z f NL (z *, d *)$ has full rank, the claim follows directly from the inverse function theorem.²²□

Lemma 4 (Instability of the Reconstruction). Assume that $M 0 = 0$ or $T 1 = + \infty$ . Then, T₁ cannot be reconstructed uniquely.

Proof. In both cases, it holds that S(z) = 0. As S(z) = 0 for any $z = (z 1, 0)$ and $z 1 \in ℝ$ , no unique reconstruction is possible.□

We are now ready to state the main result. To this end, we fix the coordinate system as $z = (T 1, M 0)$ . We start with a fixed $z * = (T_{1}^{*}, M_{0}^{*})$ with corresponding perfect data $d * = S (z *)$ . Moreover, we assume that $d \propto N (d *, σ d)$ are realizations of $d *$ with known standard deviation σ_d. We will show that the reconstructed T₁ is distributed as $T 1 \propto N (T_{1}^{*}, σ T 1)$ , where the $σ T 1$ is a product of the standard deviation of the data and the sensitivity of φ (Figure 1).

Figure 1.

Nonlinear reconstruction: measured and theoretical distribution of T₁ recovered from 65,000 noisy measurements using the nonlinear data-fit. The data were simulated using the same parameters as in Wang et al.,⁴ i.e. $T 1 = 1000 ms, M 0 = 1000, TR = 800 ms$ and flip-angles $α 1 = 10 \circ$ to $α 10 = 100 \circ$ in increments of $10 \circ$ . White Gaussian noise had a standard deviation $σ d = 10 ms$ . In this experiment, the measured predicted standard deviation both were $σ T 1 = 32.3 ms$ on an error level of $10 - 1 ms$ . The standard deviation for the linear fit was $45.3 ms$ .⁴

Theorem 1 Let Assumption 1 be satisfied, $z * = (T_{1}^{*}, M_{0}^{*})$ and $d *$ with $S (z *) = d *$ . Furthermore, let $d \propto N (d *, σ d)$ with standard deviation σ_d. Then, $T 1 \propto N (T_{1}^{*}, σ T 1)$ with

σ_{T 1}^{2} = σ_{d}^{2} ‖ \nabla φ 1 (S (z *)) ‖ 2

(5)

where

φ (d) = (φ 1 (d), φ 2 (d))

and

\nabla φ NL (d *) = 2 [\nabla z f NL (z *, d *)] - 1 \nabla z S (z *) and \nabla φ W (d *) = 2 [\frac{dz}{dz'}] [\nabla z f W (z *, d *)] - 1 WA

Here $\frac{dz}{dz'}$ denotes the diagonal matrix $diag (T_{1}^{2} / (E 1 TR), 1 / (1 - E 1))$ which accounts for the change of variables.

Proof. Equation (4) states that

[\nabla φ 1, \nabla φ 2] ⊤ = - [\nabla z f] - 1 [\nabla d f]

Matrices $\nabla z f W$ and $\nabla z f NL$ are given in the proof of Lemma 2. A direct calculation shows that for the linear case $\frac{1}{2} \nabla d f W = - WA$ and for the nonlinear case $\frac{1}{2} \nabla d f W = - \nabla z S$ . Note that the linear case is formulated in the $z' = (E 1, N)$ coordinate-system. Hence, the derivative needs to be scaled by the chain-rule and $\nabla z φ W = \frac{dz}{dz'} \cdot \nabla z' φ W$ . Equation (5) now follows directly from the error propagation theorem.²¹□

To illustrate the situation, we reconstructed T₁ from 65,000 noisy measurements. We then compared the measured standard deviation in T₁ with the one given by $σ T 1$ . Both values matched on an error level of $10 - 1 ms$ . Note, however, that it is challenging to estimate the standard deviation from the experiments reliably. Even with significantly larger sample-sizes, the experimentally derived values still varied about $5 \cdot 10 - 2$ and tended to be underestimated by the predicted value.

A linear approximation of the nonlinear weights

Theorem 1 gives an expression which describes the sensitivity of the reconstruction for a set of parameters z and weights ω. Note that the obtained expression is basically a quadratic function in ω. We define the set of admissible weights $Σ : = {ω \in ℝ L, ω \geq 0, \sum ω i = 1}$ . The normalization is necessary since the reconstruction described in (2) is invariant under scaling of the weights. Given a set of parameters $z 1, \dots, z L$ , we propose to solve the following minimization problem to obtain a set of optimized weights, which minimizes the joint variance for $z 1, \dots, z L$ :

ω Opt = argmin ω \in Σ {\sum_{j = 1}^{L} p j σ_{T 1}^{2} (z j, ω)}

(6)

Here, $p j \in ℝ_{0}^{+}$ are user-defined factors which can be used to highlight particular ranges of parameters z_j in the reconstruction. The optimized weights can be used to reconstruct T₁ by solving a linear system but with improved sensitivity of the reconstruction.

To illustrate the setting, we visualize the standard deviation landscape in dependency of the weights for $T 1 = 1000 ms$ and $M 0 = 1000$ in Figure 2. We then solved problem (6) to obtain a set of weights with minimal variance. Results are showing $ω Opt$ different from $ω NL$ , but with the same standard deviation.

Figure 2.

Map displaying $σ T 1$ for different weights $ω = (ω 1, ω 2, 1 - ω 2 - ω 3) \in ℝ 3$ and fixed T₁ and M₀. Control parameters were $T 1 = 1000 ms, M 0 = 1000, α = (5 \circ, 8 \circ, 15 \circ)$ and $TR = 5 ms$ . Highlighted by a white × are the nonlinear weights, linear weights and the optimized weights. Standard deviations for the optimized weights and for the nonlinear weights were $σ_{T 1}^{Opt} \approx σ_{T 1}^{NL} \approx 168.21 ms$ with $σ_{T 1}^{Opt} - σ_{T 1}^{NL} = 2.06 \cdot 10 - 11 ms$ . Standard deviation for the linear weight was $σ_{T 1}^{Lin} = 184.34 ms$ . For this configuration, $σ T 1$ can be reduced by 9% using the optimized or the nonlinear weighting scheme instead of the standard linear one.

Data fit for multiple locations

Since T₁ reconstruction typically aims to reconstruct spatially arranged parameter maps, we now extend the notation to cover all data. To this end, we recall that the data and unknowns are presented on a discrete m₁-by-m₂-by-m₃ grid of $n : = m 1 m 2 m 3$ points x_j, $j = 1, \dots, n$ . Thus, $T 1, M 0, E 1, N \in ℝ n, d \in ℝ Kn, z \in ℝ 2 n$ with $z j = E 1, j$ and $z n + j = N j$ .

The extension of equation (2) can be phrased in terms of where

A j = (d j / \tan α 1 1 : : d j / \tan α K 1) \in ℝ K, 2, b j = (d j / \sin α 1 : d j / \sin α K) \in ℝ K, W j = diag (ω_{j, k}^{2}, k = 1, \dots, K) .

Hence, with $z j = (z j, z n + j)$ , we now have

‖ Az - b ‖ W 2 = \sum_{j = 1}^{n} ‖ A j z j - b j ‖ W j 2 = \sum_{j = 1}^{n} \sum_{k = 1}^{K} ω_{j, k}^{2} (A_{k}^{j} z j - b_{k}^{j}) 2

= \sum_{j = 1}^{n} \sum_{k = 1}^{K} ω_{j, k}^{2} (d_{k}^{j} / \tan α k E 1, j + N j - d_{k}^{j} / \sin α k) 2

Spatial stabilization

Our experiments indicate that all the approaches outlined in ‘Data fit for a single location and weighting strategies’ section suffer from insufficient data such as zero background values or low signal-to-noise ratio especially at low flip-angles, see also Lemma 4 and Wang and Cao¹⁸ and de Pasquale et al.²³ A remedy is to further stabilize the reconstruction by including information on the spatial setup of the T₁ maps. Information on the layout of the expected parameter maps is often expressed as the property to be smooth²³ or piecewise constant.¹⁸

A common way to include these properties is by using a joint approach, where stabilization is introduced by the addition of a stabilization term R to the reconstruction problem. For some user-defined parameter $λ \in ℝ +$ , the stabilized reconstruction problem hence reads

J (z) = ‖ S (z) - b ‖ 2 + λ R (z)

For spatial stabilization, we will always use the nonlinear data-term, as this yielded the most stable reconstruction. However, it is not clear if the problem should be formulated in the (T₁, M₀)^18,23 or in the $(E 1, N)$ ¹⁴ realm. Our experiments will hence cover both choices of coordinate systems. For ease of the presentation, we will discuss the stabilization in a functional setting. With $z = (z 1, z 2) : ℝ 3 \to ℝ 2$ , we will compare the following three choices for stabilization:

R (z) = R Diff (z) : = ‖ \nabla z ‖ γ 2 : = \int γ 1 | \nabla z 1 | 2 + γ 2 | \nabla z 2 | 2 d x

(7)

R (z) = R TV (z) : = \int γ 1 | \nabla z 1 | 2 + γ 2 | \nabla z 2 | 2 d x

(8)

L₂ stabilization for the gradient of T₁ as in (7) has been proposed in a Bayesian setting for inversion time estimation.²³ An advantage of this method is that it is easy to implement, as standard nonlinear-least squares methods can be used for optimization.²⁰ A downside of this approach is that it is usually not edge preserving. The idea to include TV to stabilize T₁ estimation has been introduced in Wang and Cao.¹⁸ With TV stabilized methods, edges in the T₁ map are preserved. However, the numerical minimization requires more sophisticated methods.^24,25

Additionally, we will also cover denoising approaches, where z is reconstructed using a standard, nonlinear parameter fit. Subsequently to the reconstruction, the result is denoised by solving

Minimize J (u) : = ‖ z - u ‖ 2 + R (u)

(9)

where R again is defined as above.

Numerics

As numerical methods were not in the focus of this work, we only give a brief overview. All implementations were done in Matlab on a standard PC with 3.4 GHz and 16 GB RAM. Discrete Gradient-Operators were built using short forward-differences with Neumann boundary conditions.

We start by describing the implementation of the unstabilized reconstruction methods: for the weighted reconstruction with linear weights (2), the resulting normal equation was solved using an explicit representation of the inverse matrix. The nonlinear reconstruction was solved using a projected Gauss–Newton Method as described in Bertsekas.²⁶ To avoid instabilities in the reconstruction, we used the constraints $1 ms \leq T 1 \leq 5000 ms$ and $M 0 \geq 0$ . The optimized weights (6) were obtained using a projected Levenberg–Marquardt method tailored for optimization on the positive unit-simplex.²⁶ For the joint approaches, we used the following implementations: the joint reconstruction with diffusive stabilization (7) was solved using a standard Gauss–Newton method. The joint reconstruction with TV stabilization (8) was also solved in a Gauss–Newton fashion: the signal equation was linearized and the resulting linear TV-Denoising problems were solved using a preconditioned Primal-Dual Hybrid-Gradient algorithm.^24,25 For the denoising approaches, we used the following numerics: the diffusive denoising problem was solved explicitly using the built-in Matlab solver. The TV denoising problem was again solved using a preconditioned Primal-Dual Hybrid-Gradient algorithm.^24,25

Results

Evaluation of T₁ reconstruction methods on real data is delicate: T₁ cannot be measured directly, it depends on tissue type and also on individual characteristics.¹³ For real data, a ground-truth T₁ map is hence generally not available. To cope with this issue, we tested the proposed methods both on simulated data, where we compared our results to the ground-truth, and on real data, where we compared our results with results of the purely nonlinear reconstruction. The simulated data were set-up using a tailored phantom of the human kidney, cf. ‘An XCAT software phantom for the human kidney’ section. Real data came from a study on the function of the human kidney²⁷ and was made available by courtesy of Jarle Rørvik from the Haukeland University Hospital in Bergen, Norway.

Note that pre-alignment of variable flip-angle data is challenging, as intensities vary with different flip-angles. To overcome this issue, a model-based registration method similar to the ones described in Hallack et al.⁹ and Heck et al.¹⁴ was used for registration. Note that our modification has some mild smoothness assumptions on the parameter maps for stability.

As a metric for the success of T₁ estimation on phantom data, we used the relative error (RE), defined as the fraction of the error of the true T₁, i.e.

RE i : = \frac{| T_{1}^{i, true} - T_{1}^{i, recon} |}{T_{1}^{i, true}}

This error was measured pointwise, for regional analysis, we employed the mean RE defined as the mean of all pointwise REs.

An XCAT software phantom for the human kidney

In this section, we will describe the construction tailored software phantom for the human kidney. A simulated 2D MRI scan of the human kidney was set up as follows: the anatomy (cortex, medulla, background) was obtained from the XCAT-Phantom.²⁸ Estimates of T₁ of the mentioned structures were taken from literature.¹³ The second constant M₀ was determined by $M 0 = s 0 exp (- TE / T 2)$ for some constant s₀, which was approximately determined from our own measurements.³ More precisely, we chose for the kidney cortex $T 1 = 966 ms, T 2 = 87 ms$ and $s 0 = 3100$ . The medulla was set up with parameters $T 1 = 1412 ms, T 2 = 85 ms$ and $s 0 = 3500$ and for the background, we chose $T 1 = 110 ms, T 2 = 85 ms$ and $s 0 = 2000$ .¹³ Partial volume effects were simulated by smoothing the parameter maps M₀ and T₁ with a Gauss-filter of width 3 voxel and standard deviation of 0.7 voxel. Finally, we used the signal-equation (1) to simulate (perfect) MR signals with scan parameters taken from different publications on T₁ estimation. Sequence S1 was designed to yield robust reconstruction over a large range of T₁,³ S2 and S3 were chosen from a publication on the nonlinear reconstruction method and S4 was used to acquire the real data presented in this work. Examples of the phantom images are displayed in Figure 3.

Figure 3.

MR sequence parameters and corresponding sources from literature. Also given is the software phantom for S₄ with a noise intensity of 7%, cf. ‘An XCAT software phantom for the human kidney’ section.

To add more realistic conditions, white Gaussian noise with variance independent of the flip-angle was added to the k-space data.²⁹ Noise intensity was measured pointwise, as the fraction of noise of the undisturbed signal. For a simulation sequence determined by $α, TR$ and $TE$ , the effective noise intensity is given by the mean of all pointwise intensities over all flip angles. Similar to Wang et al.,¹⁷ in our experiments, effective noise intensities of 3% , 5% , 7% and 9% were simulated.

Results of the weighting strategies

In this section, we show how the proposed weighting strategies can improve T₁ estimation. It is well-known that the sensitivity of the reconstruction does not only depend on the reconstruction method and the expected (T₁, M₀) but also on the employed sequence parameters.^3–5 The proposed weighting strategies (Nonlinear, Linear, Optimized) were evaluated with respect to their theoretical sensitivities and with respect to the mean error on simulated as well as real data. As expected, the results depend on the employed sequence: for sequences S2 and S4, we can observe large positive impact of weighting strategies on reconstruction with reductions in the standard deviation of up to 24%, for sequences S1 and S3, improvements are only minor. In our experiments, we found that the optimized weights can perform similar to the purely nonlinear weights with a fraction of the computational overhead. The phantom data results are confirmed by the real data results, where differences between the nonlinear reconstruction and the linear reconstruction with optimized weights were in the range of 1%.

More thoroughly, we calculated $σ T 1$ for the nonlinear reconstruction method, the linear reconstruction method and the optimized reconstruction method for the sequences described in Figure 3. We found that the optimized linear weights theoretically perform close to the nonlinear weights, see Table 1. However, note that theoretical sensitivities are given with respect to a purely Gaussian noise model, whereas in MRI, we expect Rician noise.²⁹ We hence additionally tested the different weighting schemes for phantom data degraded by Rician Noise. In Table 2, we are giving reconstruction results on phantom data with respect to the mean RE. It can be seen that, also for Rician Noise, the optimized weights improve the errors in reconstructed T₁ drastically. However, as expected, the improvement depends largely on the sensitivities of the nonlinear and the linear reconstruction: if these are very similar, as it is the case for S1, only minor influence on the reconstruction can be observed (Figure 4).

Figure 4.

Displaying results of the nonlinear, linear and the optimized reconstruction for real data. Given is the mean relative error of T₁ in percent as compared with the results of the nonlinear reconstruction. Weights were obtained by solving (6) with $T 1 = (1412, 966), M 0 = (3000, 3000)$ and $f = (0.5, 0.5)$ . The images depict the reconstruction of T₁ for D₁, left kidney. Note that T₁ cannot be estimated reliably in regions with zero signal intensity, as it is the case for the lung, cf. Lemma 4. Visual inspection confirms the close agreement of the nonlinear and the optimized approach.

Table 1.

The predicted sensitivities of the different weighting schemes for the sequences described in Figure 3.

	Cortex			Medulla			Background
	$ω Lin$	$ω NL$	$ω Opt$	$ω Lin$	$ω NL$	$ω Opt$	$ω Lin$	$ω NL$	$ω Opt$
S1	18.29	17.70	17.70	40.39	37.12	37.12	2.71	1.61	1.92
S2	24.84	20.65	20.66	62.38	48.53	48.53	1.99	1.46	1.68
S3	28.61	27.57	27.57	80.28	76.09	76.09	2.10	1.99	2.00
S4	57.30	45.82	45.82	182.38	139.28	139.28	1.31	1.28	1.37

Given is the standard deviation in reconstructed T₁ with respect to $σ d = 1 ms$ , cf. (5). Factors p_j were selected proportional to the expected amount of cortex, medulla and background. Results from the nonlinear and the optimized weighting scheme are in close agreement.

Table 2.

The improvement of the mean-relative error in reconstructed T₁ by use of the weighting schemes.

		3%	5%	7%	9%
S1	$ω Lin$	5.48	9.06	12.68	16.74
	$ω NL$	5.23	8.62	12.03	15.75
	$ω Opt$	5.24	8.66	12.13	16.05
S2	$ω Lin$	4.09	6.81	9.47	12.37
	$ω NL$	3.34	5.56	7.72	10.05
	$ω Opt$	3.34	5.58	7.76	10.15
S3	$ω Lin$	10.76	18.59	43.45	96.91
	$ω NL$	10.20	17.29	25.43	35.71
	$ω Opt$	10.38	18.01	38.22	89.31
S4	$ω Lin$	11.43	51.41	89.74	182.46
	$ω NL$	8.72	22.72	24.79	34.86
	$ω Opt$	8.81	21.30	38.53	113.57

Experiments were conducted with different sequences and different noise-levels on phantom data, cf. ‘An XCAT software phantom for the human kidney’ section. Given is the mean relative error in percent, calculated in 5000 experiments and evaluated over cortex and medulla. Factors p_j were selected as described in Table 1. The nonlinear parameter fit yields the smallest and the standard linear fit the largest errors. Results of the optimized weights are in close agreement with the nonlinear weights for small noise. For high noise level, results of the optimized weights and the nonlinear weights begin to differ, possibly due to the less Gaussian character of the noise.

Analysis of stabilization strategies

In our experiments, we address multiple questions for stabilization in T₁ estimation:

Stabilization with respect to z₁ versus stabilization with respect to $z = (z 1, z 2)$ .

Stabilization in $z = (T 1, M 0)$ versus stabilization in $z' = (E 1, N)$ .

Stabilization during reconstruction (Joint) versus stabilization after reconstruction (Denoising).

Diffusive stabilization versus TV stabilization.

Experiments were performed on the software phantom described in ‘An XCAT software phantom for the human kidney’ section. Sequence parameters S1 were used for all experiments, as these yielded the most stable results, cf. ‘Results of the weighting strategies’ section. We assumed a modest noise-level of 7%. The experiments thus demonstrate the advantages of stabilization even for well-designed sequences. Optimal stabilization parameters $(γ 1, γ 2)$ were determined experimentally by sampling the parameter space. Results were obtained using Monte-Carlo simulations with 5000 reconstructions from random measurements until the mean RE was approximately stable up to $10 - 5$ .

Results are presented in Table 3. It can readily be seen that all discussed stabilization methods are capable to drastically reduce the RE as compared with no spatial stabilization. Also, stabilization in both T₁ and M₀ has advantages as compared with stabilization only in T₁. However, these effects cannot be observed for stabilization in the

(E 1, N)

coordinate system. In this case, no large difference in stabilization of

(E 1, N)

as compared with stabilization only in E₁ can be observed. This might be due to the exponential scaling, as variations in high T₁ are less accurately detectable in E₁. Comparing the joint approaches with the denoising approaches, one can see the joint approaches are only superior if both parameters are stabilized. This establishes the need for joint approaches considering both variables M₀ and T₁.

Table 3.

The impact of spatial stabilization on phantom data for sequence S1 and noise-level 7%.

$#$ Parameters	Joint				Denoising
	$R TV$		$R Diff$		$R TV$	$R Diff$	No Stab.
	1	2	1	2	1	1	No Stab.
$z = (T 1, M 0)$	6.76	5.96	7.36	6.59	6.31	6.74	12.91
$z' = (E, N)$	6.48	6.49	11.32	11.32	6.58	9.71	12.91

Results are given as mean relative error calculated from 5000 experiments on cortex and medulla and are given in percent. It can be observed that phantom data joint approaches are only superior to denoising approaches if both parameters are regularized. Also, regularization in (T₁, M₀) is superior to regularization in (E₁, M₀).

As in the case of the kidney, TV regularization turned out to yield the lowest REs on phantom data, we additionally tested the reconstruction methods on real data. Note, however, that the stabilization needs to be adjusted to fit to the structure of the data and other data might benefit from diffusive regularization. Parameters were selected manually with respect to visual plausibility of the parameter maps. Evaluation was done with respect of the mean RE on cortex on medulla. The RE was calculated voxel-wise with respect to the mean of the nonlinear approach on the regions of interest. We found the largest impact of regularization on high T₁, cf. Figure 5. This is expected as it is difficult to estimate high T₁ stably, cf. Lemma 4. However, it is difficult to adjust both parameters yielding a possible over-smoothing of M₀.

Figure 5.

Figure displaying the impact of spatial stabilization on real data. Given is the mean relative error in percent. The error was computed voxel-wise with respect to the mean of the nonlinear reconstruction over the area of interest. For stabilization, we chose $R TV$ . Comparison was performed with respect to joint approach with regularization in (T₁, M₀), joint approach with regularization in T₁ and denoising approach. Optimal stabilization parameters were determined with respect to visual plausibility of the parameter maps. As compared with the purely nonlinear approach, all stabilized results are more homogeneous. Advantages of the two parameter approach can be observed mainly in reconstruction in the medulla, where T₁ is high as compared with the cortex. An example of the reconstruction is given for dataset D3, right kidney.

Conclusion

A novel, unified framework for T₁ recovery was introduced. The framework consists of a family of weighted reconstruction techniques and includes the two most commonly used ones, namely the linear and the nonlinear approach, as special cases. Based on the inverse function theorem and the error propagation theorem, a general sensitivity analysis for all members of the proposed family was established. The sensitivities were used to determine a set of weights, which allows for linear reconstruction of T₁ with errors similar to the nonlinear approach. The optimized weights can be used for both fast and robust reconstruction of T₁, avoiding the computational overhead of nonlinear optimization. Results for synthetic as well as real data confirm improvements in T₁ estimation of up to 24% as compared with standard linear reconstruction. Matlab code to determine the optimized set of weights is made available under http://www.mic.uni-luebeck.de/people/constantin-sandmann.html.

We show that T₁ reconstruction becomes increasingly ill-posed for low signal intensities and therefore additionally analyse spatial stabilization methods. Here, it was found on phantom data that joint reconstruction and stabilization techniques are only superior to standard denoising techniques if both variables T₁ and M₀ are stabilized. Note that in some approaches for spatially stabilized T₁ recovery stabilization is applied only with respect to T₁.¹⁸ Our results suggest that results could additionally be improved by using stabilization of both variables. Also, the impact of stabilization in the $(E 1, N)$ instead of the (T₁, M₀) coordinate system was investigated. Here, it was found that stabilization in (T₁, M₀) is superior, as the change of coordinates introduces systematic biases for large T₁.

Footnotes

Acknowledgements

All real data were provided by courtesy of Jarle Rørvik and the kidney project of MedViz at the Haukeland University Hospital in Bergen, Norway. For providing the kidney segmentations, the authors would like to thank Katja Sandmann from the Sana Klinken Ostholstein, Klinik Eutin, Germany.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Sourbron

Michaely

Reiser

. MRI-measurement of perfusion and glomerular filtration in the human kidney with a separable compartment model. Invest Radiol 2008; 43: 40–48.

Tofts

Brix

Buckley

. Estimating kinetic parameters from Dynamic Contrast-Enhanced T1-weighted MRI of a diffusable tracer: standardized quantities and symbols. J Magn Reson Imaging 1999; 10: 223–232.

Cheng

H-LM

Wright

. Rapid high-resolution T1 mapping by variable flip angles: Accurate and precise measurements in the presence of radiofrequency field inhomogeneity. Magn Reson Med 2006; 55: 566–574.

Wang

Riederer

Lee

. Optimizing the precision in T1 relaxation estimation using limited flip angles. Magn Reson Med 1987; 5: 399–416.

Deoni

SCL

Peters

Rutt

. Determination of optimal angles for variable nutation proton magnetic spin-lattice, T1, and spin-spin, T2, relaxation times measurement. Magn Reson Med 2004; 51: 194–199.

Strijkers

Mulder

Willem

. MRI contrast agents: current status and future perspectives. Anti-Cancer Agent Me 2007; 7: 291–305.

Sourbron

. Technical aspects of MR perfusion. Eur J Radiol 2010; 76: 304–313.

Sourbron

Buckley

. Classic models for dynamic contrast-enhanced MRI. NMR Biomed 2013; 26: 1004–1027.

Hallack A, Chappell MA, Gooding MJ, et al. A new similarity metric for groupwise registration of variable flip angle sequences for improved T10 estimation in DCE-MRI. In: Ourselin, S and Modat M (eds) Proceedings of biomedical image registration 2014, London, 2014, pp.154–163. Cham, Switzerland: Springer International Publishing.

10.

Christensen

Grant

Schulman

. Optimal determination of relaxation times of Fourier transform nuclear magnetic resonance. Determination of spin-lattice relaxation times in chemically polarized species. J Phys Chem 1974; 78: 1971–1977.

11.

Gupta

. A new look at the method of variable nutation angle for the measurement of spin-lattice relaxation times using Fourier transform NMR. J Magn Reson 1977; 25: 231–235.

12.

Fram

Herfkens

Johnson

. Rapid calculation of T1 using variable flip angle gradient refocused imaging. Magn Reson Imag 1987; 5: 201–208.

13.

de Bazelaire

CMJ

Duhamel

Rofsky

. MR imaging relaxation times of abdominal and pelvic tissues measured in vivo at 3.0T: preliminary results 1. Radiology 2004; 230: 652–659.

14.

Heck C, Benning M and Modersitzki J. Joint registration and parameter estimation of T1 relaxation times using variable flip angles. In: Handels, H and Deserno, TM and Meinzer T-P (eds) Procedings of bildverarbeitung für die medizin 2015, Lübeck, 2015, pp.215–220. Heidelberg: Springer Berlin Heidelberg.

15.

Gourtsoyiannis

. Clinical MRI of the abdomen: why, how, when, Berlin Heidelberg: Springer, 2011.

16.

Manjón

Coupé

Martí-Bonmatí

. Adaptive non-local means denoising of MR images with spatially varying noise levels. J Magn Reson Imaging 2010; 31: 192–203.

17.

Wang

Shi

Wang

Y-XJ

. Concatenated and parallel optimization for the estimation of T1 maps in FLASH MRI with multiple flip angles. Magn Reson Med 2010; 63: 1431–1436.

18.

Wang

Cao

. Spatially regularized T1 estimation from variable flip angles MRI. Med Phys 2012; 39: 4139–4148.

19.

Zur

Stokar

Bendel

. An analysis of fast imaging sequences with steady-state transverse magnetization refocusing. Magn Reson Med 1988; 6: 175–193.

20.

Lawson CL and Hanson RJ. Solving least squares problems. Philadelphia: SIAM’s Classics in Applied Mathematics, 1995 (Originally published 1974 by Prentice-Hall).

21.

Cramér H. Mathematical methods of statistics. Bombay: Princeton University Press, 1946. Reprinted in 1961.

22.

Rudin

. Principles of mathematical analysis, 3rd ed. Auckland: McGraw-Hill, 1987.

23.

de Pasquale

Sebastiani

Egger

. Bayesian estimation of relaxation times T1 in MR images of irradiated Fricke-agarose gels. Magn Reson Imag 2000; 18: 721–731.

24.

Chambolle

Pock

. A first-order primal-dual algorithm for convex problems with applications to imaging. J Math Imaging Vis 2011; 40: 120–145.

25.

Pock T and Chambolle A. Diagonal preconditioning for first order primal-dual algorithms in convex optimization. In: Proceedings of IEEE ICCV 2011, Barcelona, 2011, pp.1762–1769. New York: IEEE.

26.

Bertsekas

. Projected Newton methods for optimization problems with simple constraints. SIAM J Control Optim 1982; 20: 221–246.

27.

Hodneland E, Kjørstad Å, Andersen E, et al. In vivo estimation of glomerular filtration in the kidney using DCE-MRI. In: Proceedings of ISPA 2011, Dubrovnik, 2011, pp.755–761. New York: IEEE.

28.

Segars

Sturgeon

Mendonca

. 4D XCAT phantom for multimodality imaging research. Med Phys 2010; 37: 4902–4915.

29.

Gudbjartsson

Patz

. The Rician distribution of noisy MRI data. Magn Reson Med 1995; 34: 910–914.