HYDRA: Hybrid Deep Magnetic Resonance Fingerprinting

TL;DR

HYDRA introduces a hybrid deep learning framework for magnetic resonance fingerprinting that enhances speed and accuracy by replacing traditional dictionary matching with neural networks, reducing errors and storage needs.

Contribution

The paper presents a novel two-stage HYDRA approach combining model-based signal restoration and deep learning for continuous parameter estimation in MRF, overcoming discretization and computational issues.

Findings

Significantly faster inference compared to traditional methods.

Produces continuous-valued tissue parameters, reducing discretization errors.

Requires less memory by eliminating large dictionaries.

Abstract

Purpose: Magnetic resonance fingerprinting (MRF) methods typically rely on dictio-nary matching to map the temporal MRF signals to quantitative tissue parameters. Such approaches suffer from inherent discretization errors, as well as high computational complexity as the dictionary size grows. To alleviate these issues, we propose a HYbrid Deep magnetic ResonAnce fingerprinting approach, referred to as HYDRA. Methods: HYDRA involves two stages: a model-based signature restoration phase and a learning-based parameter restoration phase. Signal restoration is implemented using low-rank based de-aliasing techniques while parameter restoration is performed using a deep nonlocal residual convolutional neural network. The designed network is trained on synthesized MRF data simulated with the Bloch equations and fast imaging with steady state precession (FISP) sequences. In test mode, it takes a temporal MRF signal as input and produces the corresponding tissue parameters. Results: We validated our approach on both synthetic data and anatomical data generated from a healthy subject. The results demonstrate that, in contrast to conventional dictionary-matching based MRF techniques, our approach significantly improves inference speed by eliminating the time-consuming dictionary matching operation, and alleviates discretization errors by outputting continuous-valued parameters. We further avoid the need to store a large dictionary, thus reducing memory requirements. Conclusions: Our approach demonstrates advantages in terms of inference speed, accuracy and storage requirements over competing MRF methods

Tables

Table 1. Table 1: Brief description of experiment types and settings.

Experiment	Settings
Training	Input: $𝐃$, size $K\times L=80100\times 200$. Groundtruth: $\mathrm{𝐋𝐔𝐓}$, size $80100\times 2$. k-space subsampling factor $\beta$: not available.
Testing on synthetic data	Input: $𝐗$, size $N\times L=80000\times 200$. Groundtruth: ${𝚯}^{T12}$, size $80000\times 2$. k-space subsampling factor $\beta$: not available.
Testing on anatomical data	Input: $𝐘$, size $Q\times L=16384\beta \times 200$ or $16384\beta \times 1000$. Reference: ${𝚯}^{T12}$, size $N\times 2=16384\times 2$. k-space subsampling factor $\beta$: 70%, 15% using Gaussian patterns.
Testing on anatomical data	Input: $𝐘$, size $Q\times L=16384\beta \times 1000$. Reference: ${𝚯}^{T12}$, size $N\times 2=16384\times 2$. k-space subsampling factor $\beta$: 9% using spiral trajectories.

Table 2. Table 2: Testing on synthetic dataset. Comparing parameter restoration performance, in terms of PSNR, SNR, RMSE and correlation coefficient.

	Dict. Match.	CNN 14	FNN 13	Proposed
	T1 / T2	T1 / T2	T1 / T2	T1 / T2
PSNR (dB)	59.15 / 52.31	62.96 / 49.64	58.97 / 54.96	79.30 / 72.99
SNR (dB)	55.23 / 47.15	59.05 / 44.49	55.06 / 49.81	75.38 / 67.83
RMSE (ms)	5.515 / 4.847	3.554 / 6.591	5.63 / 3.57	0.542 / 0.448
CorrCoef	1.00 / 1.00	1.00 / 1.00	1.00 / 1.00	1.00 / 1.00
time cost (s)	464.10	2.87	1.58	8.2

Table 3. Table 3: Testing on synthetic dataset involving detailed T1 / T2 examples that are not on the training grid and their intervals are much smaller than the training grid intervals. D.M. denotes dictionary matching. T1 and T2 errors are defined as the difference between estimated values and groundtruth values.

	T1 Estimation				T1 Errors
Truth	D.M.	14	13	Ours	D.M.	14	13	Ours
1005.0	1001.0	1002.9	1009.3	1004.8	-4.0	-2.1	4.3	-0.2
1005.5	1001.0	1003.3	1010.0	1005.3	-4.5	-2.3	4.5	-0.2
1006.0	1011.0	1003.6	1010.6	1005.8	5.0	-2.4	4.6	-0.2
1006.5	1011.0	1004.1	1011.2	1006.3	4.5	-2.5	4.7	-0.2
1007.0	1011.0	1004.5	1011.8	1006.8	4.0	-2.5	4.8	-0.3
RMSE	-	-	-	-	4.4	2.3	4.6	0.2
	T2 Estimation				T2 Errors
Truth	D.M.	14	13	Ours	D.M.	14	13	Ours
505.0	501.0	513.6	504.3	505.2	-4.0	8.6	-0.7	0.2
505.5	511.0	514.1	504.8	505.7	5.5	8.6	-0.7	0.2
506.0	501.0	514.7	505.3	506.2	-5.0	8.6	-0.7	0.2
506.5	511.0	515.2	505.9	506.8	4.5	8.7	-0.6	0.3
507.0	511.0	515.7	506.4	507.3	4.0	8.7	-0.6	0.3
RMSE	-	-	-	-	4.6	8.6	0.7	0.2

Table 4. Table 4: Testing on anatomical dataset with full k-space sampling. Comparing parameter restoration performance, in terms of PSNR, SNR, RMSE and correlation coefficient.

	Dict. Match.	CNN 14	FNN 13	Proposed basic	Proposed nonlocal
	T1 / T2	T1 / T2	T1 / T2	T1 / T2	T1 / T2
PSNR (dB)	56.64 / 52.04	54.06 / 49.88	54.53 / 54.36	56.59 / 60.01	56.47 / 62.56
SNR (dB)	42.20 / 27.81	39.63 / 25.66	40.09 / 30.07	42.15 / 35.76	42.03 / 38.32
RMSE (ms)	6.623 / 6.252	8.912 / 8.015	8.45 / 4.78	6.661 / 2.498	6.76 / 1.86
CorrCoef	1.00 / 1.00	1.00 / 1.00	1.00 / 1.00	1.00 / 1.00	1.00 / 1.00
time cost (s)	84.56	0.69	0.41	1.6	2.1

Table 5. Table 5: Testing on anatomical dataset with k-space subsampling ratio 70% and 15% using Gaussian patterns and 200 time frames.

	k-space subsampling factor $\beta$ = 70%
	Ma et al. 1	BLIP 3	FLOR 6	CNN 14	FNN 13	Proposed
						basic	nonlocal
PSNR (dB)	23.69 / 38.17	45.67 / 47.84	50.11 / 50.85	49.71 / 45.48	50.15 / 51.08	50.79 / 51.59	49.87 / 57.57
SNR (dB)	8.73 / 13.84	31.28 / 23.49	35.67 / 26.48	35.26 / 21.19	35.70 / 26.67	36.34 / 27.19	35.42 / 33.30
RMSE (ms)	294.32 / 30.87	23.42 / 10.14	14.01 / 7.17	14.71 / 13.31	13.99 / 6.98	12.99 / 6.57	14.44 / 3.31
time cost (s)	72.88	75.70	85.35	23.72	23.53	24.85	26.3
	k-space subsampling factor $\beta$ = 15%
	Ma et al. 1	BLIP 3	FLOR 6	CNN 14	FNN 13	Proposed
						basic	nonlocal
