Low Rank Magnetic Resonance Fingerprinting

TL;DR

This paper introduces FLOR, a novel low-rank based method for Magnetic Resonance Fingerprinting that improves quantitative tissue parameter estimation accuracy in highly undersampled MRI data.

Contribution

The paper presents a new low-rank iterative approach for MRF that leverages the low rank property of the signal to enhance parameter map accuracy over existing methods.

Findings

FLOR outperforms other methods at 5% and 9% sampling ratios.

Experimental results show improved parameter accuracy.

Both retrospective and prospective experiments validate FLOR's effectiveness.

Abstract

Purpose: Magnetic Resonance Fingerprinting (MRF) is a relatively new approach that provides quantitative MRI measures using randomized acquisition. Extraction of physical quantitative tissue parameters is performed off-line, without the need of patient presence, based on acquisition with varying parameters and a dictionary generated according to the Bloch equation simulations. MRF uses hundreds of radio frequency (RF) excitation pulses for acquisition, and therefore a high undersampling ratio in the sampling domain (k-space) is required for reasonable scanning time. This undersampling causes spatial artifacts that hamper the ability to accurately estimate the tissue's quantitative values. In this work, we introduce a new approach for quantitative MRI using MRF, called magnetic resonance Fingerprinting with LOw Rank (FLOR). Methods: We exploit the low rank property of the concatenated temporal imaging contrasts, on top of the fact that the MRF signal is sparsely represented in the generated dictionary domain. We present an iterative scheme that consists of a gradient step followed by a low rank projection using the singular value decomposition. Results: Experimental results consist of retrospective sampling, that allows comparison to a well defined reference, and prospective sampling that shows the performance of FLOR for a real-data sampling scenario. Both experiments demonstrate improved parameter accuracy compared to other compressed-sensing and low-rank based methods for MRF at 5% and 9% sampling ratios, for the retrospective and prospective experiments, respectively. Conclusions: We have shown through retrospective and prospective experiments that by exploiting the low rank nature of the MRF signal, FLOR recovers the MRF temporal undersampled images and provides more accurate parameter maps compared to previous iterative methods.

Equations

$$\mathbf{Y}=[F_u\{{\mathbf{X}}_{:,1}\},...,F_u\{{\mathbf{X}}_{:,L}\}]+\mathbf{H}$$

$${\mathbf{\Theta }}_1^j=[T{1}^j,T{2}^j,P{D}^j],1\le j\le N.$$

$${\mathbf{\Theta }}_2^i=[T{R}^i,T{E}^i,F{A}^i{]}^T,1\le i\le L$$

$$\underset{\mathbf{X},\mathbf{R}}{\mathrm{minimize}}$$

$$\text{rank}(\mathbf{X})\le r$$

$$\mathbf{X}={\mathbf{R}}_1\mathbf{D}$$

$${\mathbf{X}}^{n+1},{\mathbf{R}}_1^{n+1},{\mathbf{Z}}^{n+1}$$

$$+{\eta }_2{\Vert {\mathbf{W}}^n-\mathbf{X}+\mathbf{Z}\Vert }_F^2$$

$${\mathbf{Q}}^{n+1}$$

$${\mathbf{W}}^{n+1}$$

$$\underset{\mathbf{X}\in \mathbb{D}}{\mathrm{minimize}}$$

$${\text{NMSE}}_i=\frac{||{\theta }_i-{\hat{\theta }}_i|{|}_F^2}{||{\theta }_i-\frac{1}{N}\sum _j{\theta }_i^j|{|}_F^2}$$

$$\underset{\mathbf{X}\in \mathbb{D}}{\mathrm{argmin}}\frac{1}{2}\underset{i}{\Sigma }{\Vert {\mathbf{Y}}_{:,i}-F_u\{{\mathbf{X}}_{:,i}\}\Vert }_2^2+\lambda ||\mathbf{X}|{|}_{\ast }$$

$$\underset{\mathbf{X}\in \mathbb{D}}{\mathrm{argmin}}\{{\Sigma }_i^m{F}_i(\mathbf{X})\}$$

$${\mathbf{Z}}^k={\text{P}}_{\mathbb{D}}({\mathbf{X}}^k-{\mu }_k{g}_{i_k})$$

$${\mathbf{X}}^{k+1}=\underset{\mathbf{X}\in \mathbb{D}}{\mathrm{argmin}}{f}_{i_k}(\mathbf{X})+\frac{1}{2{\mu }_k}\Vert \mathbf{X}-{\mathbf{Z}}^k{\Vert }_F^2$$

$${\text{P}}_{\mathbb{D}}(\mathbf{X})$$

$$h(\mathbf{X})$$

$$f(\mathbf{X})$$

$$\partial h(\mathbf{X})$$

$${\text{P}}_{\mathbb{D}}(\mathbf{X})$$

$${\text{SVT}}_{{\mu }_k\lambda }({\mathbf{Z}}^k)={\mathbf{U}}_r[{\mathbf{\Sigma }}_r-{\mu }_k\lambda \mathbf{I}{]}_{+}{\mathbf{V}}_r^H.$$

