Greedy approximate projection for magnetic resonance fingerprinting with partial volumes

In this section, we introduce the notation we will use in the remainder of the paper. All the symbols used to describe the proposed method are summarised in table 1, provided at the end of section 3. We refer the reader to [23, 24] for additional details about optimisation. To have a compact notation when selecting a specific row $n \in \{1, \ldots, N\}$ of a matrix $\boldsymbol{M}\in\mathbb{C}^{N\times L}$, we use the notation $\boldsymbol{M}_{n,:}=\left(\boldsymbol{M}_{n,l}\right)_{1\leqslant l\leqslant L}$. Similarly, to select a specific column $l \in \{1, \ldots, L\}$ of this matrix, we use $\boldsymbol{M}_{:,l}=\left(\boldsymbol{M}_{n,l}\right)_{1\leqslant n\leqslant N}$. More generally, this notation is also used to select subparts of tensors. The operator ${{\rm real}}(\cdot)$ gives the real part of its complex argument, the operator ${{\rm Diag}}(\cdot)$ builds a diagonal matrix whose diagonal elements are given by its argument, and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}(\cdot){}^\dagger$ gives its adjoint. The adjoint of a linear operator $g \colon \mathbb{C}^L \to \mathbb{C}^N$ is denoted by $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}g^\dagger$. The cardinality of a countable set $\mathcal{T}$ is given by ${{\rm card}}(\mathcal{T})$. The $\newcommand{\e}{{\rm e}} \ell_p$ norm ($p\in ]0,+\infty]$) is denoted by $\|\cdot\|_p$. The $\newcommand{\e}{{\rm e}} \ell_0$ pseudo-norm [3], counting the non-zero entries of its argument, is defined as $(\forall \boldsymbol{x} \in \mathbb{R}^N) \quad \|\boldsymbol{x}\|_0 = \sum_{n=1}^N \left(\boldsymbol{x}_{n}\right){}^0$, with the convention 0⁰  =  0. By abuse of notation, the $\newcommand{\e}{{\rm e}} \ell_p$ norms and the $\newcommand{\e}{{\rm e}} \ell_0$ pseudo-norm will be used for tensors by reshaping them into vectors. Finally, the projection of a vector $\overline{\boldsymbol{x}}\in \mathbb{C}^N$ onto a non-empty closed subset $\mathcal{S}$ of $\mathbb{C}^N$ is given by $\mathcal{P}_{\mathcal{S}} (\overline{\boldsymbol{x}})= \mathop{\rm argmin}\nolimits_{\mathbf{\boldsymbol{x}\in \mathcal{S}}} \frac{1}{2}\|\boldsymbol{x}-{\overline{\boldsymbol{x}}}\|_2^2$ [23]. The same notation is used for projections of tensors.

Table 1. Table of symbols.

N	Number of voxels in the volume of interest
Q	Number of measurements per excitation and per coil
L	Number of excitation instances
$\boldsymbol{Y}\in\mathbb{C}^{Q\times L\times C}$	Measurement matrix
$\boldsymbol{M}\in \mathbb{C}^{N\times L}$	Magnetisation response of the volume of interest, introduced in (1)
h	Linear operator from $\mathbb{C}^{N\times L}$ to $\mathbb{C}^{Q\times L \times C}$ defining the acquisition process
P	Number of the tissue parameters
$\mathcal{M}\subset \mathbb{R}^{1 \times P}$	Subset describing the feasible parameter space
B	Non-linear smooth operator from $\mathcal{M}$ to $\mathbb{C}^{1\times L}$ describing the magnetic resonance experiment, introduced in section 2.2
$\boldsymbol{\Phi} \in \mathbb{C}^{D\times L}$	Discretisation of B with D elements
$\mathcal{B}^+$	Set describing $\boldsymbol{M}$ as a magnetisation response of pure voxels
$\boldsymbol{X}\in \mathbb{R}_{+}^{N\times D}$	Mixing matrix used to describe the PVE, introduced in section 2.3
$\mathcal{B}_{\mathcal{S}_+}(\boldsymbol{\Phi})$	Set describing all possible $\boldsymbol{M}$ satisfying the proposed PV model, defined in (4)
$\mathcal{S}_+$	Set describing the intersection of the sets $\mathcal{S}_1$, $\mathcal{S}_2$, $\mathcal{S}_3$ and $\mathcal{S}_4$, defined in (5)
$\mathcal{S}_1$	Positive orthant, defined in (6)
$\mathcal{S}_2$	Set describing the maximum number of active dictionary elements, defined in (7) with:
	$\mathcal{G}_{\boldsymbol{X}}$ Set indicating the columns of $\boldsymbol{X}$ contributing to the magnetisation response
	$\xi$ Minimum voxel proton density
	$\mathcal{D}_{\boldsymbol{X}}$ Set giving the voxels with significant contribution to the magnetisation response
	K Maximum number of active dictionary elements
$\mathcal{S}_{3}$	Set describing the constraint on the distinct active dictionary elements, defined in (8) with:
	$\upsilon$ Constant used to define neighbourhoods for the dictionary elements in the parameter space
	$\mathcal{N}_{\upsilon}$ Set describing the neighbour dictionary elements for given dictionary element d
$\mathcal{S}_{4}$	Set giving the minimum number of pure voxels per active dictionary element, defined in (10) with:
	$\mathcal{V}_{\boldsymbol{X}}$ Set describing the pure voxels in $\boldsymbol{X}$
	$\kappa$ Minimum number of pure voxels per active dictionary element
$\mu$	Step size in algorithm 1
$\boldsymbol{U}\in \mathbb{R}^{N \times \overline{T}}$	Sub-matrix of $\boldsymbol{X}$ (see equation (13))
$\mathbf{\boldsymbol{\Delta}}\in \mathbb{C}^{\overline{T} \times L}$	Sub-matrix of $\boldsymbol{\Phi}$ (see equation (14))
$\mathcal{Z}$	Operator from $\mathbb{R}_+^{N\times \overline{T}}$ to $\mathbb{R}_+^{N \times D}$ that maps a matrix $\boldsymbol{U}$ to the correspronding $\boldsymbol{X}$
$\boldsymbol{\lambda}\in \mathbb{C}^{N}$	Vetor used to compensate for the complex phase errors in the model
$\gamma$	Tolerance parameter for a pure voxels (see (22))
$\tau, \tau_k, \tau_\upsilon, \tau_\kappa$	Tolerance parameters for the initialisation process (see section 4)

In the context of MRF, the objective is to estimate the parameters of each voxel in the imaged volume from degraded undersampled measurements. Let $\boldsymbol{Y}\in\mathbb{C}^{Q\times L\times C}$ be the measurement matrix, where L is the excitation sequence length, C is the number of coils and Q is the number of measurements at each excitation and each coil. Let $\boldsymbol{M}\in \mathbb{C}^{N\times L}$ be the response of the imaged volume of interest with N voxels. For every $(l,c) \in \{1,\ldots,L\}\times \{1,\ldots,C\}$, the corresponding observation $\boldsymbol{Y}_{:,l,c} \in \mathbb{C}^{Q}$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{Y}_{:,l,c}=\boldsymbol{\Omega}_{:,:,l}\mathbf{F}\mathbf{S}_{:,:,c}\boldsymbol{M}_{:,l}+\boldsymbol{\eta}_{:,l,c}, \label{eq_mrfmodel} \nonumber \end{align} \tag{ 1 }$$

where $\boldsymbol{\Omega}\in\{1,0\}^{Q\times N\times L}$ is the concatenation of L selection matrices, $\mathbf{F}\in \mathbb{C}^{N\times N}$ is the 2D discrete Fourier transform, $\mathbf{S}\in\mathbb{C}^{N\times N\times C}$ is the concatenation of C spatial sensitivity coil diagonal matrices, and $\newcommand{\e}{{\rm e}} \boldsymbol{\eta}\in \mathbb{C}^{Q\times L \times C}$ is a realisation of a random i.i.d. Gaussian noise. Let ${h:\mathbb{C}^{N\times L}\rightarrow \mathbb{C}^{Q\times L \times C}}$ be the linear mapping defining the complete acquisition process such that $\newcommand{\e}{{\rm e}} \boldsymbol{Y}=h\left(\boldsymbol{M}\right) + \boldsymbol{\eta}$.