PSNR (dB)	27.94 / 32.84	35.45 / 39.25	44.95 / 46.11	43.74 / 35.98	45.03 / 45.90	45.23 / 44.44	45.39 / 51.32
SNR (dB)	13.50 / 8.61	20.99 / 14.58	30.51 / 21.89	29.23 / 12.26	30.58 / 21.32	30.76 / 19.78	30.91 / 26.99
RMSE (ms)	180.3 / 57.03	76.01 / 27.25	25.46 / 12.37	29.27 / 39.73	25.21 / 12.68	24.65 / 15.00	24.20 / 6.79
time cost (s)	106	112.8	121.7	24.54	24.36	25.67	27.31

Table 6. Table 6: Testing on anatomical dataset with k-space subsampling ratio 15% using Gaussian patterns and 1000 time frames.

	Ma et al. 1	BLIP 3	FLOR 6	CNN 14	FNN 13	Proposed
PSNR (dB)	27.53 / 33.28	35.50 / 39.10	50.90 / 50.04	41.96 / 39.21	52.62 / 49.86	52.32 / 52.79
SNR (dB)	13.09 / 9.05	21.06 / 14.87	36.44 / 25.65	27.44 / 15.05	38.17 / 25.43	37.86 / 28.35
RMSE (ms)	189.09 / 54.21	75.53 / 27.74	12.83 / 7.87	35.91 / 27.37	10.52 / 8.04	10.89 / 5.74

Table 7. Table 7: Testing on anatomical dataset with k-space subsampling ratio 9% using spiral trajectories and 1000 time frames.

	Ma et al. 1	BLIP 3	FLOR 6	CNN 14	FNN 13	Proposed
PSNR (dB)	26.66 / 30.44	29.35 / 39.47	39.32 / 44.60	35.68 / 27.74	40.26 / 44.70	41.45 / 45.41
SNR (dB)	12.22 / 6.21	15.03 / 15.22	24.88 / 20.38	21.45 / 4.48	25.84 / 20.30	27.02 / 21.04
RMSE (ms)	209.01 / 75.18	153.37 / 26.57	48.67 / 14.72	73.96 / 102.57	43.65 / 14.55	38.08 / 13.41

Equations

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

$${\mathbf{\Theta }}_{:,i}^{TRE}=[T{R}^i,T{E}^i,F{A}^i{]}^T,i\in [1,L]$$

$${\mathbf{\Theta }}_{j,:}^{T12}=[T{1}^j,T{2}^j],j\in [1,N]$$

$${\mathbf{X}}_{j,:}=f({\mathbf{\Theta }}_{j,:}^{T12},{\mathbf{\Theta }}^{TRE})$$

$${\mathbf{\Theta }}_{j,:}^{T12}=g(h(\mathbf{Y}{)}_{j,:}|{\mathbf{\Theta }}^{TRE}),j\in [1,N],$$

$${\hat{\mathbf{X}}}_{:,i}=F_u^H\{{\mathbf{Y}}_{:,i}\},\forall i$$

$$\widehat{k_j}=\arg \underset{k}{\max }\frac{\text{Re}⟨{\mathbf{D}}_{k,:},{\hat{\mathbf{X}}}_{j,:}⟩}{\Vert {\mathbf{D}}_{k,:}{\Vert }_2^2},{\hat{T}}_1^j,{\hat{T}}_2^j=\mathbf{LUT}(\widehat{k_j})$$

$${\hat{\mathbf{Z}}}_{:,i}^{t+1}={\hat{\mathbf{X}}}_{:,i}^t-\mu F_u^H\{F_u\{{\hat{\mathbf{X}}}_{:,i}^t\}-{\mathbf{Y}}_{:,i}\}$$

$$[\mathbf{U},\mathbf{S},\mathbf{V}]=\mathrm{svd}({\hat{\mathbf{Z}}}^{t+1}\mathbf{P})$$

$$\begin{array}{r}{\sigma }_j^{\prime }=\max \{{\sigma }_j-\lambda \mu ,0\},{\hat{\mathbf{X}}}^{t+1}=\mathbf{U}{\mathbf{S}}^{\prime }{\mathbf{V}}^H\end{array}$$

$$\widehat{k_j}=\arg \underset{k}{\max }\frac{\text{Re}⟨{\mathbf{D}}_{k,:},{\hat{\mathbf{X}}}_{j,:}⟩}{\Vert {\mathbf{D}}_{k,:}{\Vert }_2^2},{\hat{T}}_1^j,{\hat{T}}_2^j=\mathbf{LUT}(\widehat{k_j})$$

$$\begin{array}[]{rl}h(\mathbf{Y})=\arg\underset{\mathbf{X}}{\min}&\frac{1}{2}\sum_{i}\|\mathbf{Y}_{:,i}-F_{u}\{\mathbf{X}_{:,i}\}\|_{2}^{2}\\ \text{s.t.}&\textrm{rank}(\mathbf{X})<r\end{array}$$

$$h(\mathbf{Y})=\arg \underset{\mathbf{X}}{\min }\frac{1}{2}{\sum }_i\Vert {\mathbf{Y}}_{:,i}-F_u\{{\mathbf{X}}_{:,i}\}{\Vert }_2^2+\lambda \Vert \mathbf{X}{\Vert }_{\ast }$$

$$[\mathbf{U},\mathbf{S},\mathbf{V}]=\mathrm{svd}({\hat{\mathbf{Z}}}^{t+1})$$

$${\hat{T}}_1^j,{\hat{T}}_2^j=g({\hat{\mathbf{X}}}_{j,:})$$

$${\mathbf{\Theta }}_{j,:}^{T12}=g({\mathbf{X}}_{j,:}|{\mathbf{\Theta }}^{TRE}),j\in [1,N]$$

$${\mathbf{y}}_i=\frac{1}{C(\mathbf{x})}\sum _{\forall j}f({\mathbf{x}}_i,{\mathbf{x}}_j)g({\mathbf{x}}_j).$$

$$=\sqrt{\frac{\Vert \mathbf{X}-\hat{\mathbf{X}}{\Vert }_F^2}{N}},$$

$$=20{\log }_{10}\frac{\Vert \mathbf{X}{\Vert }_F^2}{\text{RMSE}},$$

$$=20{\log }_{10}\frac{\text{PeakVal}}{\text{RMSE}},$$

Full text

Methods: HYDRA involves two stages: a model-based signature restoration phase and a learning-based parameter restoration phase. Signal restoration is implemented using low-rank based de-aliasing techniques while parameter restoration is performed using a deep nonlocal residual convolutional neural network. The designed network is trained on synthesized MRF data simulated with the Bloch equations and fast imaging with steady state precession (FISP) sequences. In test mode, it takes a temporal MRF signal as input and produces the corresponding tissue parameters.

Results: We validated our approach on both synthetic data and anatomical data generated from a healthy subject. The results demonstrate that, in contrast to conventional dictionary-matching based MRF techniques, our approach significantly improves inference speed by eliminating the time-consuming dictionary matching operation, and alleviates discretization errors by outputting continuous-valued parameters. We further avoid the need to store a large dictionary, thus reducing memory requirements.

Conclusions: Our approach demonstrates advantages in terms of inference speed, accuracy and storage requirements over competing MRF methods.**

Key words: Magnetic Resonance Fingerprinting, Quantitative Magnetic Resonance Imaging, Deep Learning, Nonlocal Residual Convolutional Neural Network, Self-attention

This work was supported by the Royal Society International Exchange Scheme IE160348, by the European Union’s Horizon 2020 grant ERC-BNYQ, by the Israel Science Foundation grant no. 335/14, by ICore: the Israeli Excellence Center ’Circle of Light’, by the Ministry of Science and Technology, Israel, by UCL Overseas Research Scholarship (UCL-ORS) and by China Scholarship Council (CSC) and by EPSRC grant EP/K033166/1.

1 Department of Electronic and Electrical Engineering, Imperial College London, UK

2 Faculty of Mathematics and Computer Science, Weizmann Institute of Science, Israel

3 Department of Electrical Engineering, Technion – Israel Institute of Technology, Israel

4 Department of Electronic and Electrical Engineering, University College London, UK