$${\mathbf{X}}^{k+1}={\mathbf{U}}_r[{\mathbf{\Sigma }}_r-{\mu }_k\lambda \mathbf{I}{]}_{+}{\mathbf{V}}_r^H.$$

$${\mathbf{Z}}^k=({\mathbf{X}}^k-{\mu }_k{F}_u^H\{\mathbf{Y}-F_u\{{\mathbf{X}}^k\}\})\mathbf{P}$$

$${\mathbf{X}}^{k+1}={\mathbf{U}}_r[{\mathbf{\Sigma }}_r-{\mu }_k\lambda \mathbf{I}{]}_{+}{\mathbf{V}}_r^H.$$

Full text

Methods: We exploit the low rank property of the concatenated temporal imaging contrasts, on top of the fact that the MRF signal is sparsely represented in the generated dictionary domain. We present an iterative scheme that consists of a gradient step followed by a low rank projection using the singular value decomposition.

Results: Experimental results consist of retrospective sampling, that allows comparison to a well defined reference, and prospective sampling that shows the performance of FLOR for a real-data sampling scenario. Both experiments demonstrate improved parameter accuracy compared to other compressed-sensing and low-rank based methods for MRF at 5% and 9% sampling ratios, for the retrospective and prospective experiments, respectively.

Conclusions: We have shown through retrospective and prospective experiments that by exploiting the low rank nature of the MRF signal, FLOR recovers the MRF temporal undersampled images and provides more accurate parameter maps compared to previous iterative methods.

MRF, Low rank, Compressed Sensing, QMRI

I Introduction

Quantitative Magnetic Resonance Imaging (QMRI) is widely used to measure tissue’s intrinsic spin parameters such as the spin-lattice magnetic relaxation time (T1) and the spin-spin magnetic relaxation time (T2) Ehses et al. [2013]. Since tissue relaxation times vary in disease, QMRI enables the diagnosis of different pathologies, including multiple sclerosis (MS), Alzheimer, Parkinson, epilepsy and cancer Baselice et al. [2016], Antonini et al. [1993], Baselice et al. [2014], Haley et al. [2004], Mariappan et al. [1988], Roebuck et al. [2009]. In addition, the knowledge of tissue relaxation times allows generation of many clinical MR imaging contrasts (such as FLAIR and STIR) off-line, and saves a significant amount of scanning time.

Despite the advantages of QMRI, clinical MRI today mainly consists of weighted images. Values in weighted MR imaging are given in arbitrary units, since the signal strength is influenced by both intrinsic parameters (such as relaxation times and concentration of hydrogen atoms) and non-intrinsic ones. Non-intrinsic parameters include transmit and receive coils sensitivities, patient position in the scanner, vendor based scanner specific parameters, and local temperature. Weighted MRI images therefore lack quantitative information and as a result, different materials may exhibit similar or identical gray level values. In addition, weighted MRI contrast values vary between different follow-up scans of the same patient. This fact may impair disease monitoring, if based solely on those images. To date, weighted MRI scans are more common than QMRI in the clinic, due to the extremely long times often associated with QMRI using conventional techniques Liu et al. [1989], Homer and Beevers [1985], Crooijmans et al. [2011].

A plethora of methods have been proposed for QMRI. Earlier approaches are based on a series of spin echo (SE) or inversion recovery (IR) images with varying repetition times (TR) and echo times (TE) to evaluate each magnetic parameter (T1 or T2) separately. After acquisition, the curve of intensities for each pixel is matched to the expected magnetic signal, representing the appropriate magnetic tissue parameters Liu et al. [1989]. Accelerated methods for QMRI consist of the popular driven equilibrium single pulse observation of T1 (DESPOT1) Homer and Beevers [1985] or T2 (DESPOT2) Crooijmans et al. [2011] and the IR TrueFISP for simultaneous recovery of T1 and T2 quantitative maps Schmitt et al. [2004], Ehses et al. [2013]. Both techniques do not require long waiting times between excitations to reach an equilibrium state, and therefore they are significantly faster. Later works shortened the acquisition time required by those methods by under-sampling the data in both spatial and temporal domains Doneva et al. [2010], Zhao et al. [2015], Petzschner et al. [2011], Huang et al. [2012], Velikina et al. [2013]. However, obtaining accurate and high resolution QMRI in a reasonable clinical scanning time is still very challenging.

An approach for QMRI called magnetic resonance fingerprinting (MRF) has drawn increased attention in the last few years Ma et al. [2013]. MRF uses pseudo-randomized acquisitions to generate many different imaging contrasts, acquired at a high under-sampling ratio. It exploits the different acquisition parameters over time to produce a temporal signature, a “fingerprint”, for each material under investigation. By matching this unique signature to a pre-generated set of simulated patterns, the quantitative parameters can be extracted off-line. This approach saves valuable scan time compared to previous methods for accelerated QMRI, demonstrating promising efficient and reliable results.