For each voxel $n\in \{1, \ldots, N\}$, the magnetisation response $\boldsymbol{M}_{n,:}$ is modelled through the smooth non-linear mapping $B:\mathcal{M} \rightarrow \mathbb{C}^{1\times L}$ (commonly Bloch equations or extended phase graphs (EPG) model [25]) scaled by the unknown proton density $\rho_n\in \mathbb{R}_+$, $\boldsymbol{M}_{n,:}=\rho_n B(\hat{\boldsymbol{\theta}}_{n,:},\boldsymbol{\Gamma})$, where $\boldsymbol{\Gamma}\in \mathbb{R}^{A\times 1}$ represents the concatenation of A known acquisition parameters (e.g. flip angles $\alpha$, repetition times TR) chosen such that $\boldsymbol{M}_{n,:}$ is only sensitive to the P parameters $\hat{\boldsymbol{\theta}}_{n,:}\in \mathcal{M}$ under investigation, where $\mathcal{M}\subset \mathbb{R}^{1 \times P}$ denotes the subset of feasible parameters. In the remainder, we fix P  =  2 and choose $\mathcal{M}$ corresponding to T₁ and T₂ relaxation times.

The model described in the previous section considers that each voxel contains at most one element. PV voxels are introduced due to the spatial discretisation in the acquisition process. The magnetisation sequence can be described as $\boldsymbol{M} = \boldsymbol{X}\boldsymbol{\Phi}$, where $\boldsymbol{X} \in \mathbb{R}_+^{N\times D}$ is a sparse mixing matrix (each line of $\boldsymbol{X}$ represents the proton densities associated with a specific voxel, and would contain more than a nonzero value only for voxels with partial volumes), and $\boldsymbol{\Phi}\in\mathbb{C}^{D\times L}$ is the over-complete dictionary of fingerprints, introduced in [2], as a discrete sampling of the low dimensional manifold B. $\boldsymbol{\Phi}$ is constructed from D samples of $\mathcal{M}$, stored in a matrix $\boldsymbol{\theta}\in \mathbb{R}^{D\times P}$. Due to the smoothness of B, $\boldsymbol{\Phi}$ is highly coherent. Consequently, the estimation of $\boldsymbol{X}$ from highly undersampled noisy data is expected to fail without additional priors. Leveraging CS theory [3, 26–28], the sparsest matrix $\boldsymbol{X}$, fitting the measurement model, can be found by solving:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \underset{\boldsymbol{X} \in \mathbb{R}_{+}^{N\times D}}{{{\rm minimise}}}\,\, \|\boldsymbol{X}\|_0 {{\rm ~subject~to~}} \|\boldsymbol{Y}-h(\boldsymbol{X}\boldsymbol{\Phi})\|_2 \leqslant \epsilon, \label{eq_nonconvexproblem} \nonumber \end{align} \tag{ 2 }$$

where $\newcommand{\e}{{\rm e}} \epsilon>0$ is a bound chosen according to the acquisition noise level. Since this function is non-convex and non-differentiable, problem (2) is difficult to solve in practice, in particular in the context of high dimensional problems (usually, $D \sim 10^6$ and $L \sim 10^3$). Greedy methods that deal directly with the $\newcommand{\e}{{\rm e}} \ell_0$ pseudo-norm, such as [29, 30], rely on a fixed sparsity level. Besides the computational burden of dealing with the huge dictionary $\boldsymbol{\Phi}$, these methods are not suited for partial volume estimations. The main reason is that the problem is intrinsically ill-posed, resulting in noise sensitive estimations without additional constraints. Note that the inclusion of any other constraint to these methods is not straightforward. The non-convexity of the $\newcommand{\e}{{\rm e}} \ell_0$ pseudo-norm is often relaxed by the use of the $\newcommand{\e}{{\rm e}} \ell_1$-norm [31]. Nevertheless, $\boldsymbol{\Phi}$ being highly coherent, this convex relaxation cannot be used to correctly estimate the coefficients of $\boldsymbol{X}$ [17, 32]. To overcome these difficulties, similarly to the BLIP approach, we propose to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \underset{\boldsymbol{M}\in \mathcal{B}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right)}{{{\rm minimise}}}\,\, \frac{1}{2}\|\boldsymbol{Y}-h(\boldsymbol{M})\|_2^2 \label{eq_pvablipproblem} \nonumber \end{align} \tag{ 3 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{B}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right) = \left\{\boldsymbol{M}\in \mathbb{C}^{N\times L} \,|\, \boldsymbol{M}=\boldsymbol{X} \boldsymbol{\Phi}{{\rm ~with~}} \boldsymbol{X}\in \mathcal{S}_+ \right\},\label{eq_setbs+} \nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{S}_+ =\overset{4}{\underset{s=1}{\cap}} \mathcal{S}_s, \label{eq_sets+} \nonumber \end{align} \tag{ 5 }$$

and, for every $s\in \{1, \ldots,4\}$, $\mathcal{S}_s$ is a closed non-empty subset of $\mathbb{R}^{N \times D}$ used to impose feasibility constraints on $\boldsymbol{X}$. These sets are defined below.

2.3.1. Positivity constraint.

Since the proton densities of the imaged volume must be non-negative, we can restrict our solution to be in the positive orthant:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{S}_1=\mathbb{R}_+^{N\times D}. \label{eq_positivity} \nonumber \end{align} \tag{ 6 }$$

2.3.2. Constraint on the number of tissues.

Commonly MRF aims to obtain quantitative values of a small set of tissues. In practice, only $T\ll D$ elements of the dictionary $\boldsymbol{\Phi}$ are necessary to characterise $\boldsymbol{M}$. While T is unknown, we have a reasonable estimate for it. We propose to introduce a loose upper bound K, such that $T \leqslant K \leqslant D$, to limit the number of active dictionary elements. Let us define a set $\mathcal{D}_{\boldsymbol X}$ that is formed by the column indices of $\boldsymbol{X}$ with non-zero coefficients. To avoid noisy voxels, only rows with proton density greater than $\xi>0$ (chosen according to the noise level) will be considered. Formally, this set is defined as $\newcommand{\e}{{\rm e}} {\mathcal{D}_{\boldsymbol X}=\left\{d \in \{1,\ldots,D\} \, | \, (\exists n \in \mathcal{G}_{\boldsymbol X}) \; \boldsymbol{X}_{n,d} \neq 0 \right\}}$, where $\mathcal{G}_{\boldsymbol X}=\{n\in\{1,\ldots,N\} \, | \, \|\boldsymbol{X}_{n,:}\|_1>\xi\}$. The set $\mathcal{D}_{\boldsymbol X}$ indicates the columns of $\boldsymbol{X}$ contributing to the magnetisation sequence. We can limit the number of used elements of the dictionary by upper bounding the cardinality of this set by K:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{S}_2=\left\{\boldsymbol{X}\in \mathbb{R}^{N \times D} \, | \, {{\rm Card}}\!\left(\mathcal{D}_{\boldsymbol X}\right) \leqslant K \right\}. \label{eq_K} \nonumber \end{align} \tag{ 7 }$$

2.3.3. Constraint on the manifold neighbourhoods.

The tissues of interest are unique and need to be sufficiently different to be distinguished. To incorporate this prior information in the reconstruction process, we define the neighbour set associated to each element $d\in \{1,\ldots,D\}$ of the dictionary as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\mathcal{N}_v(d)=\{d'\in \{1,\ldots, D\} \backslash \{d\} |\nonumber \\ &\left(\forall p=\{1,\ldots,P\}\right) |\boldsymbol{\theta}_{d',p}-\boldsymbol{\theta}_{d,p}|<\upsilon\boldsymbol{\theta}_{d,p} \}, \label{eq_neighbours} \nonumber \end{align} \tag{ 8 }$$

where $\upsilon>0$. We define a set of all possible $\boldsymbol{X}$ such that, the parameters of each element in $\mathcal{D}_{\boldsymbol X}$ are sufficiently far from each other. Precisely, we constrict all the neighbour columns of each element in $\mathcal{D}_{\boldsymbol X}$ to be the null element $\boldsymbol{0}$ of $\mathbb{R}^N$:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{S}_3=\big\{\boldsymbol{X}\in \mathbb{R}^{N\times D} \, | \, \big(\forall d' \in \underset{d \in \mathcal{D}_{\boldsymbol X}}{\cup} \mathcal{N}_v(d) \big) \, \boldsymbol{X}_{:,d'}= \boldsymbol{0} \big\}. \nonumber \end{align} \tag{ 9 }$$

2.3.4. Constraint on the pure voxels.

Due to the additive noise in model (1), some elements of $\boldsymbol{X}$ corresponding to non-used dictionary elements take non-zero values. In order to avoid these noisy elements in the reconstructions, we impose that at least $\kappa>0$ rows (i.e. voxels) of $\boldsymbol{X}$ contain only one non-zero value for each active column of $\boldsymbol{X}$. These rows identify the pure voxels. This constraint can be formulated as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{S}_4=\{\boldsymbol{X}\in \mathbb{R}^{N\times D} \, | \, \left(\forall d \in \mathcal{D}_{\boldsymbol{X}}\right) \, \|\left(\boldsymbol{X}_{n,d}\right)_{n\in \mathcal{V}_{\boldsymbol X}}\|_0\geqslant \kappa \} \label{eq_pure} \nonumber \end{align} \tag{ 10 }$$