p.song@imperial.ac.uk, yonina.eldar@weizmann.ac.il, galmazor@technion.ac.il, m.rodrigues@ucl.ac.uk *

I. Introduction

Magnetic Resonance Fingerprinting (MRF) 1, 2, 3, 4, 5, 6, 7, 8 has emerged as a promising Quantitative Magnetic Resonance Imaging (QMRI) approach, with the capability of providing multiple tissue’s intrinsic spin parameters simultaneously, such as the spin-lattice magnetic relaxation time (T1) and the spin-spin magnetic relaxation time (T2). Based on the fact that the response from each tissue with respect to a given pseudo-random pulse sequence is unique, MRF exploits pseudo-randomized acquisition parameters to create unique temporal signal signatures, analogous to a "fingerprint", for different tissues. A dictionary matching operation is then performed to map an inquiry temporal signature to the best matching entry in a precomputed dictionary, leading to multiple tissue parameters directly.

The temporal signatures are generated by varying the acquisition parameters of a pseudo-random excitation pulse sequence, such as repetition time (TR), time of echo (TE), and radio frequency flip angle (FA) over time. The dictionary is composed of a large number of entries that are usually simulated by the Bloch equations given pseudo-random pulse sequences. Each entry represents a unique temporal signature associated with a specific tissue and its quantitative parameters, such as the T1 and T2 relaxation times. Thus, once the best matching (i.e. most correlated) entry is found, it directly leads to multiple quantitative parameters simultaneously via a lookup-table operation.

MRI physics and physiological constraints make the MR scanning procedure time-consuming. To shorten acquisition time, subsampling is commonly performed in k-space (a.k.a conjugate Fourier transform domain) in order to reduce the number of samples and accelerate imaging speed. However, such k-space subsampling results in temporal signatures that are corrupted by aliasing, blurring and noise. This hampers the accuracy associated with estimation of the tissue parameters using a dictionary matching procedure. In order to alleviate the impact of such distortion and corruption, de-aliasing operations are often exploited to restore cleaner signatures before performing signature-to-parameter mapping. Therefore, MRF reconstruction usually involves two operations: signature restoration and parameter restoration.

Inspired by the successful application of sparsity-driven image processing approaches in MRI reconstruction 9, 10, 11, 12, several works 3, 4, 5, 6 suggest to incorporate prior knowledge such as sparsity and low-rank to attenuate distortion and corruption, improving the signature restoration performance, during the initial MRF reconstruction stage. This is then followed by a dictionary matching operation, performing mapping from purified temporal signatures to tissue’s quantitative parameters. However, such dictionary matching based signature-to-parameter mapping exhibits several drawbacks 13, 14. Since the simulated dictionary and lookup-table contain a finite number of elements, they can only cover a limited number of discrete values for each type of tissue parameter. We refer to the difference between a continuous-valued tissue parameter and its closest available discrete value on a lattice as the discretization error. For example, a pair of dictionary and lookup-table that contain 101 elements will lead to a discretization error of maximum 25 ms if they cover the range of 0 ms - 5000 ms with a fixed interval of 50 ms for a specific tissue parameter, e.g. T1. To reduce the discretization error, a huge dictionary that is composed of a large number of entries is needed to represent tissues with fine granularity over the entire value range of target tissue parameters. However, storing a large dictionary becomes prohibitively memory-consuming, as the dictionary size and density often increase exponentially with the number of tissue parameters. Specifically, the number of entries in a dictionary will be $P^{s}$ for $s$ parameters each containing $P$ values, since every combination of these $s$ parameters determines a specific tissue which is characterized by a specific signature. For example, given T1, T2 relaxation times, i.e. $s=2$, if each of them contains $1000$ values, the dictionary will have $1000^{2}$ entries. In addition, finding the best matching entry becomes computationally intense for a large dictionary, considerably limiting the inference speed.

In this paper, we propose an alternative approach to dictionary matching based on deep neural networks (a.k.a. deep learning) 15, 16, which we refer to as HYDRA: HYbrid Deep magnetic ResonAnce fingerprinting. The motivation derives from the fact that a well designed and tuned deep neural network is capable of approximating complex functions, leading to state-of-the-art results in a number of tasks such as image classification, super-resolution, speech recognition, and more 17, 18, 19, 20, 21, 22, 23. Recent work 13, 14 proposed to exploit neural networks to replace the dictionary and lookup-table used in conventional MRF reconstruction approaches. These proposed neural networks suffer from two limitations: First, these approaches are based on neural network models containing only 3-layers, thus suffer from limited capacity of capturing complex mapping functions. Second, these methods focused exclusively on parameter restoration stage (the second stage in MRF reconstruction), but not on signature restoration (the first stage in MRF reconstruction). Therefore, these techniques rely on fully-sampled data instead of typically available sub-sampled k-space data.

Different from Cohen et al.’s fully-connected feed-forward neural network 13, and Hoppe et al.’s vanilla convolutional neural network (CNN) 14, the proposed HYDRA involves both a signature restoration and a parameter restoration phase. Signature restoration is implemented using a low-rank based de-aliasing method adapted from Mazor et al. 6 while parameter restoration is implemented using a deep nonlocal residual convolutional neural network developed for this purpose. Our key contributions with respect to prior work are:

• •

HYDRA is, to the best of our knowledge, the first deep network approach to combine model-based de-aliasing and learning-based parameter mapping. HYDRA eliminates the requirement for the memory and time-consuming dictionary matching operation, thus significantly improving inference speed without compromising on reconstruction performance.

• •

A 1D nonlocal residual convolutional neural network is designed to capture the mappings from temporal MRF signals to tissue parameters. Owing to residual learning and a self-attention mechanism, our network is deeper and more sophisticated than competing network models. This allows to capture complex parameter mappings more effectively, and output continuous parameters to alleviate discretization issues.

• •

The designed network is trained on synthesized MRF data simulated with the Bloch equations, but is still applicable to anatomical data. This contributes to eliminating the requirement for a large amount of real MRF data.

• •

A low-rank based de-aliasing technique is developed in order to take advantage of temporal similarity for signature restoration.

• •

The low-rank based signature restoration is organically combined with the learning-based parameter restoration to achieve fast and accurate MRF reconstruction. Such strategy enables HYDRA to handle both fully-sampled k-space data and more importantly sub-sampled k-space data.

• •

A series of numerical experiments are conducted to evaluate the proposed approach on both synthetic and anatomical data. The results demonstrate improved inference speed, accuracy and discretization errors over competing methods 1, 2, 3, 4, 5, 6, 14, 13.

The rest of the paper is organized as follows. In Section II., we formulate the MRF reconstruction problem, introduce related methods, and present our approach, involving the use of a low-rank based signature restoration procedure together with a deep network for parameter restoration. Section III. is devoted to experimental results, followed by a discussion in Section IV. and a conclusion in Section V..

II. Materials and Methods

II.A. The MRF Problem Formulation

MRF data is composed of multiple frames sampled in k-space over time. A series of such frames are vectorized and then stacked together along the temporal dimension to construct a measurement matrix $\mathbf{Y}\in\mathcal{C}^{Q\times L}$, where $Q$ is the number of k-space samples in each frame, and $L$ is the number of frames. Due to k-space subsampling, every column vector $\mathbf{Y}_{:,i}$ represents a subsampled Fourier transform of a vectorized image frame $\mathbf{X}_{:,i}$:

[TABLE]

where $F_{u}\{\cdot\}$ denotes a subsampled 2D Fourier transform.

Each column $\mathbf{X}_{:,i}$ represents a MR contrast acquired with RF sequence parameters:

[TABLE]

where $TR^{i}$ and $TE^{i}$ denote the repetition time and echo time, respectively, and $FA^{i}$ denotes the flip angle of the RF pulse during sampling the $i$-th contrast. Every row $\mathbf{X}_{j,:}$ represents a temporal signature, i.e. temporal signal evolution of a specific tissue at the $j$-th image pixel. The signature depends on the tissue’s relaxation times, such as T1 and T2, grouped as a row vector:

[TABLE]

where, $N$ denotes the number of pixels in each image frame. Note that, $j$ is the spatial index while $i$ is the temporal index throughout. Given RF sequence parameters $\boldsymbol{\Theta}^{TRE}$, and parameters $\boldsymbol{\Theta}^{T12}_{j,:}$ of a specific tissue, its temporal signature $\mathbf{X}_{j,:}$ can be derived as:

[TABLE]

where $f(\cdot)$ denotes the Bloch equations. This MR contrast matrix $\mathbf{X}$ is associated with the k-space measurements $\mathbf{Y}$ per column by the subsampled Fourier transform, and it is related to tissue parameters $\boldsymbol{\Theta}^{T12}$ per row by the Bloch equations, as illustrated in Fig. 1.

Given RF sequence parameters $\boldsymbol{\Theta}^{TRE}$ and k-space measurements $\mathbf{Y}$, the goal of MRF reconstruction is to estimate the tissue parameters $\boldsymbol{\Theta}^{T12}$. Typically, the image stack $\mathbf{X}$ is first reconstructed from $\mathbf{Y}$, referred to as signature restoration, and then mapped to tissue parameters $\boldsymbol{\Theta}^{T12}$ via dictionary matching, referred to as parameter restoration 1, 2, 3, 4, 5, 6. This process is illustrated in Fig. 2.

The dictionary is a collection of temporal signatures that are usually simulated by the Bloch equations for various typical tissues, given the pseudo-random RF pulse sequences and tissue parameters. Given an inquiry temporal signature, dictionary matching computes the inner product between the temporal signature with each dictionary entry, selecting the entry in the dictionary exhibiting the highest correlation with the inquiry one as the best matching signature. Once the best entry is found, it directly leads to multiple tissue parameters, such as T1, T2, simultaneously, via searching a lookup-table.

Let $\mathbf{LUT}\in\mathcal{R}^{K\times 2}$ denote a lookup-table composed of $K$ tissues, each containing 2 parameters, i.e., T1 and T2 relaxation times111Note that the off resonance parameter, which appeared in the original MRF paper 1, has been omitted here, since the sequence used in our experiments is derived from the FISP sequence, which is insensitive to off resonance effects 2, 6.. Let $\mathbf{D}\in\mathcal{C}^{K\times L}$ denote the corresponding dictionary simulated using Bloch equations given the RF sequence parameters $\boldsymbol{\Theta}^{TRE}$, formulated as $\mathbf{D}_{k,:}=f(\mathbf{LUT}_{k,:},\boldsymbol{\Theta}^{TRE})$. Since each temporal signature $\mathbf{D}_{k,:}$ is linked with the $k$-th tissue’s parameters $\mathbf{LUT}_{k,:}$, the choice of a large dictionary size $K$ can in principle provide enough granularity to capture a range of possible tissue values.

In conclusion, existing MRF reconstruction approaches involve two stages: signature restoration and parameter restoration, that can be succinctly written as

[TABLE]

where the function $\mathbf{X}=h(\mathbf{Y})$ represents the signature restoration operation such as sparsity or low-rank based de-aliasing and denoising methods, whereas $g(\mathbf{X}_{j,:}|\boldsymbol{\Theta}^{TRE})$ denotes the parameter restoration operation, such as dictionary matching based methods 1, 2, 3, 4, 5, 6.

Our approach aims to perform signature restoration via low-rank based de-aliasing and parameter restoration via a neural network in order to achieve improved MRF reconstruction performance. We highlight that our method only requires a simulated dictionary during network training. Once the network is trained, the dictionary is not needed anymore. In addition, our approach also eliminates a simulated dictionary for signature restoration, which is a key difference from FLOR 5, 6 during signature restoration.

II.B. Previous Methods

Dictionary Matching based MRF approaches. The original MRF reconstruction algorithm 1 is based on dictionary matching, as presented in Algorithm 1. It finds the best matching dictionary entry for the acquired temporal signature according to their inner product and then searches the lookup-table to obtain corresponding tissue parameters. Here, $F_{u}^{H}\{\cdot\}$ denotes the inverse Fourier transform operating on the zero filled k-space data where zeros are filled at the unknown frequencies and symbol $\textbf{Re}\left<\mathbf{a},\mathbf{b}\right>$ represents the real part of the inner product of two vectors $\mathbf{a}$ and $\mathbf{b}$.

Exploiting the nature of signals, by using appropriate prior knowledge, can often contribute to improved signal processing performance. In this spirit, later works suggested to incorporate sparsity in MRF reconstruction to further improve performance, inspired by successful applications of sparsity in MRI reconstruction 9, 10, 11. Davies et al. 3 proposed BLoch response recovery via Iterative Projection (BLIP) which exploits sparsity in the dictionary domain. BLIP consists of iterating between two main steps: (a) a gradient step which enforces consistency with the measurements, based on the Projected Landweber Algorithm (PLA) generalized from the iterative hard thresholding method; (b) a projection which matches each row of $\mathbf{X}$ to a single dictionary atom. Instead of exploiting sparsity in the dictionary domain, Wang et al. 4 suggested to leverage sparsity in the wavelet domain of each imaging frame, $\mathbf{X}_{:,i}$. They further replaced the Euclidean norm with the Mahalanobis distance for dictionary matching. Considering that adjacent MR image frames along the temporal dimension should exhibit high resemblance, Mazor et al. 5, 6 proposed a magnetic resonance Fingerprint with LOw-Rank prior for reconstructing the image stack and quantitative parameters, referred to as FLOR, which achieved state-of-the-art performance. The algorithm, described in Algorithm 2, relies on two priors: a low rank prior on the the matrix $\mathbf{X}$, and the fact that the rows of $\mathbf{X}$ lie in the column space of the dictionary $\mathbf{D}$.

Learning-based MRF approaches. The above techniques all use dictionary matching to perform mapping from temporal signatures to tissue parameters. Therefore, these methods suffer from drawbacks such as discretization error, slow inference speed and memory-consuming storage. In order to alleviate these issues, recent works 13, 14 propose to exploit neural networks to replace dictionaries and lookup-tables used in conventional MRF reconstruction approaches. Cohen et al. suggest a fully-connected feed-forward neural network (FNN) 13. Since the input layer of the FNN is fully connected with the input temporal signature, the number of neurons in the input layer corresponds to the length of the input temporal signature. This makes the network structure less flexible, as a FNN network trained on temporal signatures with a certain length is not applicable to temporal signatures with a different length. In addition, the fully-connected structure results in rapid increase in the number of parameters along with the growth of depth, making the network more susceptible to overfitting. Hoppe et al. 14 propose a 3-layer vanilla CNN for parameter restoration. Both 13 and 14 focus exclusively on learning the signature-to-parameter mapping from a pair of dictionary and lookup-table simulated using the Bloch equations. During the validation, they assume that clean temporal signatures are available as input into the trained networks. However, since temporal signatures obtained from k-space subsampled MRF data are always contaminated by aliasing and noise, their approaches, when applied directly in such k-space subsampling situations, suffer from heavy artifacts introduced during the signature restoration phase, leading to poor performance.

II.C. Proposed Methods

The proposed hybrid deep magnetic resonance fingerprinting (HYDRA) approach, summarized in Algorithm 3, consists of two stages: signature restoration and parameter restoration, (see also (5)). As illustrated in Fig. 3, a low-rank based de-aliasing method is used to restore signatures, and then a 1D nonlocal residual convolutional neural network is used to map each restored signature to corresponding tissue parameters.

In particular, given $\boldsymbol{\Theta}^{TRE}$ and k-space samples $\mathbf{Y}$, in our proposed approach, the function $\mathbf{X}=h(\mathbf{Y})$ in (5) represents a signature restoration operation using low-rank based de-aliasing techniques without requiring a dictionary. The function $\boldsymbol{\Theta}^{T12}_{j,:}=g(\mathbf{X}_{j,:}|\boldsymbol{\Theta}^{TRE})$ in (5) represents a parameter restoration operation that exploits a trained neural network to map each restored signature $\mathbf{X}_{j,:}$ to corresponding tissue parameters $\boldsymbol{\Theta}^{T12}_{j,:}$ directly. In the subsequent sections, we provide a detailed description of both stages of our technique.