MRF utilizes the fact that each tissue responds differently to a given quasi-random pulse sequence. By varying the acquisition parameters (e.g. repetition time (TR), echo time (TE), and radio frequency flip angle (FA)), unique signals are generated from different tissues. After acquisition, a pattern recognition algorithm is used to match the acquired signal from each voxel to an entry from a dictionary of possible tissue candidates. The dictionary entries are created by simulating the tissue’s response to the sequence for a range of T1 and T2 parameter values, using the Bloch equations. The resulting dictionary contains the temporal signatures of various simulated materials, given the pseudo-random pulse sequence. The quantitative parameters, such as the tissue’s T1 and T2 relaxation times, can be retrieved from the data by matching the signature acquired to the most correlated entry in the dictionary.

In MRI, data is acquired in the Fourier domain of the spatial image (a.k.a. k-space). The acquisition time of a high resolution, single contrast 3D MRI lasts a substantial amount of time. Since MRF is based on rapid acquisition of hundreds of different contrasts, severe under-sampling is performed in k-space to obtain the temporal resolution required for MRF. Figure 1 demonstrates the effect of fully sampled versus under-sampled data, acquired with spiral trajectories and recovered using the inverse non uniform fast Fourier transform (NUFFT) Fessler and Sutton [2003]. It can be seen that the under-sampled data is blurred and introduces aliasing artifacts. Figure 2 illustrates the noise and under-sampling artifacts of a representative brain voxel intensity as function of time, where the data is acquired with an MRF sequence based on fast imaging with steady state precession (FISP) Jiang et al. [2015]. Clearly, under-sampling also introduces a substantial level of noise in the time domain. In addition, MRF uses a dictionary with discrete values, while QMRI values are continuous. This leads to quantization error, depending on the values represented in the dictionary.

While in the original MRF paper Ma et al. [2013] these imaging artifacts are not handled explicitly, recent works have implemented advanced reconstruction techniques to overcome under-sampling artifacts. Approaches based on exploiting the sparsity of the signal in some transform domain in a compressed sensing (CS) Eldar and Kutyniok [2012], Eldar [2015] framework are examined by Davies et al. Davies et al. [2014] and Wang et al. Wang et al. [2016]. Zhao et al. Zhao et al. [2016a] formulated MRF as a maximum likelihood (ML) problem, and developed an iterative reconstruction approach for each time point. While these techniques showed improved results compared to the original MRF method, they do not exploit the temporal similarity between adjacent time-points, which is intrinsic to the dynamic acquisition used in MRF.

A common approach to exploit redundancy exists in dynamic MRI, based on modeling the acquired data as low-rank. This modelling was successfully applied for various dynamic MRI applications, such as cardiac imaging Feng et al. [2013] and functional MRI Chiew et al. [2015]. In the context of MRF, early works use low-rank MRF to compress the dictionary for faster reconstruction McGivney et al. [2014], Cline et al. [2016]. This saves reconstruction time, but does not necessarily improve the quality of the reconstructed maps or the acquisition time. The first introduction of a low-rank constraint for improved reconstruction in MRF was proposed by Zhao et al. Zhao [2015], Zhao et al. [2016b] followed by a sub-space constrained low-rank approach introduced by us Mazor et al. [2016]. Extensions of these ideas include adding a sparse term to the low-rank-based reconstruction Liao et al. [2016] (a.k.a robust PCA Candès et al. [2011]) and representing the data as low-rank in the k-space domain Doneva et al. , Doneva et al. [2017]. Recently, a few approaches that utilize prior knowledge of the dictionary together with a low-rank constraint have been published. Zhao et al. Zhao et al. [2017] presented an efficient algorithm that performs a singular value decomposition (SVD) on the dictionary and embeds the right singular vectors in the solution, to obtain better estimation of the temporal signatures. A similar approach was presented by Assländer et al. Assländer et al. [2017], who embed the left singular vectors in the solution. These methods show that exploiting the redundancy via a low-rank based solution improves the results compared to a sparsity approach. However, the obtained reconstructed maps still suffer from quantization error, due to the nature of a matched-filter based solution that matches a single dictionary atom to a single pixel. In addition, most of these methods are based on a fixed rank, set in advance, which may be difficult to determine in advance.

In this paper we extend our initial idea presented in our conference paper Mazor et al. [2016] and enforce a low-rank constraint in the image domain together with constraining the solution to the dictionary subspace. In particular, we exploit the low-rank property of the temporal MRF domain, via an iterative scheme that consists of a gradient step followed by a projection onto the subspace spanned by the dictionary elements in order to constrain the structure of the tissue behaviour simulated in the dictionary. The estimated images are then decomposed using SVD and the singular values are soft-thresholded to reduce the nuclear norm value in every step. Our approach, called magnetic resonance Fingerprinting with LOw Rank constraint (FLOR), incorporates three main advantages that were only partially introduced in previous work:

• •

FLOR formulates the problem as a convex problem. The solution is then rigorously developed based on the incremental subgradient proximal method Nedic and Bertsekas [2001]. This technique is known to convergence to the global minimum, regardless of the initial starting point.

• •

FLOR is based on a nuclear-norm solution, and does not require fixing the rank in advance. This leads to a solution that adapts the rank according to the nature of the specific dataset.

• •