where $\mathcal{V}_{\boldsymbol X}=\left\{n\in \left\{1,\ldots,N\right\} \, | \, \|\boldsymbol{X}_{n,:}\|_0=1\right\}$.

To solve problem (3), we use an iterative projected gradient descent method [33]. At each iteration $i\in \mathbb{N}$, this method updates $\boldsymbol{M}^{(i+1)}$ by computing a gradient step followed by a projection step:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{algo:grad_proj_mod} \boldsymbol{M}^{(i+1)} = \mathcal{P}_{\mathcal{B}_{S_+}\left(\boldsymbol{\Phi}\right)}\left(\boldsymbol{M}^{(i)}-\mu h^\dagger\big(h (\boldsymbol{M}^{(i)}) -\boldsymbol{Y} \big) \right), \nonumber \end{align} \tag{ 11 }$$

where $\mu >0$. In [5], it is shown that choosing $\mu\approx N/Q$ is theoretically justifiable. However, in order to ensure the stability of the iterative projected gradient descent algorithm and accelerate convergence, in [5, 34] the authors proposed to choose $\mu$ using a backtracking method. In order to handle efficiently the constraint $\mathcal{B}_{S_+}\left(\boldsymbol{\Phi}\right)$, we propose to compute inexactly the projection onto this set in (11). The resulting method, named greedy approximate projection for MRF (GAP-MRF), is described in algorithm 1. It can be noticed that the GAP-MRF method and BLIP are solving similar problems, using the same algorithmic structure. In this context, as in [5], a condition on both L and the undersampling ratio N/Q might be derived for recovery guarantee. However, the investigation of such condition is beyond the scope of this article.

Algorithm 1. GAP-MRF global iterations.

1: Input: $\boldsymbol{Y} \in \mathbb{C}^{Q\times L \times C}$, $\zeta <1$, $\boldsymbol{M}^{(0)} \in \mathbb{C}^{N \times L}$

2: Iterations:

3: for $i=0, 1,\ldots$ do

4:   $\mu = 2N/Q$, $\nu = 0$

5:   while $\mu > \nu$ do

6:       $\mu = \mu/2$

7:       Gradient Step:

8:       $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\overline{\boldsymbol{M}}^{(i)}= \boldsymbol{M}^{(i)}-\mu h^\dagger\left(h\left(\boldsymbol{M}^{(i)}\right)-\boldsymbol{Y}\right)$

9:       Projection Step:

10:      $\boldsymbol{M}^{(i+1)} \approx \mathcal{P}_{\mathcal{B}_{S_+} \left(\boldsymbol{\Phi}\right)} \left(\overline{\boldsymbol{M}}^{(i)} \right)$

11:      Backtracking step

12:      $\nu =\zeta \frac{\|\boldsymbol{M}^{(i+1)}-\boldsymbol{M}^{(i)}\|_2^2}{\|h\left(\boldsymbol{M}^{(i+1)}-\boldsymbol{M}^{(i)}\right)\|_2^2}$

13:   end while

14: end for

For every $\overline{\boldsymbol{M}}\in \mathbb{C}^{N\times L}$, we have:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \mathbb{\mathcal{P}}_{\mathcal{B}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right)}\left(\,\overline{\boldsymbol{M}}\right) &=\mathop{\rm argmin}\limits_{\boldsymbol{M}\in\mathcal{B}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right)} \frac{1}{2} \|\boldsymbol{M}-\overline{\boldsymbol{M}}\|_2^2\nonumber \\ &=\mathop{\rm argmin}\limits_{\boldsymbol{M}=\boldsymbol{X}\boldsymbol{\Phi},\, \boldsymbol{X}\in \mathcal{S}_+} \frac{1}{2} \|\boldsymbol{X}\boldsymbol{\Phi}-\overline{\boldsymbol{M}}\|_2^2 \nonumber \\ &= \big(\mathop{\rm argmin}\limits_{\boldsymbol{X}\in \mathcal{S}_+} \frac{1}{2} \|\boldsymbol{X}\boldsymbol{\Phi}-\overline{\boldsymbol{M}}\|_2^2 \big) \boldsymbol{\Phi}. \label{eq_trueprojection} \nonumber \end{align} \tag{ 12 }$$

Note that $\mathcal{S}_2,\mathcal{S}_3$ and $\mathcal{S}_4$ can be handled through the definition of $\boldsymbol{\Phi}$. Let $\boldsymbol{M}=\boldsymbol{X}\boldsymbol{\Phi} \in \mathcal{B}_{\mathcal{S}_+}\! (\boldsymbol{\Phi})$ and $\overline{T}\in \{1,\ldots,K\}$ (K is the upper bound defined in (7)). Let $\boldsymbol{U}\in \mathbb{R}^{N\times \overline{T}}$ be a subpart of $\boldsymbol{X}$ with non-zero columns and $\boldsymbol{\Delta}\in\mathbb{C}^{\overline{T}\times L}$ the corresponding subpart of $\boldsymbol{\Phi}$ such that $\boldsymbol{M}=\boldsymbol{U}\boldsymbol{\Delta}$. Then we have

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Proj_XS_PhiS} \mathcal{P}_{\mathcal{B}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right)} (\overline{\boldsymbol{M}})= \big(\mathop{\rm argmin}\limits_{\boldsymbol{U}\in\mathbb{R}_+^{N\times \overline{T}}} \frac{1}{2} \|\boldsymbol{U} \boldsymbol{\Delta} - \overline{\boldsymbol{M}}\|_2^2 \big)\boldsymbol{\Delta}. \nonumber \end{align} \tag{ 13 }$$

In (13), the dictionary $\boldsymbol{\Delta}$ is defined as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\Delta}=\mathop{\rm argmin}\limits_{\overline{\boldsymbol{\Delta}}\in \mathcal{C}} \big(\min_{\overline{\boldsymbol{U}}\in \mathbb{R}_+^{N\times \overline{T}}} \frac{1}{2}\|\overline{\boldsymbol{U}}\,\overline{\boldsymbol{\Delta}}-\overline{\boldsymbol{M}}\|_2^2 \big), \label{eq_Delta} \nonumber \end{align} \tag{ 14 }$$

where $\mathcal{C}$ is the set given by

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{C} = \Big\{\boldsymbol{\Delta} \in \mathbb{C}^{\overline{T} \times L} \, | \, (\exists \boldsymbol{X}\in \mathcal{S}_+) \; \boldsymbol{X} = \mathcal{Z} \big(\overline{\boldsymbol{U}} \big) \nonumber \\ {{\rm with~}} \overline{\boldsymbol{U}} = \underset{\boldsymbol{U}\in\mathbb{R}_+^{N\times \overline{T}}}{{{\rm argmin}}} \frac{1}{2}\|\boldsymbol{U}\boldsymbol{\Delta} - \overline{\boldsymbol{M}} \|_2^2 \Big\} \label{eq_dictionaries} \nonumber \end{align} \tag{ 15 }$$

with $\mathcal{Z} \colon \mathbb{R}_+^{N\times \overline{T}} \to \mathbb{R}_+^{N \times D}$ defined such that $\mathcal{Z}(\boldsymbol{U})\boldsymbol{\Phi}=\boldsymbol{U}\boldsymbol{\Delta}$.

As mentioned earlier, $\boldsymbol{\Phi}$ is an over-complete dictionary which makes the exact projection practically impossible to compute. Recent advances in reconstruction methods have introduced neural networks to efficiently approximate the projection or proximal operators within model based iterative algorithms [35–38]. A major challenge with such methods is obtaining sufficient accurate training data. In consequence, they can only accelerate the techniques where a prior computational solution to provide ground truth already exists. To overcome this difficulty, we propose a greedy approach to approximate the projection by finding a reduced dictionary $\widetilde{\boldsymbol{\Delta}}\in \mathbb{C}^{\overline{T}\times L}$ and its corresponding mixing matrix $\widetilde{\boldsymbol{U}}\in \mathbb{C}^{N\times \overline{T}}$, with $\overline{T}\leqslant K$, such that $\boldsymbol{U}\boldsymbol{\Delta} \approx \widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{\Delta}}$. Then the projection in step 10 of algorithm 1 can be approximated as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \mathbb{\mathcal{P}}_{\mathcal{B}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right)} (\overline{\boldsymbol{M}})\approx \big(\mathop{\rm argmin}\limits_{\widetilde{\boldsymbol{U}} \in \mathbb{R}_+^{N\times \overline{T}}} \frac{1}{2} \| \widetilde{\boldsymbol{U}} \widetilde{\boldsymbol{\Delta}} - \overline{\boldsymbol{M}} \|_2^2 \big) \widetilde{\boldsymbol{\Delta}}. \label{eq_pvaprojection} \nonumber \end{align} \tag{ 16 }$$