II.C.1. Low-rank based signature restoration

Since MRF data consists of multiple frames exhibiting temporal similarity, the imaging contrasts matrix $\mathbf{X}$ is typically a low-rank matrix 6. Therefore, $h(\cdot)$ leverages a low-rank prior for denoising and de-aliasing, formulated as

[TABLE]

where the parameter $r$ is the rank of the matrix, a fixed pre-chosen parameter. Since typically $r$ is not known in advance, we consider a relaxed regularized version:

[TABLE]

where $\|\mathbf{X}\|_{*}$ denotes the nuclear norm 24 of $\mathbf{X}$, defined as the sum of the singular values of $\mathbf{X}$, and $\lambda$ is the Lagrangian multiplier manually selected for balancing data fidelity and the rank. Problem (10) can be solved using the incremental subgradient proximal method 25, similar as to FLOR 6. The procedure for solving (5) is shown in Algorithm 3.

One of differences from FLOR 6 is the fact that we removed the operation of projecting the temporal signal onto a dictionary. This allows to eliminate the requirement for a simulated dictionary in the signature restoration stage, which also alleviates the memory consumption issue. In addition, the computational complexity is reduced by $N\cdot L^{2}$ floating-point operations in each iteration, where $L$ is the dimension of a dictionary element, and $N$ is the number of pixels in each image frame. On the other hand, the gained benefits are at the price of requiring more iterations to converge. Another difference from FLOR 6 is that we exploit a network, instead of dictionary matching, for signature-to-parameter mapping.

II.C.2. Learning-based parameter restoration

Once the imaging contrasts matrix $\mathbf{X}$ is recovered from the k-space samples $\mathbf{Y}$, each temporal signature $\mathbf{X}_{j,:}$ is input into the trained network for parameter restoration, formulated as:

[TABLE]

where $g(\cdot)$ denotes the trained network, $\boldsymbol{\Theta}^{TRE}$ denotes the fixed RF sequence parameters. We next describe the network structure, training and testing procedures.

Network structure.

The proposed network has a 1D nonlocal residual CNN architecture with short-cuts for residual learning and nonlocal operations for achieving a self-attention mechanism. As illustrated in Fig. 3, it starts with two 1D convolutional layers before connecting with 4 residual / non-local operation blocks, and finally ends with a global-average-pooling layer followed by a fully-connected layer. Every residual block is followed by a non-local operation block. Four such blocks are interspersed with each other.

Each residual block contains a max-pooling layer with stride 2, two convolution layers and a shortcut that enforces the network to learn the residual content. The filter size is set to be equal to 21 throughout convolutional layers. The number of channels, a.k.a feature maps, in the first two convolutional layers is set to 16 and then is doubled in subsequent four residual blocks until 128 in the final residual block. The size of feature maps in the next block halves in contrast with the previous one due to max-pooling. In this way, we gradually reduce temporal resolution while extract more features to increase content information.

Inspired by the self-attention scheme and nonlocal neural networks 26, 27, non-local operations are incorporated into the designed network to achieve the attention mechanism, in order to capture long-range dependencies with fewer layers. In contrast to the progressive behavior of convolutional operations that process one local neighborhood at a time, the non-local operations compute the response at a position as a weighted sum of the features at all positions in the feature maps. Formally, the nonlocal operation is formulated as 27:

[TABLE]

Here, $i$ is the index of an output position and $j$ is the index that enumerates all possible temporal positions, $\mathbf{x}$ is the input temporal signal or its features and $\mathbf{y}$ is the output signal of the same size as $\mathbf{x}$. A pairwise function $f$ computes a scalar between $i$ and all $j$ to represent the affinity relationship of these two positions. The unary function $g$ computes a representation of the input signal at position $j$. The response is normalized by a factor $C(\mathbf{x})$. There exists a few instantiations for function $f$ and $g$. For simplicity, the unary function $g$ is chosen as a linear embedding: $g(\mathbf{x}_{j})=\mathbf{W}_{g}(\mathbf{x}_{j})$, where $\mathbf{W}_{g}$ is a weight matrix to be learned. Regarding the affinity matrix $f$, we adopt the embedded Gaussian to compute similarity in an embedding space, which is formulated as: $f(\mathbf{x}_{i},\mathbf{x}_{j})=e^{{\theta(\mathbf{x}_{i})}^{\top}\phi(\mathbf{x}_{j})}.$ Here, $\theta(\mathbf{x}_{i})=\mathbf{W}_{\theta}\mathbf{x}_{i}$ and $\phi(\mathbf{x}_{j})=\mathbf{W}_{\phi}\mathbf{x}_{j}$ are two learned embeddings. The normalization factor is set as $C(\mathbf{x})=\sum_{\forall j}f(\mathbf{x}_{i},\mathbf{x}_{j})$. For a given $i$, $\frac{1}{C(\mathbf{x})}f(\mathbf{x}_{i},\mathbf{x}_{j})$ becomes the softmax computation along the dimension $j$, which leads to $\mathbf{y}=\textit{softmax}(\mathbf{x}^{\top}\mathbf{W}_{\theta}^{\top}\mathbf{W}_{\phi}\mathbf{x})g(\mathbf{x})$, which is the self-attention form 26.

The non-local behavior in (12) is due to the fact that all positions $(j)$ are considered in the operation. As a comparison, a convolutional operation sums up the weighted input in a local neighborhood 27. It implies that the non-local operation directly captures long-range dependencies in the temporal dimension via computing interactions between any two points, regardless of their positional distance. In this way, the network is able to extract global features and take advantage of the full receptive field in each layer.

The global-average-pooling layer is used to average each feature map in order to integrate information in every channel for improved robustness to corrupted input data. This global-average-pooling layer also reduces the number of parameters significantly, thus lessening the computation cost as well as preventing over-fitting. The last fully-connected layer outputs estimated parameters – T1 and T2 relaxation times. The designed network contains around 0.27 million parameters. The weights are initialized using He-normal-distribution 28. The max-norm kernel constraint 29 is exploited to regularize the weight matrix directly in order to prevent over-fitting. The designed network can also be adapted for various MRF sequences, such as the original MRF sequence – inversion-recovery balanced steady state free-precession (IR-bSSFP) sequence, that depends also on the intrinsic df parameter. It is possible to adjust the number of outputs to adapt to more parameters, such as proton density, B0.

To summarize, our network is motivated and inspired by recent successful applications of convolutional neural networks and variants. Convolutional neural networks have been proved to be a powerful model to capture useful features from signals and images. By introducing convolution, local receptive field and weight sharing design, a CNN is capable of taking advantage of local spatial coherence and translation invariance characteristics in the input signal, thus become especially well suited to extract relevant information at a low computational cost 17, 18, 19, 20, 21, 22. On the other hand, the residual network architecture 18, 19 provides an effective way to design and train a deeper model, since it alleviates the gradient vanishing or exploding problems by propagating gradients throughout the model via short-cuts, a.k.a skip connections. By leveraging non-local operation based attention mechanism, neural networks are endowed with capability of extracting global features and capturing long-range dependencies.

Network training.

The designed network is trained on a synthesized dictionary $\mathbf{D}$ and corresponding lookup-table $\mathbf{LUT}$ to learn the signature-to-parameter mappings $\mathbf{LUT}_{k,:}=g(\mathbf{D}_{k,:}|\boldsymbol{\Theta}^{TRE})$.