The subspace constraint in FLOR is not limited to dictionary items, but rather allows a solution that is spanned by the dictionary. This enables better reconstruction of the temporal imaging contrasts. It also allows generation of quantitative parameters that do not necessarily exist in the simulated dictionary, thereby reducing the quantization error of the resulting maps.

While there are previous publications that introduce one or two of the advantages pointed above (e.g. Zhao et al. Zhao et al. [2017] describes a subspace constraint that is not limited to the dictionary items), our work incorporates all of them together in a convenient optimization framework.

Our reconstruction results are based on sampling with variable density spiral trajectories, using 5% and 9% sampling ratios, for retrospective and prospective experiments, respectively. We compare our results to the methods developed by Davies et al. Davies et al. [2014] and Zhao Zhao [2015], and show that FLOR provides quantitative parameter maps with higher accuracy or correspondence to literature compared to those methods.

This paper is organized as follows. Section II describes the MRF problem and provides a review of common reconstruction methods, followed by our low-rank based approach. Section III compares our results to previous MRF algorithms, using retrospective and prospective under-sampled MRF data of a human subject. Sections IV and V discuss experimental results using retrospective and prospective under-sampled MRF data, followed by conclusions.

II Method

II.1 Problem formulation

MRF data consists of multiple frames, acquired in the image’s conjugate Fourier domain (a.k.a k-space), where each frame results from different acquisition parameters. We stack the measurements into a $Q\times L$ matrix $\mathbf{Y}$, where $L$ is the number of frames and $Q$ is the number of k-space samples in each frame. Every column in $\mathbf{Y}$ is an under-sampled Fourier transform of an image frame, $\mathbf{X}_{:,i}$:

[TABLE]

where $F_{u}\{\cdot\}$ denotes an under-sampled 2D Fourier transform and $\mathbf{H}$ denotes a zero mean complex Gaussian noise. The row $\mathbf{X}_{j,:}$ represents the temporal signature of a single pixel (assumed to correspond to a single tissue). The signature depends on the tissue’s relaxation times, T1 and T2, and its proton density (PD), grouped as a row vector:

[TABLE]

Each column, $\mathbf{X}_{:,i}$ represents a response image acquired at a single time point with different acquisition parameters, stacked as a column vector:

[TABLE]

where TR and TE are the repetition time and time to echo and FA represents the flip angle of the RF pulse. Therefore, $\mathbf{X}_{j,:}=f(\mathbf{\Theta}_{1}^{j},\mathbf{\Theta}_{2})$, where $f\{\cdot\}$ represents the Bloch equations. Note that we omit the off resonance parameter (which appeared in $\mathbf{\Theta}_{1}$ in the original MRF paper Ma et al. [2013]), since the sequence used in our retrospective experiments is derived from the FISP sequence, which is insensitive to off resonance effects Jiang et al. [2015].

The goal in MRF is to recover, from the measurements $\mathbf{Y}$, the imaging contrasts $\mathbf{X}$ and the underlying quantitative parameters of each pixel defined in (2), under the assumptions that every pixel in the image contains a single type of tissue and that $\mathbf{\Theta}_{{2}}$ is known.

Recovery is performed by defining a dictionary that consists of simulating the signal generated from $M$ tissues using the Bloch equations (represented as $M$ different combinations of T1 and T2 relaxation times), when the length-$L$ acquisition sequence defined in (3) is used. As a result, we obtain a dictionary $\mathbf{D}$ of dimensions $M\times L$ ($M>L$ as the number of simulated tissues is greater than the sequence length). The PD is not simulated in the dictionary, as it is the gain used to match the Bloch simulation performed on a single spin to the signal obtained from a pixel containing multiple spins. It can be easily determined after the T1 and T2 maps are known. After successful recovery of $\mathbf{X}$, each row in $\mathbf{X}$ is matched to a single row in the dictionary, and T1 and T2 are estimated as those used to generate the matched dictionary row. Each dictionary signature has its own unique T1 and T2 values stored in a look up table (LUT), represented as the matrix $\mathbf{LUT}$ of dimensions $M\times 2$.

II.2 Previous Methods

The approach suggested in the original MRF paper Ma et al. [2013] is described in Algorithm 1, and uses matched filtering to match dictionary items to the acquired data. In the algorithm, $F^{H}\{\cdot\}$ is the 2D inverse NUFFT operator. The parameters $k_{j}$ are the matching dictionary indices, $j$ is a spatial index and $i$ is the temporal index, representing the $i$th frame in the acquisition. The parameter maps are extracted from $\mathbf{LUT}$, which holds the values of T1 and T2 for each $k_{j}$. This approach does not incorporate sparse based reconstruction, which has been proven to be very successful in MRI applications based on under-sampled data Lustig et al. [2008], Weizman et al. [2015, 2016].

Davies et al. Davies et al. [2014] suggested a method incorporating sparsity of the data in the dictionary domain (i.e. each pixel is represented by at most one dictionary item), referred to as the BLoch response recovery via Iterative Projection (BLIP) algorithm. This approach is based on the Projected Landweber Algorithm (PLA) which is an extension of the popular iterative hard thresholding method. BLIP (described here as Algorithm 2) consists of iterating between two main steps: A gradient step that enforces consistency with the measurements, and a projection that matches each row of $\mathbf{X}$ to a single dictionary atom.