As mentioned in [5], it is a common practice to allow the proton density to be complex-valued in order to absorb phase terms correcting for timing and coil sensitivity errors. We incorporate a vector $\boldsymbol{\lambda} \in \mathbb{C}^{N}$ to compensate for these errors. Let $\widetilde{\mathcal{B}}_{\mathcal{S}_+}(\boldsymbol{\Phi})$ be the set of magnetisation sequences of the form $\boldsymbol{M} = {{\rm Diag}} (\boldsymbol{\lambda})\boldsymbol{X} \boldsymbol{\Phi}$ such that $\boldsymbol{X}\in \mathcal{S}_+$ and $\boldsymbol{\lambda}\in\mathbb{C}^{N}$ satisfies $(\forall n \in \{1, \ldots, N\}) \, |\boldsymbol{\lambda}_{n}|=1$. The approximate projection with the phase compensation is given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbb{\mathcal{P}}_{\widetilde{\mathcal{B}}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right)} (\overline{\boldsymbol{M}}) \approx {{\rm Diag}}(\boldsymbol{\lambda}) \widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{\Delta}}, \nonumber \end{align} \tag{ 17 }$$

where $(\boldsymbol{\lambda},\widetilde{\boldsymbol{U}})$ are obtained by solving:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \underset{\boldsymbol{\lambda}\in \mathbb{C}^{N}, \widetilde{\boldsymbol{U}}\in\mathbb{R}^{N\times \overline{T}}_+}{{{{\rm minimise}}}} \frac{1}{2}\|{{\rm Diag}}(\boldsymbol{\lambda}) \widetilde{\boldsymbol{U}} \widetilde{\boldsymbol{\Delta}} - \overline{\boldsymbol{M}} \|_2^2 \nonumber \\ {{\rm subject~to~}} (\forall n \in \{1, \ldots, N\}) \, |\boldsymbol{\lambda}_{n}|=1. \label{eq_complexapvprojection} \nonumber \end{align} \tag{ 18 }$$

It is worth mentioning that in (16) and (18), all the rows of $\widetilde{\boldsymbol{U}}$ can be computed independently in parallel.

On the one hand, forward–backward based algorithms [7, 39, 40] can be used to solve problem (16) (in particular, in our simulations, we use the built-in Matlab function of non-negative least-squares, that is an implementation of [41]). On the other hand, to solve problem (18), to jointly estimate $\boldsymbol{\Lambda}$ and $\boldsymbol{U}$, block coordinate approaches must be considered (e.g. Gauss–Seidel approaches [18], alternating forward–backward methods [19–22]). Note that in comparison with the traditional MRF methods which densely sample the manifold, our approach reduces the memory requirements, by using the dictionary $\widetilde{\boldsymbol{\Delta}}$ containing at most K elements, without the inaccuracies related to the manifold discretisation.

The GAP-MRF algorithm takes advantage of the dictionary coherence and the constraints imposed on $\boldsymbol{X}$ (described in section 2.3) to approximate the projection onto $\mathcal{B}_{\mathcal{S}_+}\left(\boldsymbol{\Phi}\right)$ in line 10 of algorithm 1. As described in section 3.2, this projection can be approximated at each iteration $i\in \mathbb{N}$, by solving (16), which necessitates to estimate the dictionary $\widetilde{\boldsymbol{\Delta}}^{(i)}$. We propose to estimate it using a greedy approach, leveraging both the knowledge of $\overline{\boldsymbol{M}}^{(i)}$ and the properties of the sets $\mathcal{S}_{2}$, $\mathcal{S}_{3}$ and $\mathcal{S}_{4}$ (note that the constraint $\mathcal{S}_1$ is handled directly in (16)). The proposed approach is described in details in this section.

The process to obtain $\widetilde{\boldsymbol{\Delta}}^{(i)}$ consists in three main steps leveraging the set of pure voxels. The first step consists in approximating the parameters of the pure voxels ($\mathcal{S}_4$ constraint) using the projection onto the set $\mathcal{B}^+$ defined as:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle &\mathcal{B}^+= \big\{\boldsymbol{M}\in \mathbb{C}^{N\times L}\,|\, (\forall n\in\{1,\ldots,N\}) \; \boldsymbol{M}_{n,:} = \rho \mathbf{m} ,\nonumber \\ &\; {{\rm with~}} \rho\in\mathbb{R}_+ {{\rm ~and~}} \mathbf{m} \in B\left(\mathcal{M},\boldsymbol{\Gamma} \right) \big\}. \label{eq_blochset} \nonumber \end{align} \tag{ 19 }$$

The objective of the second step is to find K regions of interest ($\mathcal{S}_2$ constraint) of the manifold by exploiting its smoothness. Finally, in the third step, the parameters that are too close to each other are discarded ($\mathcal{S}_3$ constraint) by using a non-maximum suppression based method [42]. This method acts on the number of voxels that corresponds to each parameter and keeps only the elements which have enough pure voxels to satisfy the $\mathcal{S}_4$ constraint. This process is summarised in the dictionary estimation step on figure 2. The remaining blue blocks in the diagram are used to update the variables in the greedy approximate projection. More precisely, we compute the mixing matrix $\widetilde{\boldsymbol{U}}^{(i)}$ and the magnetisation sequence $\boldsymbol{M}^{(i+1)}$ using equation (16) with the resulting dictionary $\widetilde{\boldsymbol{\Delta}}^{(i)}$. Then, we update the pure voxel set $\mathcal{V}_{\boldsymbol{X}}$ using the mixing matrix $\widetilde{\boldsymbol{U}}^{(i)}$. Finally, the dictionary $\boldsymbol{\Phi}$ is refined by randomly sampling around the parameters $\widetilde{\boldsymbol{\theta}}^{(i)}$. The complete method is described in algorithm 2 and explained in the following paragraphs.

Zoom In Zoom Out Reset image size

Algorithm 2. Greedy approximate projection.

1: Input: $\overline{\boldsymbol{M}}^{(i)},\boldsymbol{\Phi}^{(i)},\boldsymbol{\theta}^{(i)},\mathcal{V}_{\boldsymbol X}^{(i)},\boldsymbol{\Sigma}^{(i)},K,\boldsymbol{\Gamma},\kappa,\upsilon,\gamma,\beta,\xi,n_s$

2: Dictionary Estimation:

3:   Projection onto $\mathcal{B}^+$

4:     for $n=1,2,...,N$ do

5:     $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\displaystyle\hat{d}_n=\underset{d}{\rm argmax} \,\,\, {{\rm real}} (\overline{\boldsymbol{M}}^{(i)}_{n,:} \boldsymbol{\Phi}^{\dagger(i)}_{d,:}) / \|\boldsymbol{\Phi}^{(i)}_{d,:}\|_2$

6:     $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\displaystyle \hat{\rho}_n=\max ({{\rm real}} (\overline{\boldsymbol{M}}^{(i)}_{n,:} \boldsymbol{\Phi}^{\dagger(i)}_{\hat{d}_n,:}) / \| \boldsymbol{\Phi}^{(i)}_{\hat{d}_n,:} \|_2^2,0)$

7:     $\displaystyle\hat{\boldsymbol{\theta}}_{n,:}=\boldsymbol{\theta}^{(i)}_{\hat{d}_n,:}$

8:     end for

9:   Clustering

10:     $\displaystyle \mathcal{I}=\{n\in\mathcal{V}^{(i)}_{\boldsymbol X} \,|\, \hat{\rho}_n>\xi\}$

11:     $\displaystyle [\boldsymbol{\theta}^{\mathcal{S}_1\cap\mathcal{S}_2},\boldsymbol{c} ] = {{\rm k\mbox{-}means}} (\hat{\boldsymbol{\theta}}_{\mathcal{I},:},K)$

12:  Non-Maximum Suppression

13:     $\displaystyle\widetilde{\boldsymbol{\theta}}^{(i)}={{\rm NonMaximumSuppression}} (\boldsymbol{\theta}^{\mathcal{S}_1\cap\mathcal{S}_2},\boldsymbol{c},\upsilon,\kappa)$

14:     $\displaystyle\widetilde{\boldsymbol{\Delta}}^{(i)}=B(\widetilde{\boldsymbol{\theta}}^{(i)},\boldsymbol{\Gamma})$