The training dataset is synthesized as follows. First, we determine the range of tissue parameters. For example, one may set T1 relaxation times to cover a range of [1, 5000] ms and T2 relaxation times to cover a range of [1, 2000] ms with an increment of 10 ms for both. Thus, the T1 and T2 values constitute a grid with dimension 500 $\times$ 200, in which each point represents a specific combination of T1 and T2 values, and hence characterizes a specific tissue. Points corresponding to T1 < T2 have been excluded as such combinations have no physical meaning. All the valid points are stacked together to generate a lookup-table. For instance, the above setting for T1 and T2 leads to a lookup-table of dimension $80100\times 2$. The RF pulse sequences used in our work are fast imaging with steady state precession (FISP) pulse sequences with parameters that have been used in previous publications in the field of MRF 6, 2, 8. Given the lookup-table and RF pulse sequences, dictionary entries can be synthesized by solving the Bloch equations using the extended phase graph formalism 30, 31.

When the training dataset is ready, the dictionary entries are used as input signals and corresponding lookup-table entries serve as the groundtruth. All the dictionary entries are input into the designed network batch by batch which outputs estimated parameters. The root mean square errors (RMSE) of the outputs are calculated with respect to corresponding groundtruth. The resulting RMSE loss is then backpropagated from the output layer to the first layer to update the weights and bias by using Adam 32 as the optimization algorithm. More training details are provided in the subsequent experiment section. Once the training procedure is completed, given an inquiry signal evolution $\mathbf{X}_{j,:}$, it is able to map such a time sequence directly to corresponding tissue parameters, as formulated in (11), implying that no dictionary or lookup-table are required during the inference. Since we only need to store the trained network which is a compact model, it consumes less memory than storing the dictionary and lookup-table.

We emphasize that even though the network is trained on a grid of tissue values, it is expected to capture the mapping function from temporal signatures to tissue parameters. Thus the trained network is capable of outputting tissue values not existing in the grid of training values. Detailed results can be found in Fig. 7 and Table 2. This feature is favorable, as it implies that well designed and trained networks have an ability to overcome discretization issues. The overall procedures for solving (10) and (11) are shown in Algorithm 3.

III. Experimental Results

In this section, we conduct a series of experiments to evaluate our approach, comparing it with other state-of-the-art MRF methods 1, 3, 6, 14, 13.

The experiments are categorized into a few types: training, testing on synthetic data, testing on anatomical data using variable density Gaussian sampling patterns and spiral sampling patterns at different sampling ratios and number of time frames, as described in Table 1. For the network training, synthesized temporal signatures, i.e. simulated dictionary entries of $\mathbf{D}$ shown as Fig. 4, are used as input signals and corresponding parameter values in the lookup-table $\mathbf{LUT}$ serve as the groundtruth. The proposed network is trained to capture the signature-to-parameter mappings. For testing on synthetic data, synthesized temporal signatures in $\mathbf{X}$ are used as input signals and corresponding parameter values in $\boldsymbol{\Theta}^{T12}$ serve as the groundtruth. The aim is to test the parameter restoration performance only. For testing on anatomical data, the k-space measurements $\mathbf{Y}$ which are derived from the Fourier transform (for Gaussian patterns) or non-uniform FFT (for spiral trajectories) 33 of $\mathbf{X}$, are used as input and corresponding parameter values in $\boldsymbol{\Theta}^{T12}$ serve as reference. When there is no k-space subsampling, the aim is to test the parameter restoration performance only. When there exists k-space subsampling, the aim is to test the overall performance, including both signature restoration and parameter restoration. More detailed descriptions are provided in each subsection.

III.A. Training

As mentioned in Section II.C.2., the designed network is trained on a pair of synthesized dictionary $\mathbf{D}$ and lookup-table $\mathbf{LUT}$, simulated using Bloch equations and FISP pulse sequences 2, 6.

The FISP pulse sequence used in our experiments was designed with parameters $\boldsymbol{\Theta}^{TRE}_{:,i}=[TR^{i},TE^{i},FA^{i}]^{T},i\in[1,L]$ that have been used in previous publications in the field of MRF 6, 2, 8. The echo time $TE^{i}$ was constant of 2ms. The repetition time $TR^{i}$ was randomly varied in the range of 11.5 - 14.5 ms with a Perlin noise pattern. All the flip angles $FA^{i},i\in[1,L]$ constituted a sinusoidal variation in the range of 0 - 70 degrees to ensure smoothly varying transient state of the magnetization, as shown in Figure 5.

For the range of tissue parameters, T1 relaxation times are set to cover a range of [1, 5000] ms and T2 relaxation times to cover a range of [1, 2000] ms with an increment of 10 ms for both. Such parameter ranges cover the relaxation time values that can be commonly found in a brain scan 34. All the valid combinations of T1 and T2 values are stacked together, generating a lookup-table $\mathbf{LUT}$ of dimension $K\times 2$ where $K=80100$. Given the lookup-table and RF pulse sequences, dictionary entries are synthesized by solving the Bloch equations using the extended phase graph formalism, leading to a dictionary of dimension $K\times L$ where $L=200$ or $1000$ is the number of time frames.

When the training dataset is ready, the dictionary entries are used as input signals and corresponding lookup-table entries serve as the groundtruth to train the designed network, as mentioned in Section II.C.2.. The model was trained for 50 epochs. It takes around 30 seconds for running one epoch on average, thus around 25 minutes for completing 50 epochs, on a NVIDIA GeForce GTX 1080 Ti GPU. In each training epoch, 20% of the training samples are separated aside for validation dataset. The learning rate decays from 1e-2 to 1e-6 every 10 epochs. Each batch was experimentally set to contain 256 time-sequences in order to balance the convergence rate and weights updating rate well. For comparison purposes, we also implemented Hoppe et al.’s CNN referring to 14, and Cohen et al.’s FNN referring to 13 with the same structure and parameters as specified in their papers. Then we use the same GPU and training dataset to train their networks with specified learning rate and number of epochs until convergence.

We adopt a few widely used metrics, such as root mean square error (RMSE), signal-to-noise ratio (SNR) and peak signal-to-noise ratio (PSNR) to evaluate the image quality quantitatively. The definitions of RMSE, SNR and PSNR are given as follows:

[TABLE]

where matrices $\mathbf{X}$ and $\widehat{\mathbf{X}}$ denote the ground truth signal and its reconstructed version, respectively, $N$ denotes the total number of elements in the signal and $\|\cdot\|_{F}$ denotes the Frobenius norm. PeakVal stands for the pixel peak value in an image, e.g., 1 for a normalized signal.

III.B. Testing on synthetic dataset

In this subsection, we evaluate the performance of HYDRA on a synthetic testing dataset. The procedures of constructing a synthetic testing dataset is similar to the construction of the training dataset: 500 different T1 values are randomly selected from 1 - 5000 ms, while 200 different T2 values are randomly selected from 1 - 2000 ms, using random permutation based on uniformly distributed pseudorandom numbers. All the valid combinations from the selected T1 and T2 values are stacked together, generating a parameter matrix $\boldsymbol{\Theta}^{T12}$ of dimension $80000\times 2$ with $N=80000$. The RF pulse sequences are the same as in the training stage. Given the parameter matrix and RF pulse sequences, input signal signatures are synthesized by solving the Bloch equations using the extended phase graph formalism, leading to a signature matrix $\mathbf{X}$ of dimension $N\times L=80000\times 200$, with each row representing a temporal signature corresponding to a specific combination of T1 and T2 values. The signature matrix $\mathbf{X}$ and parameter matrix $\boldsymbol{\Theta}^{T12}$ constitute the synthetic testing dataset, with $\mathbf{X}$ as input and $\boldsymbol{\Theta}^{T12}$ as the groundtruth.