BLIP and a few other works that are based on it Wang et al. [2016] do not incorporate the temporal similarity across time points, which is a fundamental characteristic of the MRF sequence. In addition, there is a high degree of similarity across signatures in $\mathbf{D}$. As a result, the imaging contrasts matrix $\mathbf{X}$ is typically a low-rank matrix.

Low-rank based modelling for dynamic MRI in general Feng et al. [2013], Chiew et al. [2015] and MRF in particular McGivney et al. [2014]-Assländer et al. [2017] has been proposed in the past. To demonstrate the low-rank property of $\mathbf{X}$, we used T1, T2 and PD maps of size $128\times 128$ (acquired using DESPOT Deoni et al. [2005]) as an input to a simulation of a FISP sequence Jiang et al. [2015], using $L=500$ TRs. In addition, we used random TR and FA values that have been used in previous publications in the field of MRF Ma et al. [2013], Jiang et al. [2015]. Note that the general assumption of $\mathbf{X}$ being a low-rank matrix holds as long as temporal similarity exist between time-frames in $\mathbf{X}$, and multiple voxels in the image belong to a single tissue, regardless of the specific acquisition parameters. Figure 3 shows the singular values of $\mathbf{X}$. It can be seen that $\mathbf{X}$ is indeed low-rank, as most of the data is represented in the highest 15 singular values.

This low-rank property of $\mathbf{X}$ can be exploited for improved reconstruction using the following optimization problem:

[TABLE]

where $\mathbf{R}_{1}$ is a matrix that matches each pixel ($\mathbf{X}_{j,:}$) with the dictionary signatures. In many previous implementations of low-rank for MRF, a matching of a single dictionary atom to a single pixel is enforced, which means that the rows of $\mathbf{R}_{1}$ are one-sparse vectors. The parameter $r$ is the rank of the matrix, and is usually defined as a fixed pre-chosen parameter. Typically $r$ is not known in advance and determining it upfront arise difficulty and may add error to the reconstruction scheme.

Zhao et al. Zhao [2015] suggested an approximation for problem (4), using an ADMM formulation Boyd et al. [2011] as follows:

[TABLE]

where the low rank constraint is applied via the function $\psi(\mathbf{Z})$, defined as the $p$ norm ($p<1$) of the singular values of $\mathbf{Z}$ to the power of $p$. The matrices $\mathbf{Q}$ and $\mathbf{W}$ are the Lagrange multipliers. The algorithm, coined Model Based Iterative Reconstruction MRF (MBIR-MRF) Zhao [2015] is described in Algorithm 3.

II.3 Proposed Method

The constraint presented in previous approaches Zhao [2015] on $\mathbf{R}_{1}$ to have one sparse rows that contain the corresponding PD values for each row of $\mathbf{X}$, is justified by the assumption that only a single dictionary item should match an acquired signature. However, in practice, we found that superior results (in terms of spatial resolution and correspondence to ground truth) are obtained by relaxing this constraint, and allowing $\mathbf{X}$ to be comprised of multiple dictionary elements at each step of the optimization algorithm, where at the final stage each voxel is matched to a single tissue. This allows for non-simulated signatures to be described by a linear combination of simulated ones. In addition, the relaxation enables formulating the problem as a convex problem, and saves the pattern recognition search time during reconstruction. The matching between $\mathbf{X}$ and the dictionary is done only at the final stage, after $\mathbf{X}$ is fully recovered by using a matched filter (MF), in order to extract the parameter maps. For brevity we write the constraint $\mathbf{X}=\mathbf{RD}$ as $\mathbf{X}\in\mathbb{D}$ where $\mathbb{D}=\{\mathbf{X}:\mathcal{N}(\mathbf{X})\supseteq\mathcal{N}(\mathbf{D})\}$, and we consider the next regularized form:

[TABLE]

for some fixed regularization parameter $\lambda$.

Problem (6) is not convex due to the rank constraint. We therefore relax this constraint by replacing the rank of $\mathbf{X}$ with the nuclear norm $\left\|\mathbf{X}\right\|_{*}$, defined as the sum of the singular values of $\mathbf{X}$ Cai et al. [2010]. This results in the relaxed problem:

[TABLE]

In order to solve (7) we use the incremental subgradient proximal method Sra et al. [2012] as described in Appendix A.

Due to the convex modelling of the problem, we also introduce an improvement that significantly reduces convergence time. The improvement uses the acceleration approach suggested by Nesterov Nesterov [1983] for minimizing a smooth convex function, and its extension for non smooth composite functions of Beck and Teboulle Beck and Teboulle [2009], Palomar and Eldar [2010]. Our final algorithm is detailed in Algorithm 4 and referred to as magnetic resonance Fingerprint with LOw Rank (FLOR), where the parameter $\lambda$ is chosen experimentally. Note that by setting $\lambda=0$, enforcing $\mathbf{R}$ to have one-sparse rows and eliminating the acceleration step, FLOR reduces to BLIP Davies et al. [2014].