15: Approximate Projection onto $\displaystyle\mathcal{B}_{\mathcal{S}_+}$

16:     $\displaystyle\widetilde{\boldsymbol{U}}^{(i)}=\underset{\overline{\boldsymbol{U}}\in \mathbb{R}_+^{N\times T}}{{\rm argmin}} \,\, \frac{1}{2}\|\overline{\boldsymbol{U}}\widetilde{\boldsymbol{\Delta}}^{(i)}-\overline{\boldsymbol{M}}^{(i)}\|^2_2$

17:     $\displaystyle\boldsymbol{M}^{(i+1)}=\widetilde{\boldsymbol{U}}^{(i)}\widetilde{\boldsymbol{\Delta}}^{(i)}$

18: Pure Voxel Set Update

19:     $\mathcal{G}_{\boldsymbol{X}}^{(i)}=\{n\in \{1,\ldots,N\} \,|\, \|\widetilde{\boldsymbol{U}}^{(i)}_{n,:}\|_1>\xi \}$

20:     $\mathcal{V}_{\boldsymbol X}^{(i+1)}=\{n\in\mathcal{G}_{\boldsymbol{X}}^{(i)}\,|\, \max (\widetilde{\boldsymbol{U}}^{(i)}_{n,:})\geqslant \gamma \|\widetilde{\boldsymbol{U}}^{(i)}_{n,:}\|_1 \}$

21: Parameter Re-sampling

22:     $\displaystyle\boldsymbol{\theta}^{(i+1)}={{\rm ParameterReSampling}}(\widetilde{\boldsymbol{\theta}}^{(i)},\boldsymbol{\Sigma}^{(i)},n_s)$

23:     $\displaystyle\boldsymbol{\Phi}^{(i+1)}=B(\boldsymbol{\theta}^{(i+1)},\boldsymbol{\Gamma})$

24:     $\displaystyle\boldsymbol{\Sigma}^{(i+1)}=\boldsymbol{\Sigma}^{(i)}\beta$

25: Output: $\boldsymbol{\theta}^{(i+1)} ,\boldsymbol{\Phi}^{(i+1)},\boldsymbol{\Sigma}^{(i+1)},\mathcal{V}_{\boldsymbol X}^{(i+1)}$ and $\boldsymbol{M}^{(i+1)}$

3.3.1. Projection onto ${\mathcal{B}^+}$.

At iteration $i\in \mathbb{N}$, we have:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \displaystyle \mathcal{P}_{\mathcal{B}_{\mathcal{S}_+ (\boldsymbol{\Phi})}} \big(\overline{\boldsymbol{M}}_{\mathcal{V}_{\boldsymbol{X}}^{(i)},:} \big) = \mathcal{P}_{\mathcal{B}^+ (\boldsymbol{\Phi})} \big(\overline{\boldsymbol{M}}_{\mathcal{V}_{\boldsymbol{X}}^{(i)},:} \big), \label{eq_purevoxels} \nonumber \end{align} \tag{ 20 }$$

where $\mathcal{B}^+$ is the set defined in equation (19), and $\displaystyle \overline{\boldsymbol{M}}_{\mathcal{V}_{\boldsymbol{X}}^{(i)},:} = (\overline{\boldsymbol{M}}_{n,:})_{n \in \mathcal{V}_{\boldsymbol{X}}^{(i)}}$, $\mathcal{V}_{\boldsymbol{X}}^{(i)}$ corresponding to an estimate of the pure voxel positions in $\boldsymbol{X}^{(i)}$ at iteration i (the true set $\mathcal{V}_{\boldsymbol{X}}$ corresponding to the pure voxels of the original $\boldsymbol{X}$ being unknown). At the first iteration, we choose $\mathcal{V}_{\boldsymbol{X}}^{(0)} = \{1, \ldots,N \}$, and it is updated during the greedy process (see algorithm 2, step 20). Note that that the region containing pure voxel does not need to be known a priori (it is automatically estimated), and does not need to be large to be detected (depending on the noise).

From (20), we can estimate the parameters $\hat{\boldsymbol{\theta}}$ and the proton density $\hat{\rho}$ of the voxels in $\mathcal{V}_{\boldsymbol{X}}^{(i)}$ using the projection onto $\mathcal{B}^{+}$ with a dictionary $\boldsymbol{\Phi}^{(i)}$ (see steps 4–8 of algorithm 2). $\boldsymbol{\Phi}^{(i)}$ is an adaptive dictionary that is refined at each iteration to reduce the computational cost, the simulations suggest that the accuracy of the reconstructions is preserved. Since there are at least $\kappa$ pure voxels for each active element in $\boldsymbol{\Phi}$ and the value of the proton density is at least $\xi$, we expect that the voxel parameters in $\mathcal{V}^{(i)}_{\boldsymbol X}$ with $\hat{\rho}>\xi$ will form clusters around the true values of the dictionary elements, an example can be seen in figure 3 (Left).

Zoom In Zoom Out Reset image size

3.3.2. Clustering.

3.3.3. Non-maximum suppression.

3.3.4. Inexact projection onto ${\mathcal{B}_{\mathcal{S}_+}}$.

3.3.5. Pure voxel set update.

In order to avoid noisy voxels, we re-define the set $\mathcal{G}_{\boldsymbol X}$, introduced in section 2.3.2, for $\widetilde{\boldsymbol{U}}^{(i)}$. Note that $\newcommand{\bi}{\boldsymbol} \mathcal{Z}\big(\widetilde{\boldsymbol{U}}^{(i)}\big)$ is a matrix of the size of $\boldsymbol{X}$ filling the missing values of $\widetilde{\boldsymbol{U}}^{(i)}$ with zeros, and thus we can re-define the set $\mathcal{G}_{\boldsymbol{X}}^{(i)}$ in terms of $\widetilde{\boldsymbol{U}}^{(i)}$ as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{G}_{\boldsymbol{X}}^{(i)}= \{n\in \{1,\ldots,N\} \,|\, \|\widetilde{\boldsymbol{U}}^{(i)}_{n,:}\|_1>\xi \}. \nonumber \end{align} \tag{ 21 }$$

Then, we update the pure voxel set as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{V}_{\boldsymbol X}^{(i+1)}= \{n\in\mathcal{G}_{\boldsymbol{X}}^{(i)}\,|\, \max(\widetilde{\boldsymbol{U}}^{(i)}_{n,:}) \geqslant \gamma \|\widetilde{\boldsymbol{U}}^{(i)}_{n,:}\|_1\}, \label{eq_gamma} \nonumber \end{align} \tag{ 22 }$$

where $0<\gamma<1$ is a relaxation factor used to compensate both for the noise and for the fact that the true dictionary elements are not guaranteed to be present. Note that the parameter $\gamma$ is defined as a proportion of the total proton density in the voxels, and it is used as a threshold to determine if a voxel is pure or not.

3.3.6. Parameter re-sampling.

The choice of $\xi$, setting the minimum proton density, is related to the background noise, the ideal $\xi$ is a value between the background noise and the signal in the volume of interest. If $\xi$ is too small, empty voxels will affect the clustering process. If it is too big, the tissue voxels will not be considered in the clustering process.

As mentioned before, the dictionary $\boldsymbol{\Phi}^{(i)}$ is updated through the iterations to reduce the complexity of the algorithm. We fix $\boldsymbol{\Phi}^{(0)}$ to all possible combinations of 20 values of T₁ and 20 values of T₂, equally spaced in $\mathcal{M}$.

Concerning the number of random samples n_s, on the one hand if we choose it too big, we increase the complexity of our pure voxel projection. On the other hand if we set n_s too small, more iterations will be needed to find the elements of interest. In all our simulations (simulated and in vivo data) we fix n_s  =  10.

For the diagonal elements of the covariance matrix (i.e. $\boldsymbol{\Sigma}_{1,1}^{(0)}$ and $\boldsymbol{\Sigma}_{2,2}^{(0)}$) associated to the resampling of the dictionary, if they are chosen too big, the parameter sampling will be far from the parameters of interest, increasing the number of iterations required to find them. If they are too small, the algorithm may not find the parameter of interest. $\boldsymbol{\Sigma}^{(0)}$ should be chosen based on the parameter separation of $\boldsymbol{\Phi}^{(0)}$. In all the reconstructions we fix $\boldsymbol{\Sigma}_{1,1}^{(0)}=40$ and $\boldsymbol{\Sigma}_{2,2}^{(0)}=10$.

Similarly, for the decreasing parameter $\beta$ of the covariance matrix (see step 24 in algorithm 2), if it is chosen too big, the algorithm will need more iterations to find the correct elements while if it is too small the algorithm may not explore the true parameters. We fix $\beta=0.9$ in the considered scenarios.