We input the synthetic testing signatures $\mathbf{X}$ into Hoppe et al.’s CNN 14, Cohen et al.’s FNN 13, and the network of HYDRA to compare the outputs with groundtruth T1 and T2 values in $\boldsymbol{\Theta}^{T12}$. We also compare with dictionary matching methods 1, 2, 3, 4, 5, 6 which exploit the same dictionary $\mathbf{D}$ and lookup-table $\mathbf{LUT}$ to find the best matching entry for each signature in $\mathbf{X}$ and then estimate parameter values by searching the lookup-table. As shown in Table 2, Table 3, Fig. 6 and Fig. 7, the estimated parameter values using the proposed network obtained outstanding agreement with the groundtruth, yielding higher PSNR, SNR and smaller RMSE than the dictionary matching method 1, 2, 3, 4, 5, 6, as well as competing networks 14, 13.

In particular, to illustrate in detail how well neural networks tackle the discretization issue inherent to dictionary matching, we show the testing performance on continuous-valued T1, T2 parameters which have small intervals, e.g. 0.5ms, that is 20 times smaller than the training grid intervals 10ms, between neighboring values in Table 3. Since these values and their corresponding MRF signatures do not exist in the training dictionary and lookup-table222 As mentioned in the experiment setting in Section III.A., in the training dataset, T1 relaxation times are set to cover a range of [1, 5000] ms and T2 relaxation times to cover a range of [1, 2000] ms with an increment of 10 ms for both, that is, T1 values = {1, 11, 21, $\cdots$, 4991}, and T2 values = {1, 11, 21, $\cdots$, 1991}. , the dictionary matching methods report a T1 and T2 value – the closest discretized value present in the dictionary – that can be quite distinct from the groundtruth. In contrast, the various neural network approaches can potentially learn an underlying mapping from the temporal signatures to the respective T1 and T2 values, leading to estimates that are much closer to the groundtruth. Interestingly, our approach outperforms previous networks 14, 13 as shown in Table 3 and Fig. 6. Evidently, neural networks demonstrate much better robustness to discretization issues, leading to improved parameter restoration in comparison to dictionary based methods.

Another impressive advantage of HYDRA is the fast inference speed. HYDRA takes only 8.2 s to complete the mapping operation for eighty thousand temporal signatures, that is, 53$\times$ faster than dictionary matching. Furthermore, the inference speed of HYDRA is subject to the network topology. That is, once the network structure is fixed, the complexity is fixed. In contrast, the complexity of dictionary matching is limited by the dictionary density. This implies that our advantage will be more prominent in comparison with competing techniques using a dictionary with higher density.

III.C. Testing on anatomical dataset

In this subsection, we evaluate our approach on an anatomical testing dataset. We construct the dataset from brain scans that were acquired with GE Signa 3T HDXT scanner from a healthy subject.333The experiment procedures involving human subjects described in this paper were approved by the Institutional Review Board of Tel-Aviv Sourasky Medical Center, Israel. Since there are no groundtruth parameter values for the T1 and T2 parameter maps, we obtain gold standard data by acquiring Fast Imaging Employing Steady-state Acquisition (FIESTA) and Spoiled Gradient Recalled Acquisition in Steady State (SPGR) images, at 4 different flip angles (3∘,5∘,12∘ and 20∘), and implementing corrections 35 followed by DESPOT1 and DESPOT2 36 algorithms. The constructed gold standard T1, T2 parameter maps have a dimension of $128\times 128$ for each map, accordingly leading to a parameter matrix $\boldsymbol{\Theta}^{T12}$ of size $16384\times 2$ by stacking vectorized T1, T2 maps together. Based on the parameter matrix $\boldsymbol{\Theta}^{T12}$ and pre-defined RF pulse sequences, we generate temporal signatures using Bloch equations, the same mechanism as generating the synthetic testing dataset, leading to a signature matrix $\mathbf{X}$ of dimension $N\times L=16384\times 200$. The signature matrix $\mathbf{X}$ and parameter matrix $\boldsymbol{\Theta}^{T12}$ constitute the anatomical testing dataset, with $\mathbf{X}$ as input and $\boldsymbol{\Theta}^{T12}$ as the gold standard reference.

Note that, since the gold standard T1, T2 maps exhibit spatial structures in the image domain, the resulting signature matrix $\mathbf{X}$ can be regarded as a stack of $L=200$ vectorized image frames, where each frame exhibits specific spatial structures. Therefore, it makes sense to perform Fourier transform and k-space subsampling for each column of $\mathbf{X}$ to get k-space measurements $\mathbf{Y}$. This is the key difference between the anatomical dataset and the synthetic dataset.

We first explore the case with full k-space sampling in order to evaluate the parameter restoration performance of HYDRA. Then, we consider situations with k-space subsampling in order to evaluate both the signature restoration and the parameter restoration performance of HYDRA.

III.C.1. Full k-space sampling

In the first case, the fully-sampled k-space measurements $\mathbf{Y}$, derived from the Fourier transform of $\mathbf{X}$, are used as input to obtain the estimated $\boldsymbol{\Theta}^{T12}$. This is equivalent to inputting $\mathbf{X}$ into the network of HYDRA, or performing dictionary matching based on $\mathbf{X}$ directly, since the inverse Fourier transform of the fully-sampled measurements $\mathbf{Y}$ is exactly the same as $\mathbf{X}$. The aim is to test the parameter restoration performance only. In the experiment, corresponding parameter values in $\boldsymbol{\Theta}^{T12}$ serve as the gold standard reference. For comparison, dictionary matching methods 1, 2, 3, 4, 5, 6 exploit the same dictionary $\mathbf{D}$ and lookup-table $\mathbf{LUT}$ as in our training stage to find the best matching entry and estimate parameter values for each signature in $\mathbf{X}$.

Visual and quantitative results are shown in Fig. 8, Fig. 9 and Table 4. It can be seen that our basic version of HYDRA outperforms dictionary matching 1, 2, 3, 4, 5, 6, yielding better visual and quantitative performance, e.g., 7.9 dB SNR gains for T2 mapping. The RMSE of T2 mapping is also reduced to 2.498 from 6.252, accordingly. Our nonlocal version of HYDRA achieves even better performance, leading to 10 dB SNR gains with RMSE as small as 1.86. This is owing to the advantage that the trained network is a powerful function approximator, which is able to output well-estimated parameter values based on learnt mapping functions, even though these values do not exist in the training dictionary and lookup-table. In contrast, dictionary matching only matches signatures to discrete parameters existing in the training dataset. In other words, if there are no exact matching dictionary element and parameter values for an inquiry MRF signature, it will find adjacent values as approximations, thus introducing discretization error. On the other hand, the advantage of HYDRA over dictionary matching on T1 mapping is not as significant as on T2 mapping quantitatively. But the visual improvements are evident. A similar trend is observed when comparing our network with competing networks such as Hoppe et al.’s CNN 14 and Cohen et al.’s FNN 13. In addition, HYDRA takes around 2 s to accomplish the mapping for 16384 signatures, 40$\times$ faster than dictionary matching 1, 2, 3, 4, 5, 6.

III.C.2. k-space subsampling using Gaussian patterns

In k-space subsampling situations, the developed low-rank based de-aliasing method is applied to restore the signature matrix $\mathbf{X}$ from the measurements matrix $\mathbf{Y}$. Then, the reconstructed $\mathbf{X}$ is used as input into the network for parameter mapping to obtain the corresponding tissue parameter values. In the experiments, the sub-sampling factor $\beta$ is set to be 70% and 15%. For $\beta$ = 15%, 15% k-space data is acquired by a series of 2D random Gaussian sampling patterns, shown in Fig. 10, leading to a k-space measurement matrix $\mathbf{Y}$ of size $Q\times L=2458\times 200$. Similarly, $\beta$ = 70% gives rise to a k-space measurement matrix $\mathbf{Y}$ of size $Q\times L=11469\times 200$. A larger $\lambda$ enforces lower rank for the restored signature matrix $\mathbf{X}$ to strengthen the de-aliasing effect, while a smaller $\lambda$ encourages $\mathbf{X}$ to have a subsampled Fourier transform that approximates the k-space measurements matrix $\mathbf{Y}$ better. Therefore, we tried a range of values from 1 to 20 for $\lambda$ and experimentally select the best one $\lambda=5$. Since the low-rank based signature restoration involves gradient descent steps, a larger step size $\mu$ accelerates gradient descent speed, but tends to result in oscillation or even divergence, while a smaller $\mu$ leads to a slower convergence. We experimentally find that $\mu=1$ gives a good balance. The same k-space measurements $\mathbf{Y}$ are also used by dictionary-matching based methods 1, 3, 6 for comparison, and the same signature restoration approach is used to convert $\mathbf{Y}$ onto $\mathbf{X}$ for learning based methods 14, 13. The aim is to evaluate the overall performance on both signature restoration and parameter restoration.