Figure 4 shows the reconstruction error of FLOR as the number of iterations varies with and without the acceleration step. Note that the CPU time of both algorithms is similar.

II.4 Possible extension

Conventional MRF algorithms use MF for the magnetic parameter extraction. MF introduces a quantization error since map values are continuous, as opposed to discrete dictionary values. A possible extension of FLOR is to add values to the dictionary by linear interpolation, in regions where a few candidates from the dictionary match a single signature from the data. We then select the dictionary signatures that exhibit a high correlation value (the ones above a certain threshold) and average their matching T1 and T2 values. This improvement expands the possible solutions to include ones that do not exist in the dictionary, and therefore exhibits improved accuracy compared to the conventional matching. The major benefit from this extension is reduced quantization errors that arise from conventional MF used in MRF. This extension, coined FLOR II, is examined in the first part of our experimental results in the next section.

III Experimental Results

This section describes two MRI experiments that were carried out using brain scans of a healthy subject. The first experiment is based on well known quantitative maps that were used, in a purely simulation environment to generate an MRF experiment with retrospectively sampled data. While this experiment is a simulation based on real quantitative maps, it allows accurate comparison of the results of the different algorithms using a well defined reference.

In the second experiment, we used prospective sampled real MRF data that was used to generate the results in Ma et al. Ma et al. [2013]. While this experiment lacks a gold-standard for accurate error evaluation, it allows comparison between different algorithms in a realistic multi-coil acquisition. To compare between different algorithms, for prospective sampling, where no ground truth is available, we examined the performance of the various algorithms as a function of the total number of excitations, where correspondence to values provided in literature for various brain tissues is used for validation. In both experiments, variable density spiral trajectories were used for sampling.

For quantitative error analysis, we calculated the normalized MSE (NMSE) between each quantitative map estimation and the reference map, defined as:

[TABLE]

where $\theta_{i}$, $\hat{\theta}_{i}$ represent a reference map (such as T1,T2 or PD) and its corresponding reconstructed map (respectively), $N$ is the number of pixels in the map and $j$ is a spatial index.

In the first experiment, forward and inverse non-uniform Fourier transforms were applied using SPURS, which is a fast approach published recently Kiperwas et al. [2016]. For the second experiment, we used the NUFFT package Fessler and Sutton [2003], to adhere with the reconstruction results of the original MRF paper Ma et al. [2013].

III.1 Experiment 1: Retrospective undersampling of real data

The data for this experiment was acquired with a GE Signa 3T HDXT scanner. The procedures involving human subjects described in this experiment were approved by the Institutional Review Board of Tel-Aviv Sourasky Medical Center, Israel. We generated our reference data by acquisition of Fast Imaging Employing Steady-state Acquisition (FIESTA) and Spoiled Gradient Recalled Acquisition in Steady State (SPGR) images, at 4 different flip angles ($3^{\circ}$ ,$5^{\circ}$,$12^{\circ}$ and $20^{\circ}$), implementing the fast and well known DESPOT1 and DESPOT2 Deoni et al. [2005] algorithms, after improvements as described in Liberman et al. Liberman et al. [2014], to generate T1,T2 and PD quantitative maps, each of size $128\times 128$ pixels. While it is well known that the gold standard method for T1 measurement is the inversion recovery spin echo with varying TIs and for T2 measurement is the spin echo sequences with varying TEs, in this experiment DESPOT was used as a reference thanks to its availability and its relatively fast acquisition time. The FISP pulse sequence has been applied for simulating acquisition of the reference. It was simulated with constant TE of 2ms, random TR values in the range of 11.5-14.5 ms, and a sinusoidal variation of FA (RF pulses) in the range of 0-70 degrees Jiang et al. [2015].

To simulate noisy undersampled MRF samples, we added complex Gaussian zero-mean noise to the k-space data to obtain an SNR of 67dB in the undersampled measurement domain. Data was then under-sampled to acquire only 876 k-space samples in each TR with spiral trajectories. In particular, we used 24 variable density spirals with inner region size of 20 and FOV of 24. In every time frame, each spiral is shifted by 15 degrees. Figure 5 demonstrates the first spiral trajectory. We define the under-sampling ratio by the number of the acquired samples in the k-space domain divided by the number of pixels in the generated image. This leads to an undersampling ratio of ~5% in this experiment. For comparison, the under-sampling ratio of the original MRF paper Ma et al. [2013] is ~9%, since for each single spiral 1450 data points were acquired.

We generated the dictionary using Bloch equations, simulating T1 values of [100:20:2000,

2300:300:5000] ms and T2 values of [20:5:100,110:10:200,300:200:1900] ms. This range covers the relaxation time values that can be found in a healthy brain scan Vymazal et al. [1999]. The tuning parameters were experimentally set as $\mu=1$ and $\lambda=5$, after $\lambda$ was tested in the range between 0 and 30. Data was fed as an input to BLIP, MBIR-MRF and the improved FLOR algorithm (described as Algorithms 2,3 above and Algorithm 4 with the additional extension of interpolating the parameter maps). In addition, we performed reconstruction using 100% of the data (without the addition of noise) via conventional MRF (Algorithm 1), for comparison purposes and to evaluate the error caused by the effect of discretized dictionary. All the iterative algorithms were run until the difference between consecutive iterations was below the same threshold.