The choice of the pure voxel tolerance $\gamma$ is related to the noise and the accuracy of the dictionary during the iterations of the algorithm. If it is too big, the elements of interest could be eliminated through the iterations since pure voxels may be considered as PV voxels, if it is too small, the PV voxels may be considered as pure affecting the clustering process. We found in our simulations that $\gamma=0.85$ is a suitable choice.

The choice of the different parameters K, $\upsilon$ and $\kappa$ has been investigated during preliminary work. In particular, we observed a significant increase in the residual $\|\boldsymbol{Y}-h(\boldsymbol{M})\|_2$ when K is not sufficiently large. For $\upsilon$ and $\kappa$, we see a significant increase in the residual when they are chosen too large (i.e. merging proton density maps of the true tissues), and an increase of noisy proton density maps when they are chosen too low.

We propose to automatically choose K, $\upsilon$ and $\kappa$ by analysing the residual. Precisely, we choose a tolerance value on the residual, denoted by $\tau>0$. This value, indicates the minimum contribution of an element of the dictionary in the residual. If $\tau$ is chosen too big, our solution will contain noisy elements. While if it is chosen too small our elements of interest will be removed from the reconstruction.

The global GAP-MRF method, including the initialisation process, is described in algorithm 3. It describes the process to choose the parameters K, $\upsilon$ and $\kappa$. Firstly, the estimation of K is described in steps 2–10. Fixing all the other parameters, K is estimated by running multiple times the GAP-MRF iterations given in algorithm 1. We assume that we have a suitable estimate of K when the stopping criteria given in step 10 of algorithm 3 is reached. The same process is adopted for the estimation of $\upsilon$ described (steps 11–19) and $\kappa$ (steps 20–28). For these two estimates, we allow for a small tolerance ($\tau_\upsilon>0$ and $\tau_\kappa>0$, respectively), for robustness purposes. Note that each new run of algorithm 1 uses the previous estimated of $\boldsymbol{M}$, $\boldsymbol{\Delta}$ and $\boldsymbol{\theta}$, in order to accelerate the global method.

Algorithm 3. GAP-MRF global method.

1: Input: $\boldsymbol{Y}$, $\boldsymbol{\Gamma}$, $\boldsymbol{\Phi}$, $\boldsymbol{\theta}$, $\xi$, $\tau$, $\boldsymbol{\Sigma}^{(0)}_{1,1}=40$, $\boldsymbol{\Sigma}^{(0)}_{2,2}=10$, $\zeta=0.99$, $\mathcal{V}_{\boldsymbol X}^{(0)}=\{1,\ldots,N\}$, $\beta=0.9$, ns  =  10,

$\boldsymbol{M}^{(0)}=\boldsymbol{0}$, $(\tau_K, \tau_\upsilon, \tau_\kappa) = (10, 0.02, 10)$

2: Estimation of K:

3:   Input: $(\gamma, K, \upsilon, \kappa)=(0,0,0,0)$, $\widetilde{\boldsymbol{\Delta}}^{(0)}=\{\}$, $\widetilde{\boldsymbol{\theta}}^{(0)}=\{\}$, j   =  0.

4:   Do

5:    $\boldsymbol{\Phi}^{(0)}=\left[\boldsymbol{\Phi} ,\widetilde{\boldsymbol{\Delta}}^{(\,j)}\right]$, $\boldsymbol{\theta}^{(0)}=\left[(\boldsymbol{\theta}){}^{T} ,(\widetilde{\boldsymbol{\theta}}^{(\,j)}){}^T\right]^T$

6:    $K=K+ \tau_{K}$

7:    $\left[\boldsymbol{M}^{(\,j+1)},\widetilde{\boldsymbol{\Delta}}^{(\,j+1)},\widetilde{\boldsymbol{\theta}}^{(\,j+1)}\right] = {{\rm Algorithm~1}} (\boldsymbol{Y},\zeta,\boldsymbol{M}^{(\,j)})$

8:    j   =  j   +  1

9:   while $\|\boldsymbol{Y}-h(\boldsymbol{M}^{(\,j-1)})\|_2-\|\boldsymbol{Y}-h(\boldsymbol{M}^{(\,j)})\|_2 > \tau$.

10:   Output: $K^\star = K$, j   =  j   −  1

11: Estimation of $\upsilon$:

12:   Input: $(\gamma, K, \upsilon, \kappa)=(0.85, K^\star, 0, 0)$

13:   Do

14:    $\boldsymbol{\Phi}^{(0)}=\widetilde{\boldsymbol{\Delta}}^{(\,j)}$, $\boldsymbol{\theta}^{(0)}=\widetilde{\boldsymbol{\theta}}^{(\,j)}$

15:    $\upsilon=\upsilon+ \tau_{\upsilon}$

16:    $\left[\boldsymbol{M}^{(\,j+1)},\widetilde{\boldsymbol{\Delta}}^{(\,j+1)},\widetilde{\boldsymbol{\theta}}^{(\,j+1)}\right] = {{\rm Algorithm~1}}\left(\boldsymbol{Y},\zeta,\boldsymbol{M}^{(\,j)}\right)$

17:    j  =  j  +  1;

18:   while $\|\boldsymbol{Y}-h(\boldsymbol{M}^{(\,j)})\|_2-\|\boldsymbol{Y}-h(\boldsymbol{M}^{(\,j-1)})\|_2> \tau$.

19:    Output: $\upsilon^\star = \upsilon- 2\tau_{\upsilon}$, j   =  j   −  1

20: Estimation of $\kappa$:

21:   Input: $(\gamma, K, \upsilon, \kappa)= (0.85, K^\star, \upsilon^\star, 0)$

22:   Do

23:    $\boldsymbol{\Phi}^{(0)}=\widetilde{\boldsymbol{\Delta}}^{(\,j)}$, $\boldsymbol{\theta}^{(0)}=\widetilde{\boldsymbol{\theta}}^{(\,j)}$

24:    $\kappa=\kappa+ \tau_{\kappa}$

25:    $\left[\boldsymbol{M}^{(\,j+1)},\widetilde{\boldsymbol{\Delta}}^{(\,j+1)},\widetilde{\boldsymbol{\theta}}^{(\,j+1)}\right] = {{\rm Algorithm~1}}\left(\boldsymbol{Y},\zeta,\boldsymbol{M}^{(\,j)}\right)$

26:    j  =  j  +  1;

27:   while $\|\boldsymbol{Y}-h(\boldsymbol{M}^{(\,j)})\|_2-\|\boldsymbol{Y}-h(\boldsymbol{M}^{(\,j-1)})\|_2> \tau$.

28:    Output: $\kappa^\star = \kappa-2\tau_{\kappa}$

29: GAP-MRF Global Iterations:

30:   Input: $(\gamma, K, \upsilon, \kappa)= (0.85, K^\star, \upsilon^\star, \kappa^\star)$, $\boldsymbol{\Phi}^{(0)}=\widetilde{\boldsymbol{\Delta}}^{(\,j-1)}$, $\boldsymbol{\theta}^{(0)}=\widetilde{\boldsymbol{\theta}}^{(\,j-1)}$

31:   $\left[ \overline{\boldsymbol{M}}, \overline{\boldsymbol{\Delta}}, \overline{\boldsymbol{\theta}}\right] = {{\rm Algorithm~1}}\left(\boldsymbol{Y}, \zeta, \boldsymbol{M}^{(\,j-1)} \right)$

32:   $\displaystyle \overline{\boldsymbol{U}} = \mathop{\rm argmin}\limits_{\boldsymbol{U}\in\mathbb{R}_+^{N\times \overline{T}}} \frac{1}{2} \| \boldsymbol{U} \overline{\boldsymbol{\Delta}} - \overline{\boldsymbol{M}} \|_2^2$

33:   Output: $\overline{\boldsymbol{M}}$, $\overline{\boldsymbol{\Delta}}$, $\overline{\boldsymbol{\theta}}$ and $\overline{\boldsymbol{U}}$