Quantitative performance is shown in Table 5. Note that the advantage of learning-based methods over dictionary matching degrades when the subsampling factor increases. This is due to the fact that the restored signatures from highly subsampled k-space data exhibit deviations and distortions, thus leading to poorer input for the trained networks. In spite of this, the proposed approach outperforms the dictionary matching based methods 1, 3 with significant gains, and also yields better or comparable performance as the state-of-the-art methods FLOR 6, CNN 14 and FNN 13. In addition, it takes around 23s for low-rank based signature restoration and less than 3s for network based parameter restoration. Thus, the total time cost is around 26s, almost 4.8$\times$ faster than FLOR 6. Furthermore, the speed of our method is 60$\times$ faster than FLOR 6 for parameter restoration.

We compared the performance with/without nonlocal operations in our developed network. The results in Table 4 and 5 show that the proposed network with nonlocal operations based self-attention scheme outperforms the basic counterpart. In particular, the nonlocal version achieves 6 dB gains in terms of SNR over the basic version for T2 mapping. Such significant improvement demonstrates the benefits of capturing long-range dependencies and global features using the nonlocal operation based attention scheme.

We also investigated the performance with respect to the number of time frames. In particular, we increased $L$ from 200 to 1000 and kept other experiment settings the same as before. The quantitative results are shown in Table 6. It is noticed that given more time frames, all the methods show better performance. Moreover, the performance of learning-based methods, including CNN 14, FNN 13 and HYDRA, improve more than model-based techniques 1, 3, 6. In particular, our approach outperforms competing algorithms quantitatively in terms of PSRN, SNR, and RMSE, as well as demonstrates visual advantage, as shown in Figure 11 and Figure 12.

III.C.3. k-space subsampling using spiral trajectories

We carried out additional experiments with widely used non-Cartesian sampling patterns – variable density spiral trajectories 37, 6. A set of spiral trajectories used in the experiments are shown in Figure 13. They have FOV of 24 and rotation angle difference of 7.5 degrees between any two adjacent spirals to spread out the alias artifacts. Given such spiral trajectories, data were subsampled to acquire 1488 k-space samples in each time frame, leading to a subsampling ratio of 9% which is defined by the number of acquired samples in the k-space domain divided by the number of pixels in a frame. This setting closely matches the original MRF paper 1 where each single spiral trajectory samples 1450 k-space points (leading to a subsampling ratio around 9%) and any two adjacent spiral trajectories have a rotation angle of 7.5 degrees.

In the case of spiral subsampling, during the signature restoration, SParse Uniform ReSampling (SPURS) algorithm 38 was exploited to implement nonuniform Fourier transform between k-space domain and image domain, as SPURS has proved to achieve smaller approximation errors while maintaining low computational cost comparing with other resampling methods, such as nonuniform-FFT algorithm 33 and regularized Block Uniform ReSampling (rBURS) 39. In the experiments, 1000 density variable spiral trajectories were used for k-space subsampling, leading to 1000 time frames. The quantitative and qualitative reconstruction results demonstrate that HYDRA outperforms competing methods with smaller estimation errors, as shown in Table 7, Figure 15 and Figure 14.

IV. Discussion

IV.A. Relation to previous works

Our low-rank based signature restoration method is adapted from FLOR 6 by removing the operation of projecting the temporal signal onto a dictionary. Thus, the signature restoration does not require a simulated dictionary, and saves computational cost. Although recent works 13, 14 exploit neural networks to perform parameter mapping, replacing dictionaries and lookup-tables used in conventional MRF reconstruction approaches, our technique is different from these methods 13, 14. We design a deep nonlocal residual CNN for capturing signature-to-parameter mapping which is organically combined with low-rank based de-aliasing techniques for signature restoration. In this way, our algorithm can bypass some of the issues associated with other techniques: (1) The input dimension issue. The proposed approach can ingest temporal signatures with different lengths without the need to change the structure of the network. This is due to the fact that we rely on convolutional neural networks (CNNs) rather than fully-connected neural networks (FNNs) such as the model used in 13. (2) The k-space subsampling issue. The proposed approach involves a hybrid of a neural network with a low-rank based de-aliasing approach. Thus it is able to deal with correlations both over time and space via exploiting low-rank regularization and convolution operation. This enables our work to handle k-space subsampling situations. (3) The complex mappings issue. By exploiting a residual network structure, our method can be successfully extended to deeper levels and thus obtain a better capacity to learn complex signature-to-parameter mapping functions. (4) Distortion and corruption issue. Due to the subsampling in k-space, the restored temporal signatures suffer from local distortion and corruption. Such deviation may lead to performance degeneration in the second stage. By incorporating non-local operations in the network design, our method is able to capture global features and find most relevant components for inference, thereby reducing interference of local distortion and corruption.

IV.B. Computational complexity

HYDRA involves two main stages: the low-rank based signature restoration stage and the network based parameter restoration stage. Even though the time cost for parameter restoration is longer than previous methods 13, 14, the time cost in the this stage is only a small fraction of the total time consumption, as the computational complexity is dominated by the signature restoration stage. In other words, the computational burden of HYDRA lies in the SVD calculation in the first stage. Hence, fast SVD methods can be employed to dramatically improve the efficiency of signature restoration.

IV.C. Model storage requirements

Regarding the storage requirement (in double precision), HYDRA needs only 2.1 megabytes to store the network with 0.5 million parameters, while it requires 108 megabytes to store a simulated dictionary of size $80100\times 200$ and 551 megabytes for size $80100\times 1000$. Note that the dictionary volume will grow exponentially with the number of parameters, but the space required for storing a network is not strictly limited by the dictionary density once the topology of the network is fixed, thus significantly alleviating the storage burden inherent to the exponential growth of multi-dimensional dictionaries.

IV.D. Impact of providing continuous T1/T2 values

Providing continuous T1/T2 values is an advantage of neural network based parameter mapping over dictionary matching. This property may find promising applications in some practical scenarios, for exampling, monitoring sensitive changing of pathology condition over time, such as multiple sclerosis 40, 41, stroke 42, and treatment responses 43, 44, where the differences in T1 and T2 values between healthy and diseased tissues or between disease stages could be very small 45. On the other hand, to fulfil this potential of network based MRF techniques, prerequisites on the accuracy and precision of MRI measurements are needed. Taking T1/T2 quantification as an example, even for the inversion recovery spin echo (IR-SE) / multiple single-echo spin echo MRI sequences which are considered as the gold standard for T1/T2 quantification, there exist variations of 2% - 9% on the measured relaxation times 45. Such anatomical measurement uncertainties and model imperfections may weaken the advantage and clinical impact of providing continuous T1/T2 values using network based MRF techniques to some extent. Therefore, improving the accuracy of gold standard approaches in the future would contribute to making the most of the potential of neural networks in the MRF domain.

V. Conclusion

We proposed a hybrid deep MRF approach which combines low-rank based signature restoration with learning-based parameter restoration. In our approach, a low-rank based de-aliasing method is used to restore clean signatures from subsampled k-space measurements. Then, a 1D deep nonlocal residual CNN is developed for efficient signature-to-parameter mapping, replacing the time-consuming dictionary matching operation in conventional MRF techniques. Our approach demonstrates advantages in terms of inference speed, accuracy and storage requirements over competing MRF methods as no dictionary is needed for recovery.