The MATLAB code for reproducing the experiment provided in this section can be be found at: http://webee.techni- on.ac.il/Sites/People/YoninaEldar/software_det18.php. In this code, spiral sampling trajectories design was based on Lee et al. Lee et al. [2003].

Figure 6 shows the resulting maps for the recovery of T1, T2 and PD obtained with the various algorithms against the reference (left). The corresponding error maps of each method versus the reference are shown in Fig. 7. To allow detailed view of the reconstruction results for the reader, Fig. 6 shows a zoomed region for each map.

It can be seen that both FLOR and MBIR-MRF outperform BLIP reconstruction results, when using 5% of sampled data by utilizing the low rank property. In addition, FLOR provides a lower error compared to MBIR-MRF. The details in the FLOR maps are comparable to those obtained by the original MRF algorithm using 100% of the noise-free data. Due to the very low sampling ratio in our experiments (measured as the number of samples divided by the number of pixels in the image), conventional MRF using 5% of the data did not provide valuable reconstruction results and is therefore omitted in this analysis.

We next implemented the MF improvement described in Section II.D. The results are shown in Fig. 8, with corresponding error maps in Fig. 9. These figures compare the recovery maps of FLOR without (FLOR I) and with (FLOR II) the proposed improvement. It can be seen that FLOR II improves the results of FLOR I and produces a smoother solution which better fits the reference maps.

III.2 Experiment 2: In vivo prospective sampling experiment

The experiment in this section was carried out using the data of the original MRF paper Ma et al. [2013], which consisted of 48 spiral trajectories shifted by 7.5 degrees, where 1450 samples were acquired in each trajectory, leading to an underampling ratio of 9%. The data was acquired on a 1.5-T whole body scanner (Espree, SIEMENS Healthcare) using a 32-channel head receiver coil.

Due to the lack of gold standard maps for this data, we are unable to provide quantitative error results (e.g. NMSE). Therefore, in this experiment we compare between the various algorithms by examining reconstruction results using 400 TRs (representing 40% of scanning time), to quantitative values of brain tissues from the literature. Since the results obtained in the original MRF experiment (using 1000 TRs) mostly correspond to quantitative values from the literature, the maps generated using 1000 TRs using the original MRF algorithm are provided in Fig. 10, for reference.

The results of T1, T2 and PD maps for BLIP, MBIR-MRF and FLOR appear in Fig. 11. Since IR-bSSFP sequence has been used in this experiment, off-resonance frequency has also been computed and shown. We used 109 different values in the range between -250 and 240 Hz. It can be seen that for T1, all iterative algorithms provide similar results, and T1 values of grey matter (GM), white matter (WM) and cerebrospinal fluid (CSF) regions correspond to similar values that appear in the literature (see Table 1 in Ma et al.Ma et al. [2013]) and in Fig. 10. While T2 results exhibit visible differences between the various methods, WM and GM values for all methods correspond to values that appear in the literature. However, both BLIP and MBIR MRF exhibit T2 values for CSF that are lower than those reported in the literature. This can be seen in Fig. 12, where the color scale for T2 is adjusted to 500-2000ms (T2 values for CSF are around 2000ms). Shortened T2 values in CSF were also reported in the original MRF experiment with 1000 TRs (and were justified as out-of-plane flow in this 2D experiment). In our case, using the same acquired data, it can be seen in Fig. 12 that FLOR provides CSF values that better correspond to literature values, when compared to the other methods.

IV Discussion

IV.1 Relation to previous works

Although works that exploit the low rank structure of MRF sequences have been published in the past by others McGivney et al. [2014], Cline et al. [2016], Zhao [2015], Zhao et al. [2016b], Liao et al. [2016], Doneva et al. [2017], Zhao et al. [2017], Assländer et al. [2017] and also by us Mazor et al. [2016], our solution is unique mainly in the combination of convex modelling and the ability to enable a solution with quantitative values that do not exist in the dictionary. Our solution is based on soft-thresholding the singular values Cai et al. [2010], which is mathematically justified in Appendix A.

Moreover, we compare our algorithm to both CS-based and low-rank based methods for MRF and demonstrate superior results. While BLIP treats the original MRF problem as an $\ell_{0}$ optimization problem, FLOR first solves the relaxed problem of (6) and only then uses MF to extract the magnetic parameters. It leads to some beneficial properties such as convergence guarantees, and the ability to use the acceleration step as described in Algorithm 4, which is also novel in MRF reconstruction methods.

IV.2 Computational complexity