We create a simulated phantom according to [46], with five tissues: adipose, WM, GM, muscle and cerebrospinal fluid (CSF). More precisely, to introduce the PVE, we use blocks of $2 \times 2$ voxels to form a lower resolution phantom containing PV voxels. The resulting volume is resized to $256\times 256$ voxels, and we generate the magnetisation sequence from this volume using the Bloch equations. All the reconstructions are performed with the same resolution. In the first column of figure 4 proton density maps and the voxel distribution of the simulated phantom are shown. Using this representation we can see the structure of the tissues of interest. Traditionally in qMRI, individual parameter maps are evaluated since only a tissue per voxel is considered but in a PV scenario this is not meaningful since several parameter maps would be needed and visually do not show the tissue structures. We also compute the dominant tissue (highest proton density in the voxel) parameter maps for a traditional evaluation. The phantom dominant tissue parameter maps can be seen in the first column of figure 5. Note that for the construction of the phantom, we only consider in-plane PV, while in reality through-plane PV and in-plane PV occurs. Both kind of PV are modelled the same way and should not make any difference in the reconstructions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to evaluate the algorithms, we use the signal-to-noise ratio (SNR in dB) defined as ${{\rm SNR}}(\boldsymbol{U}_{:,t},\widetilde{\boldsymbol{U}}_{:,t})=10\log\left(\sum_{n=1}^N \left(\boldsymbol{U}_{n,t}\right){}^2/\sum_{n=1}^N (\boldsymbol{U}_{n,t}-\widetilde{\boldsymbol{U}}_{n,t}){}^2 \right)$, where $t\in\{1,\ldots,T\}$ is the index of the evaluated tissue, $\boldsymbol{U}$ is the mixing matrix ground truth and $\widetilde{\boldsymbol{U}}$ is the estimation. Similarly for the magnetisation sequence SNR, we sum for all values in the matrix. To construct the matrix $\widetilde{\boldsymbol{U}}$, a tolerance of $15\%$ from the ground truth parameter values is used (i.e. for T₁  =  530 and T₂  =  77 milliseconds (ms) all the dictionary elements that fall for T₁ in the range of [450.5–609.5] ms and simultaneously for T₂ in the range of [65.45–88.55] ms are considered). In order to evaluate if the tissues are correctly identified, we define the success rate (SR) index as the proportion of voxels where the number of elements are correctly identified and its corresponding parameters fall within the $15\%$ of the true parameters. The same definition of SR is used for both pure and PV voxels (considering only the corresponding phantom voxels). Due to noise, there could be small values in $\widetilde{\boldsymbol{U}}$ that could significantly affect the SR. In consequence, we choose not to consider values that are smaller than 30, given that the range of the proton densities is from 80 to 400.

In this experiment, we investigate the behaviour of both the BLIP and the GAP-MRF algorithms while changing the input noise. We fix the magnetisation sequence length L  =  1000. The dictionary for BLIP is defined as in [5] with $D=16\,170$. The results correspond to an average (with standard deviation) over 10 runs of each choice of iSNR.

The results of the proton density maps are shown in figure 6 (Left). GAP-MRF significantly outperform BLIP when the iSNR is greater than 30 dB. We can notice that GAP-MRF estimates correctly the number of true atoms when the iSNR is 30 dB or greater. The reconstruction of adipose tissue is more affected by the noise since there are significantly less pure voxels of this tissue. GAP-MRF magnetisation sequence reconstruction is significantly more accurate than BLIP reconstruction, because BLIP does not consider the PVE and also because of the dictionary inaccuracy. GAP-MRF magnetisation sequence SNR has a linear behaviour with respect to the iSNR. In figure 6 (Center), the SR with respect to the iSNR can be seen. We can observe that the SR is significantly affected by the iSNR.

Zoom In Zoom Out Reset image size

The results for the dominant tissue parameter maps SNR can be seen in figure 6 (Right). GAP-MRF outperforms BLIP reconstructing the dominant tissue parameter maps. It is important to mention than GAP-MRF is more affected by noise because the linear combination of dictionary elements overfits the noise.

We show an example of the proton density maps for each tissue in figure 4 when the iSNR is 30 dB. By visual inspection, we can observe that the GAP-MRF method outperforms the BLIP method for PV reconstructions for moderate noise scenarios. The values of BLIP in table 2 are given in a range because multiple parameters were assigned to the corresponding ground truth tissue. On the contrary, GAP-MRF has a single value because only one value was assigned to the corresponding ground truth tissue. In this example, for BLIP and GAP-MRF respectively, the SNR values are as follows: 9.70 dB and 11.94 dB for Adipose, 9.14 dB and 19.52 dB for WM, 17.66 dB and 39.29 dB for Muscle, 8.31 dB and 31.87 dB for GM, 5.72 dB and 52.60 dB for CSF and for the magnetisation sequence 23.84 dB and 48.18 dB. The SR: 0.9944 and 0.9745 for pure voxels, and for PV voxels 0 and 0.9465. The GAP-MRF correctly estimates the manifold regions of interest. BLIP has a residual map formed by all the elements that are not sufficiently close to the true elements (see last row of figure 4). Note that the residual map is quite similar to the distribution of the PV voxels shown in figure 7, this shows that the parameter mismatch is due to the PVE. In the GAP-MRF reconstructions, the WM and adipose tissue are slightly mixed due to the noise since their parameters are close one to each other. By choosing a better $\boldsymbol{\Gamma}$ we can make the atoms of the dictionary more distant in the $\newcommand{\e}{{\rm e}} \ell_2$-norm sense, this would provide noise robustness to the reconstructions. The dominant tissue parameter maps are shown in figure 5. The T₁ and T₂ maps reconstructed by BLIP show a smooth transition from one tissue to another due to the partial volume. On the contrary, GAP-MRF reconstructions show abrupt transitions in the T₁ and T₂ maps delimiting the tissues. This is expected since each tissue is modelled with a unique set of parameters. We can observe that the dominant tissue proton density reconstruction of GAP-MRF is significantly affected by the noise. Nevertheless, thanks to the constraint $\mathcal{S}_+$ handled by the proposed method, the T₁ and T₂ parameter maps are accurate.

In this subsection, we compare the proposed GAP-MRF algorithm with the BLIP algorithm, for different number of excitation instances L. The iSNR is set to 50 dB. The dictionary for BLIP is defined as in [5] with $D=16\,170$. The results correspond to an average (with standard deviation) over 10 runs of each choice of L.

Figure 8 (Left) shows the evaluation of the proton density maps for each tissue (Adipose, WM, Muscle, GM, and CSF) and the magnetisation sequence. Note that GAP-MRF results are taken directly from the matrix $\widetilde{\boldsymbol{U}}$ without using any post-processing. We can observe that GAP-MRF outperforms BLIP in reconstructing $\boldsymbol{U}$. This can be explained by the fact that BLIP is restricted to the input dictionary, while our method estimates the dictionary. In addition, we can observe that the SNR values of magnetisation sequence reconstructed with BLIP slightly decreases while L increases (while it is not the case for the proton density maps). This is expected since the linear combination of short fingerprints are less distinctive, hence it is easier to approximate it with fingerprints of other elements (allowing BLIP to fit better PV voxels with other elements). This behaviour is not observed with the proposed GAP-MRF method for which accurate proton density map estimates result in accurate magnetisation sequence reconstructions.

Zoom In Zoom Out Reset image size

Figure 8 (Center) gives the SR for both pure and PV voxels. Since BLIP can only reconstruct one element per voxel, its SR for PV is always equal to 0. In a low noise scenario, GAP-MRF can identify the correct voxel elements even for short sequences as can be seen in figure 8 (Center), where the SR of GAP-MRF for pure and PV voxels is 1. An important remark is that due to PV, the dictionary sampling and the number of excitation instances, the BLIP algorithm can mis-reconstruct pure voxels even in a low noise scenario.

For the dominant tissue parameter maps in the low noise scenario, GAP-MRF outperforms BLIP as shown in figure 8 (Right). The proton density map of BLIP is affected by the PV since it is not able to distinguish between the voxel tissues. The T₁ and T₂ maps are affected by the PV voxels and the dictionary inaccuracies.

BLIP reconstructions show a variation on T₁ and T₂ for the same tissue while GAP-MRF reconstructions are accurate. The GAP-MRF has the additional advantage that it simultaneously estimates the manifold regions of interest, resulting in better reconstructions.

In this subsection, we show the results obtained with the proposed approach and BLIP method, considering a dataset from a GE HDx MRI system with an 8 channel receive only head RF coil (GE Medical Systems, Milwaukee, WI). The acquisition scheme uses a variable density spiral with 377 interleaves using FISP based $\alpha$ [48] and a constant TR  =  10ms. The excitation sequence length is L  =  1000. In this experiment, we have FOV  =  $22.5\times22.5$ cm² with a 5 mm slice thickness. The EPG model is used for the reconstructions with an inversion time (TI) of 18 ms and an echo time TE  =  1.902 ms. The scanned objects are the tubes 1, 5 and 9 of the EUROSPIN phantom. We reconstruct the parameter maps with two spatial resolutions: the first one at $180\times 180$ with an undersampling ratio of N/Q  =  44.8753, and the second one at $40\times 40$ with an undersampling ratio of N/Q  =  20.6869 to introduce the PV. Note that for the $40\times 40$ reconstruction, only the Fourier samples corresponding to the target resolution are used. Reconstructing for higher spatial resolution would introduce high frequency artefacts as shown in [49]. An acquisition without the tubes is performed to estimate $\sigma_{\boldsymbol{Y}}$ and compute a lower bound on the iSNR. More precisely, using the triangle inequality, since $\newcommand{\e}{{\rm e}} \| \boldsymbol{Y}\|_2 \geqslant \| \boldsymbol{\eta}\|_2$, we have $\newcommand{\e}{{\rm e}} {{\rm iSNR}}\geqslant 20\log\left((\|\boldsymbol{Y}\|_2-\|\boldsymbol{\eta}\|_2)/(\sqrt{QLC}\sigma_{\boldsymbol{Y}})\right)=64.73~{{\rm dB}}$, where $\boldsymbol{Y}$ corresponds to the measurements with the tubes, $\newcommand{\e}{{\rm e}} \boldsymbol{\eta}$ corresponds to the measurements without the tubes, and the value $\sigma_{\boldsymbol{Y}}$ is the standard deviation of $\newcommand{\e}{{\rm e}} \boldsymbol{\eta}$.

The box in red shows a PV voxel artificially created by reconstructing a lower resolution image. As predicted by the corresponding high resolution maps, this voxel is formed by a linear combination of the Tubes 1 and 9.

In figure 9, we show a comparison of the reconstructions with two different spatial resolutions. The T₁ and T₂ lower resolution maps of BLIP show a variation introduced by the PVE. A clear example of the PVE is the voxel in the red box where two tissues appear, the BLIP method shows a parameter mismatch. Note that the parameters predicted by BLIP suggest that the voxel contains the same substance as Tube 5, contrary to the true composition (Tubes 1 and 9). GAP-MRF reconstructions do not show this behaviour since we take the PV into account in the model. Note that GAP-MRF is more sensitive to noise as shown in the simulations, this may explain small artefacts in the proton density maps.

Zoom In Zoom Out Reset image size

The T₁ values in table 3 are in agreement with the values of the EUROSPIN phantom. The T₂ values are higher than expected. As seen in figure 9, BLIP results show the same increased T₂, suggesting that the errors may be related to the acquisition parameters specifically to the constant TR as shown in [50].

Table 3. Comparison between the parameters obtained with GAP-MRF and the EUROSPIN phantom values.

	Phantom values		$180\times 180$		$40\times 40$
	T1	T2	T1	T2	T1	T2
Tube 1	$200\pm 6$	$52\pm 1.6$	197.0	93.9	195.4	96.0
Tube 5	$450\pm 13.5$	$94\pm 2.8$	455.9	159.8	459.5	168.1
Tube 9	$754\pm 22.6$	$116\pm 3.5$	766.3	199.2	757.5	199.9

In figure 10, the normalised proton density maps reconstructed by GAP-MRF with two different resolutions can be seen. As highlighted by the red box, the PV voxel in the low resolution reconstruction has values different than 0 in the maps corresponding to Tube 1 and 9, this is in agreement with the high resolution maps.

Zoom In Zoom Out Reset image size

This dataset was acquired by self-experimentation on our team members (all experts in MRI). Since these experiments are not intended to be qualified as a clinical investigation, they do not require any formal IRB approval according to the German act on medical devices (Medical Device Act, MDA). The self-experiments were performed on a device that has already met the requirements of the assessment procedure of conformity, certifying its safety and functionality for the intended purpose (aka 'CE marking of MR scanner'). The experiments were neither invasive nor stressful, therefore, they fully comply with internal GE and German/EU regulations. The scanning for this dataset was performed on a GE HDx MRI system with an 8 channel receive only head RF coil (GE Medical Systems, Milwaukee, WI). The acquisition scheme uses a variable density spiral with 89 interleaves using FISP based $\alpha$ and TR as in [48]. The excitation sequence length is L  =  1000. In this experiment, we have FOV  =  $22.5\times22.5$ cm² and the spatial resolution is $180\times 180$ voxels, with a 5 mm slice thickness. The undersampling ratio is N/Q  =  89.53. The EPG model is used for the reconstructions with a TI of 18 ms and a TE of 2 ms. The reconstruction for BLIP and GAP-MRF was accelerated with the SVD compression in the time domain described in [13, 49] using 30 eigenvectors.

In figure 11, we can observe the resulting proton density maps provided by the GAP-MRF algorithm and the table 4 shows a comparison between the parameters reported in [51] and the parameters obtained by GAP-MRF. The WM, GM and Fat parameters obtained by GAP-MRF slightly differ from those reported for MRF sequences but the values are in agreement with the parameters of other qMRI methods reported in [51]. The muscle parameters are far from the expected values. This could be due to the small number of pure voxels that are not sufficient to accurately estimate the parameters. We believe that choosing better acquisition parameters $\boldsymbol{\Gamma}$ to make the elements of the dictionary more distant in the $\newcommand{\e}{{\rm e}} \ell_2$-norm sense can significantly improve the accuracy of the parameters. Also, inaccuracies in the model such as calibration or motion in the acquisition can produce artefacts in the reconstruction. In order to show the importance of the phase compensation in the real data, we present the results without phase compensation. As seen in figure 12, due to the phase errors there are voxels within the brain without proton density.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Comparison between the parameters obtain with GAP-MRF for the brain dataset with spiral sampling and the reported values in [51].

	Values reported in [51]		GAP-MRF
	T1	T2	T1	T2
WM	781 $\pm$ 61	65 $\pm$ 6	758.7	42.1
GM	1193 $\pm$ 65	109 $\pm$ 11	872.4	67.3
CSF			1658.5	799.8
Muscle	1100 $\pm$ 59	44 $\pm$ 9	1218.0	23.2
Fat	253 $\pm$ 42	68 $\pm$ 4	325.5	68.1

In figure 13, the T₁ and T₂ maps reconstructed by BLIP show a smooth transition from one tissue to another (similar to the simulated phantom). Moreover, the proton density map reconstructed by BLIP does not provide any information on the tissue distribution. On the contrary, GAP-MRF reconstructions show abrupt transitions in the T₁ and T₂ maps of the dominant tissues. In addition, the proton density map shows more structure than BLIP, but not all the tissue structures are appreciated compared to the normalised proton density maps. Note that the voxels with higher proton density values indicate the pure voxels, and the voxels with reduced values, which are observed at tissue interfaces, indicate the dominant tissue that occupying only a fraction of the voxel.

Zoom In Zoom Out Reset image size

The scanning for this dataset has been performed on a 3T GE MR750w scanner with a 12 channel receive only head RF coil (GE Medical Systems, Milwaukee, WI). The study was approved by the local ethics committee. The used acquisition scheme was 16-shot EPI-MRF on a healthy volunteer using a variable flip angle $\alpha$ ramp, ranging from $1^{\circ}$ to $70^{\circ}$. The excitation sequence length is L  =  500. The repetition time TR was set to 16 ms. In [50], it was shown to be as effective at estimating the MRF parameters but had better sensitivity than the FISP sequence in [48]. The acquisition bandwidth (BW)  =  5 kHz and the Field of View (FOV)  =  $22.5\times22.5$ cm². The spatial resolution is $128\times 128$ voxels, with a 5 mm slice thickness. The undersampling ratio is N/Q  =  16. The EPG model is used for the reconstructions with an inversion time (TI) of 18 ms and an echo time (TE) of 3.5 ms. The acquisition time for the slice was 9 s. A reference scan with null Gy gradient was performed for phase correction of EPI raw data.

Table 5. Comparison between the parameters obtain with GAP-MRF for the brain dataset with EPI sampling and the reported values in [51].

	Values reported in [51]		GAP-MRF
	T1	T2	T1	T2
WM	781 $\pm$ 61	65 $\pm$ 6	762.6	67.2
GM	1193 $\pm$ 65	109 $\pm$ 11	1116.6	107.1
CSF			2391.1	856.2

In figure 14, we can observe the resulting proton density maps provided by the GAP-MRF algorithm and the table 5 shows a comparison between the parameters reported in [51] for MRF FISP sequences and the parameters obtained by GAP-MRF. CSF values are not reported for the MRF FISP sequence. The WM parameters are similar to the ones reported in [51] and the GM T₁ is slightly lower than the reported one. We believe that the lack of pure voxels (due the spatial resolution) made the approach unable to find the other tissues.

Zoom In Zoom Out Reset image size

In figure 15, the T₁ and T₂ maps reconstructed by BLIP shows a smooth transition from one tissue to another. On the contrary, GAP-MRF reconstructions show abrupt transitions in the T₁ and T₂ maps of the dominant tissues.

Zoom In Zoom Out Reset image size