FLOR is divided into two main components: The first recovers the imaging contrasts, and the second extracts the parameter maps from the recovered contrasts. The computational burden of FLOR lies in the low-rank projection step, or specifically, in the SVD calculation. This step does not exist in BLIP nor the original MRF reconstruction. However, there are efficient fast techniques to calculate the SVD Drmac and Veselic [2005], required by FLOR. Moreover, unlike BLIP, and other low rank based algorithms such as MBIR-MRF, FLOR does not require the pattern recognition calculation at every iteration. Another time consuming step that exists in all algorithms is the non uniform Fourier transform. By using the acceleration step, FLOR reduces significantly the number of iterations required for convergence and therefore saves computational cost.

In addition, while previous implementations of CS-based reconstruction algorithms mainly use the inverse NUFFT (iNUFFT) algorithm, in our retrospective experiemtns we use SPURS Kiperwas et al. [2016]. Based on our observations, SPURS provides improved image reconstruction with the same computational complexity compared to iNUFFT.

V Conclusions

We presented FLOR, a method for high quality reconstruction of quantitative MRI data using MRF, by utilizing the low-rank property of MRF data. Due to the fact that we exploit low-rank on top of the well known sparsity of MRF in the dictionary matching domain, we are able to obtain high quality reconstruction from highly under-sampled data. Our method is based on a convex minimization problem, leading to a solution in the dictionary subspace that overcomes its quantization error.

We provide results that are comparable to fully sampled MRF, using only 5% of the data in a simulation environment. In addition, comparison against CS-based and low-rank based methods for MRF shows the added value of our approach in generating quantitative maps with less artifacts. Our results also consist of real-data, in-vivo experiments that exhibit FLOR superiority also for realistic multi-coil data acquisition. Future work will examine more sophisticated patch wise recoveries.

Acknowledgements

The authors wish to thank the Tel Aviv center for brain functions at Tel Aviv Sourasky Medical Center for providing the data required for experiment 1. We also wish to thank Dan Ma and Prof. Mark Griswold for providing the real data used in their experiments. This work was supported by the Ministry of Science, by the ISF I-CORE joint research center of the Technion and the Weizmann Institute, and by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 646804-ERC-COG-BNYQ, Assaf Tal acknowledges the support of the Monroy-Marks Career Development Fund, the Carolito Stiftung Fund, the Leona M. and Harry B. Helmsley Charitable Trust and the historic generosity of the Harold Perlman Family. The authors have no relevant conflicts of interest to disclose.

Appendix A

The basic implementation of FLOR, as described in Algorithm 4 in the paper, aims to solve the following optimization problem:

[TABLE]

where $F_{u}$ is the partial Fourier transform operator, $\mathbf{X}$ has dimensions $N^{2}\times L$ and $\mathbb{D}=\{\mathbf{X}:\mathcal{N}(\mathbf{X})\supseteq\mathcal{N}(\mathbf{D})\}$.

FLOR solves (9) using the incremental proximal method Sra et al. [2012], which treats problems of the form:

[TABLE]

where $F_{i}(\mathbf{X})=f_{i}(\mathbf{X})+h_{i}(\mathbf{X})$. The function $f_{i}(\mathbf{X})$ is convex and non-differentiable, $h_{i}(\mathbf{X})$ is a convex function and $\mathbb{D}$ is a non-empty, closed, and convex subspace. The general step in solving (10) is given by [27, (4.12)-(4.13)]:

[TABLE]

where $g_{i_{k}}\in\partial h_{i_{k}}(\mathbf{X}^{k})$, $\mu_{k}$ is a positive step size, and $\text{P}_{\mathbb{D}}$ is the projection operator onto ${\mathbb{D}}$ defined as

[TABLE]

The optimization problem, defined in the update step of $\mathbf{X}^{k+1}$, is referred to as the proximal gradient calculation of the non-differentiable $f_{i_{k}}$, under the constraint $\mathbf{X}\in\mathbb{D}$.

Our problem in (9) corresponds to $m=1$ in (10) and

[TABLE]

Therefore,

[TABLE]

and,

[TABLE]

The solution of (11b) for $f(\mathbf{X})=\lambda\|\mathbf{X}\|_{*}$ without the constraint $\mathbf{X}\in\mathbb{D}$, is the singular value soft-thresholding operator (SVT) Cai et al. [2010] defined as:

[TABLE]

Here $\mathbf{\Sigma}_{r}$ is a diagonal matrix with the non-zero singular values of $\mathbf{Z}^{k}$ on its diagonal, $\mathbf{U}_{r}$ and $\mathbf{V}_{r}$ are the $r$ left and right singular vectors of the SVD of $\mathbf{Z}^{k}$, associated with the $r$ non-zero singular values, and $[x]_{+}=\text{max}(0,x)$. In our case, since $\mathbf{Z}^{k}\in\mathbb{D}$ (as follows from (11a)) and the SVT calculation keeps the operand in the same subspace, the constraint $\mathbf{X}\in\mathbb{D}$ can be omitted. Therefore, (11b) reduces to

[TABLE]

Combining (14), (17) and (15), the incremental subgradient-proximal method for solving (9) consists of two updates in each iteration:

[TABLE]

This constitutes the core of Algorithms 4. In our framework, the step sizes are set to constant, $\mu_{k}=\mu$, and $\lambda$ is chosen experimentally.