Systematic review of reconstruction techniques for accelerated quantitative MRI

To systematically review the techniques that address undersampling artifacts in accelerated quantitative magnetic resonance imaging (qMRI).


INTRODUCTION
Quantitative MRI (qMRI) measures tissue properties such as diffusion, relaxation times, myelin water fraction, and perfusion. 1 qMRI produces maps of calibrated physical parameters expressed in physical units that are reproducible across the scanner. 2 Conventionally, qMRI reconstruction consists of a two-step procedure. The first step reconstructs multiple images from the acquired k-space data. For each image to be reconstructed, the corresponding k-space data is acquired with unique acquisition settings such as echo time, repetition time, or diffusion weightings. These different acquisition settings lead to differently weighted reconstructed images, that is, images with different contrasts. In what follows, these images will be referred to as contrast images. The second step estimates the qMRI parameters by fitting a signal model to the reconstructed contrast images.
Acquisition of multiple contrast images traditionally requires a long scan time, which limits the clinical application of qMRI. Reducing this time has been an important research focus over recent years. Sampling only part of the k-space for each of the contrast images to be reconstructed, that is, undersampling, is one way to shorten the scan time. However, substantial undersampling of the k-space leads to ill-posed reconstruction problems for the individual contrast images, causing artifacts in the reconstructed images.
Acquiring data from the same object with different contrasts produces shared information across the contrasts that can be exploited to improve the qMRI reconstruction of the undersampled data. The different qMRI reconstruction techniques differ from each other in the way they exploit shared information across the contrasts. The multitude of reconstruction techniques and their respective advantages and limitations give rise to the question of to what extent they are different and whether some of these techniques can be advantageously combined. Moreover, as techniques are typically presented for specific applications, the suitability to different applications can be hard to appreciate. Furthermore, the reconstruction techniques are usually introduced using different terminologies, which complicates their comparison. Review studies on qMRI reconstruction techniques are typically specific to a class of techniques such as compressed sensing, 3 MR fingerprinting [4][5][6] or deep-learning. 7,8 Alternatively, they are specific to a particular class of parameters of interest, such as relaxometry, 7,[9][10][11][12] or diffusion 13 parameters, or to a particular application domain, such as the knee, 14 or cardiac and abdominal imaging. 15 To facilitate the comparison among different reconstruction techniques used for accelerated qMRI, the current work presents a systematic review that provides a technical overview of the existing techniques, without confining to any specific parameter or application domain, and categorizes them based on methodology. Moreover, all categories are described in a unified mathematical framework.
This systematic review was conducted following the PRISMA 2020 guidelines. 16 The methodology used for the systematic review is introduced in Section 2. Next, Section 3 describes the mathematical formulation of conventional two-step qMRI. Section 4 provides the mathematical formulation of the categorized reviewed techniques and describes how each technique is adjusted to address specific reconstruction issues. In Section 5, a graphical representation of the categorized techniques is provided. Finally, Section 6 provides a discussion of the trends observed in the reviewed techniques, limitations of this work, and directions for future research.

METHODS
The systematic review was performed in four phases, which are discussed in the four subsections below and depicted in Figure 1.

Identification
A literature search was conducted in the following databases: Embase, Medline, Web of Science Core Collection, Coherence Central Register of Controlled Trials, and Google Scholar. The search query was compiled with the help of a librarian specialized in performing systematic reviews. The study was registered in Prospero with id CRD42021242450. The search query used was based on a combination of four main segments, combined with an AND statement. The first segment included words related to "MRI," "MR fingerprinting," and "quantitative MRI." The second segment included words associated with different "reconstruction" techniques to select the papers that deal with image reconstruction. The third segment referred to the qMRI tissue parameters such as T 1 , T 2 , and diffusion parameters. The last segment was designed to select the works that deal with image reconstruction and parameter estimation. The exact search query for each database is provided in Table A2 in the Appendix. For the search conducted in Google Scholar, the first 200 results were initially selected. Among those 200 results, the ones that were not found previously were added to the results of the other databases. The final search was performed on July 7, 2022 and was confined to articles published in English, without restriction on the publication date. This led to the identification of 4190 papers.

Screening
The search results were examined based on their title and abstracts by two authors (Banafshe Shafieizargar and Riwaj Byanju) independently. Papers were selected when meeting all of the following inclusion criteria: 1. Focusing on qMRI parameter estimation.

2.
Accelerating data acquisition by undersampling the k-space. 3. Exploiting the joint information between contrasts to avoid undersampling artifacts. 4. Novel reconstruction and/or parameter estimation technique or improvement of an existing technique other than improvement based on better hardware, in comparison to conventional qMRI reconstruction.

F I G U R E 1
Flow chart showing the selection of articles.

F I G U R E 2
Categorization of the techniques reviewed in the current study.

5.
Including a detailed description of the proposed image reconstruction and/or parameter estimation technique.
If it was unclear whether a paper met a particular criterion from the information in the abstract, then the criterion was considered fulfilled. This phase led to the selection of 391 papers.

Full-text reading
The full-text of each paper selected during the screening phase was read by one of the two authors who also performed the screening (Banafshe Shafieizargar and Riwaj Byanju). If during or after full-text reading it was concluded that not all criteria of Section 2.2 were met, the paper was excluded. This led to the exclusion of 99 papers so finally, 292 papers were selected for the current study.

Categorization
In this final phase, the selected works were categorized based on the year of publication, parameter of interest, methodology, application domain, and accessibility (i.e., a public link to the implemented code and/or acquired dataset). The methodologies proposed in the selected papers were categorized according to several properties, which are listed in Figure 2. These properties will be further defined in Sections 3 and 4.

CONVENTIONAL QMRI RECONSTRUCTION
This section describes the mathematical formulation of conventional qMRI, introducing a notation that will be used throughout this work. qMRI reconstruction aims at the estimation of parameter maps from measured k-space data and is conventionally performed in two steps: "image reconstruction" and "parameter estimation." To facilitate reading, a list of symbols used in this work is provided in Table A1 in the Appendix.

Image reconstruction
In a multicontrast multichannel experiment with n c contrasts indexed by c ∈ {1, … , n c } and n s coil sensitivity maps indexed by s ∈ {1, … , n s }, the measured MRI signal Y c,s (v) at a point v ∈ Ω, where Ω ⊂ R 3 denotes the three-dimensional (3D) Fourier space, also known as k-space, can be described as: whereĨ(r) represents the transverse magnetization at spatial positionr ∈ R 3 in image space,C s (r) is the coil sensitivity of the sth coil at positionr, andẼ is a zero-mean complex-valued Gaussian noise contribution. Note that, since the measurementsỸ c,s (v) are made in discrete space, the reconstruction of the continuous magnetizationĨ(r) is ill-posed. In practice, however, images of the transversal magnetization are reconstructed on a finite 3D grid {r } n 1 with r ∈ R 3 , ∀ ∈ {1, · · · , n } where n denotes the number of grid points. In this case, the relationship described by Equation (1) can be approximated by a discrete Fourier Transform: Y ∈ C n c ×n s ×n k represents the acquired k-space points with n k points sampled for each contrast and coil channel, indexed by k ∈ {1, · · · , n k }, C ∈ C n s ×n represents the coil sensitivity maps, I ∈ C n c ×n denotes the contrast images, and N ∈ C n c ×n s ×n k is a zero-mean, complex-valued Gaussian noise contribution. For an undersampled acquisition, the set of acquired k-space points does not satisfy the Nyquist criteria. Note that Equation (2) can be written using an encoding operator  that combines the coil sensitivity maps and the (nonuniform) Fourier Transform for the sampled points: where Y, I, and N are (implicitly) reshaped into vectors when operators are applied. Note that the operator can be applied to both Cartesian and non-Cartesian (such as spiral or radial) sampling trajectories by only changing the set of k-space locations (v). Conventionally, for a fully sampled multichannel acquisition, contrast images are reconstructed by applying an inverse discrete (or density compensated nonuniform) Fourier transform to the k-space data of each coil channel, followed by a combination of images corresponding to each coil channel. 17 Application of this technique to undersampled k-space coil data may produce aliased images. However, for suitable undersampling patterns, it is possible to reconstruct artifact-free images using an iterative approach through parallel imaging (PI) reconstruction. To ease the notation below, we use a named index to select an element in the corresponding dimension of an object, without explicitly listing the other dimensions. So, I c ∈ C n is a vector with the n voxels in the c th contrast. Assuming the image domain coil sensitivity maps are available, the PI image reconstruction, also known as SENSE, 18 can be formulated as follows: In case the explicit coil sensitivity maps are not available, auto-calibrating PI techniques such as GRAPPA 19 are performed to reconstruct individual contrast images. In GRAPPA, the missing samples in undersampled k-space data are predicted using a linear combination of the acquired neighboring k-space data across all the coils. A fully sampled patch of k-space data, also known as a calibration region, is used to learn the relationship between the acquired and missing samples using reconstruction weights, known as kernels. The fully sampled k-space is then followed by the Fourier transform and a combination of images corresponding to each coil channel. The image reconstruction approaches described in this section are referred to as FTc.

Parameter estimation
The second step of the conventional qMRI reconstruction approach is parameter estimation using model fitting.
In this technique, n q tissue parameter maps, represented by X ∈ R n q ×n , are estimated from the contrast images. Let X ∈ R n q represent the vector containing the tissue parameters for the jth voxel and I ∈ C n c is a vector that contains the values of the n c different contrast images, referred to as signal evolution, in the th voxel. Next, let the relationship between the contrast images and the underlying tissue parameters of interest in each voxel be modeled as: with  (X) ∶ R n q ×n → R n c ×n a parametric function that represents a qMRI model such as a relaxometry, diffusion, or perfusion model. Note that this function will not only depend on X, but also on the image acquisition settings, which are assumed to be known.
In the parameter estimation step, the tissue parameter maps X can be estimated by fitting the parametric model to the reconstructed contrast images using a goodness-of-fit criterion, such as the sum of squared differences, least squares criterion, or the (penalized) likelihood function. [20][21][22][23] This process can be denoted by  −1 . When  evaluates each voxel independently, the tissue parameter maps can be estimated voxelwise: whereÎ refers to the reconstructed images obtained in the image reconstruction step. Note that depending on the application, the parameters can be estimated from reconstructed complex-valued images (e.g., flow parameter estimation 24 ), or reconstructed magnitude images (e.g., conventional diffusion tensor imaging 25 ).

RESULTS
The papers considered in this review aim at estimating parameter maps X from undersampled k-space data Y. These works either modify the conventional image reconstruction step (Section 3.1), the conventional parameter mapping step (Section 3.2), both, or directly estimate qMRI parameters X from undersampled Y. Works that include an intermediate image reconstruction step are categorized as indirect reconstruction, whereas techniques that skip the image reconstruction step and estimate parameters directly from k-space data are categorized as direct reconstruction. Both indirect and direct reconstruction are further divided into different subcategories as shown in Figure 2.
In the current section, for each subcategory, first, a theoretical overview is provided that describes the basic technique used in that subcategory. Next, the advancements proposed in the reviewed papers within this subcategory are described. Note that this is not necessarily a complete literature review on such subcategory. The complete list of the papers that were selected in the categorization phase, in combination with the properties that describe them, can be found in Table A3 in the Appendix and a more detailed version of this table with additional information on the acquisition sequence and trajectory, and optimization solvers is available online.

Indirect reconstruction
This category presents techniques that follow a two-step procedure of image reconstruction followed by parameter estimation.

Image reconstruction
This section categorizes the techniques that jointly reconstruct all the contrast images.

Multicontrast PI
The multicontrast PI subcategory includes PI techniques that have been extended to include multicontrast redundancy. Both image and k-space domain-based techniques have been developed and are described below: Image domain PI. In SENSE reconstruction, the coil sensitivity maps can be estimated from a separate scan or a calibration region. Alternatively, studies have been performed to eliminate the need to acquire a calibration region by jointly estimating the coil sensitivity maps with the contrast images. 26,27 In this case, the iterative reconstruction described by Equation (4) is formulated as a nonlinear inverse problem with both I c and C as unknowns. This problem has then been solved by enforcing smoothness of the coil sensitivities in a regularized reconstruction, 26 or using kernel calibration. 27 Multicontrast SENSE reconstruction can be combined with reconstruction techniques for other acceleration techniques such as simultaneous multislice imaging (SMS), where multiple slices in the field-of-view are acquired simultaneously by using multi-band RF pulses 28 or the "UNFOLD" method that involves sampling of different k-space points from contrast to contrast and applies filters along the contrast dimension. 29 k-space domain PI. K-space domain-based PI techniques, such as GRAPPA, have been extended to include the prediction of missing k-space samples from data with different contrasts. [30][31][32][33][34][35][36][37] Additionally, GRAPPA has been applied in undersampled 3D acquisitions to reconstruct the missing k-space points along the partition encoding dimension. 38 Furthermore, GRAPPA has been extended to include SMS acquisitions. In this case, GRAPPA is used to unwrap the intermixed signal of the simultaneously acquired slices. [39][40][41][42][43][44][45] Finally, GRAPPA has been followed by an iterative reconstruction where distortion and noise in the reconstructed images are reduced. 46

Regularized reconstruction
Techniques in the regularized reconstruction subcategory use prior knowledge of the sparsity in specific domains to improve image reconstruction. Regularized reconstruction can be formulated based on Bayes' theorem as a combination of the likelihood function and prior knowledge in the form of regularizer terms. Assuming the operator transforms the images to be reconstructed to a desired sparse domain, then the reconstruction problem with n i sparsifying transforms, indexed by i, can be formulated as the following optimization problem:Î = arg min where i is the regularization weight of the ith sparsifying transform and l ∈ N defines the regularization norm. Note that, the transform operator can be applied across spatial as well as contrast dimensions. Transform domain. Finding a suitable transform domain where the signal is sparse has been one of the challenges in performing regularized reconstruction. Different transform operators and priors have been proposed to find a sparse representation of the images. This includes the wavelet transform, 27,47-52 the total variation transform, 26,50,52-68 group sparsity where images are divided into multiple sparse regions, 69 weighted quadratic prior that aims to suppress the noise and reconstruction artifacts based on the intensity differences between neighboring voxels, 56 gradient across the contrast dimension, 53,70-74 second-order discrete derivative in the contrast dimension, 75,76 principal component analysis-based transform, 75,[77][78][79][80] image ratio constraints, where the ratio between a low-resolution image and the reconstructed image is used as a constraint, 50 and learned sparsifying transform from the measurements. 81 Apart from these, alternative ways to use regularizers and transform domains have been proposed. First, using a mixed l 1 − l 2 norm through distributed compressed sensing has been proposed that promotes sparsity along the contrast and spatial dimensions simultaneously. 52, 82,83 Secondly, priors based on the distribution of the contrast images, for example, a Gaussian process prior are proposed. 84 Finally, the reconstruction problem can be constrained using a parametric model: 59 Low-rank priors. The matrix of contrast images I, in which each row is a vectorized contrast image and each column shows the signal evolution for one voxel, can be assumed to be a low-rank matrix (referred to as Casorati low-rankness). Such assumptions can be incorporated as prior information by formulating the regularized reconstruction in the form of a rank constraint: where defines the weight of the rank constraint. However, due to the rank constraint, the optimization problem described by Equation (9) is not convex. Therefore, the nuclear norm of the matrix of contrast images ||I|| * , that is, the sum of its singular values, has been often used as convex relaxation of the rank function, 61,79,80,82,[85][86][87][88][89][90][91][92][93][94][95][96][97][98][99] giving the reconstruction problem:Î = arg min The nuclear norm regularization can also be applied to phase data in applications where phase correction needs to be integrated into the image reconstruction. 100 Instead of the nuclear norm, a log-sum relaxation form of the rank constraint based on the weighted nuclear norm minimization has been proposed to simplify the optimization of Equation (9). 101 As another option, a patch-based low-rank constraint can be used. 52,102-106 Let  be an operator that vectorizes a patch of contrast images around voxel , corresponding to the selection of a set of columns of I, then the reconstruction problem can be formulated as: Combination with other techniques. Regularized reconstruction has been used in the Generalized SLIce Dithered Enhanced Resolution technique, which is a newly developed RF encoding technique. 107 In this technique, the reconstruction has an extra processing step where the super-resolution reconstruction model relates the acquired thick-slice RF-encoded volumes to the desired thin-slice volumes. A joint regularized optimization is performed on the images that are initially reconstructed using k-space-based PI techniques. 36,43,108 Different regularizers have been used in the extra enhancement step, including the Huber function that preserves the sharp edges of the contrast images, 108 the Tikhonov regularizer, 36 and the finite difference transform. 43 Optimization algorithms. Solving the combined optimization of data fidelity and regularizer terms can be difficult and often a more relaxed version of the cost function has been used which is computationally easier to optimize. One of the commonly used algorithms in optimization of the regularized cost functions is the Alternate Direction Method of Multipliers. 52,53,57, 69,93,95,100,101,103 Other techniques include the gradient descent algorithm, 72-74 the Fast Composite Splitting Algorithm 61,82 and majorization-minimization-based algorithms such as the fast iterative shrinkage-thresholding algorithm. 106

Subspace constrained reconstruction
The central assumption behind papers in the subspace-constrained reconstruction subcategory is that the signal evolution across the contrasts lies in a low-dimensional subspace. 109 This differs from the low-rank prior introduced in Section 4.1.1.2 by explicitly, and often a priori, specifying a rank for the signal evolution subspace and often also a subspace basis. In this case, the matrix I can be approximated with a linear relationship as I ≈ , with ∈ C n c ×n b containing the n b basis functions, and ∈ C n b ×n their voxel-wise weights. The image reconstruction problem can then be simplified by estimating voxel-wise weights instead of contrast images I as follows: Subspace computation. The subspace is typically obtained by singular value decomposition of a dictionary of possible signal evolutions.
Model-based dictionaries D ∈ C n c ×n d consist of n d dictionary items, each with the expected signal evolution for the given acquisition settings and specific qMRI parameters.  More details about the generation of the dictionaries are presented in Section 4.1.2.1.
It has been shown that using patches of I to create subspace functions, improves the compression of the dictionary by accounting for local similarities among signal evolutions. [137][138][139][140][141][142] Furthermore, the subspace initially generated using the model-based dictionary, can be updated in an iterative reconstruction. In each iteration of this reconstruction, the subspace is updated using the reconstructed images. [143][144][145] Additionally, to better represent the nonlinearity in the signal evolution, it has been proposed to generate the subspace using kernel-based principal component analysis, in which the subspace is learned from low-resolution training data. [146][147][148] Alternatively to dictionary generation, the subspace basis functions can be estimated along with the voxel-wise weights during the reconstruction step by using the "blind compressed sensing" principle: 149 where || || 2 F < 1 represents a unit Forbenius norm applied to make the optimization problem well-posed. The optimization (13) can be solved using a variable splitting and augmented Lagrangian optimization scheme. 149 Multiple subspace basis functions. During data acquisition, subject motion, including rigid, respiratory, and cardiac motion, can influence the acquired signal. To account for subject motion, some techniques based on partially separable functions have been proposed. In these techniques, different subspace basis functions related to the contrasts and motion can be separately designed and used during the reconstruction. 145,[150][151][152][153][154][155][156] Specifically, the contrast images I can be decomposed into multiple tensors using a higher-order singular value decomposition, such as the Tucker decomposition. For instance, I is separated into components representing the signal evolution across contrastš, and motion as: where m is the index for the motion states, b 1 and b 2 indices the basis functions of̌∈ C n c ×n b 1 and ∈ C n m ×n b 2 , respectively, anď∈ C n b 1 ×n b 2 ×n represents a tensor that contains spatial weights for each basis. 145,151,156,157 During the reconstruction the spatial weightšare estimated. Equation (14) is limited to two bases but this number can be increased to include more aspects of the signal evolution. 130,152,[158][159][160][161] By including SMS in the acquisition protocol, the scan time for such techniques can be further reduced. 160 Along similar lines, it has been proposed to separate the spatial weight into a magnitude and a unit phase component that represents the phase variations. These weights are then estimated using a block coordinate descent method. 162 Additional steps to address other aspects of the acquisition. First, subspace constrained reconstruction has been modified to correct for phase evolution across the readout jointly with image subspace reconstruction. 136,163 Secondly, it has been used for compensation of B0 inhomogeneities in echo planar time imaging sequences. 116 Finally, it has been proposed to separate fat and water in subspace constrained reconstruction. 117,138 Additional prior in the form of a Regularizer. Regularizers can be introduced to add prior information, such as sparsity, 164 to the optimization problem defined in Equation (12), yielding: Subspace constrained reconstruction has been used jointly with a spatially sparse representation of using the wavelet transform, [110][111][112]119,120,163 total variation, 117,119,123,124,129,143,145,155,162,165 finite difference, 164 Fourier representation of images in the temporal direction, 165 low-rank priors, 112,115,126,128,[135][136][137][138][139]141,142,145,[166][167][168][169][170][171][172] and deep learning priors. 173 Alternative ways to exploit the low-rank nature of the signal. The low-rank nature and sparsity of the signal evolution can be exploited by assuming it can be decomposed into its low-rank component L and sparse component S, that is, I = L + S. 60,[174][175][176][177][178] This results in the reconstruction problem: The low-rank subspace of the signal has been used to fill the missing k-space samples. 113,[179][180][181][182][183][184][185] Specifically, the subspace Ψ constructed in the image domain can be used to fill in missing k-space samples. 113 The k-space data can also be transformed into structured matrices, for example, block-Hankel 181,182,[184][185][186] and Toeplitz, 183,187 to exploit the low-rank properties. In this case, the optimization problem can be rewritten as: where  is an operator that first converts the contrast images into the Fourier domain and then into a structured matrix such as block-Hankel, or Toeplitz. The reconstruction problem has been originally posed with a minimization constraint on I, min rank[I], however, since this is an NP-hard problem, it is replaced with a Schatten-p quasi norm (0 ≤ l < 1) 187 or the nuclear norm. 181 It should be noted that reconstruction techniques based on structured low-rank matrices are computationally demanding due to the large size of the structured matrices. To address this issue, a more computationally efficient implementation of these techniques based on modified iterative reweighted least squares has been proposed. 180 To improve the reconstruction further, after solving the reconstruction problem described in Equation (17), a convolutional neural network with a U-NET architecture has been used. 181,188 Another way to exploit the low-rank nature of the contrast images is by adding a regularization term penalizing the distance of the images to the subspace. 47,62,63,[189][190][191][192][193][194] In this case, the optimization problem can be defined as: Optimization algorithms. Most of the works use alternate direction method of multipliers with variable splitting for solving subspace constrained reconstruction problems. 112,115,116,118,121,122,126,128,130,[139][140][141][142]144,149,[162][163][164][165][166]169,183,187 Other solvers include the nonlinear conjugate gradient descent with backtracking line search, 120 the orthogonal matching pursuit method, 143,190 the linear conjugate gradient algorithm 123 and the fast iterative shrinkage-thresholding algorithm. 135

View-sharing reconstruction
Acquisition protocols have been proposed that create densely sampled k-space frames by combining the sampled points of different subsampled contrasts. While reconstruction of the individual contrast images from their corresponding subsampled k-space data results in aliased or low-resolution images, a smart combination of the multicontrast k-space data can reduce these issues. However, combining k-space points of different contrasts in one k-space may lead to artifacts in image space due to the difference in contrast across the k-space. Thus, additional steps, such as filtering, may be required. 195 View-sharing with shared high spatial frequency information. The principle behind techniques in this subcategory is that the changes in the acquisition of one contrast to the other mainly affect the low spatial frequencies in the central k-space regions. Therefore, while the central regions are fully sampled for every contrast, the noncentral regions, containing the high spatial frequencies, are undersampled. For each contrast, a k-space frame is formed by complementing the fully sampled central region with high-frequency data of multiple contrasts, after which these k-space frames are reconstructed into images. [196][197][198][199][200][201][202] The techniques in the view-sharing subcategory have been further improved by applying additional filters to help maintaining the desired information. For instance, "k-space weighted image contrast" filters are used that help reducing streaking artifacts in radial sampling acquisitions. [199][200][201][202] Similarly, the application of multiple radial mask filters is examined that exploit the low spatial-frequency nature of image-to-image contrast changes. 195 View-sharing with information shared across similar contrasts. Techniques in the view-sharing subcategory can also assume a slow variation among the contrasts such that the data from several similar contrasts can be combined into a single k-space and reconstructed into a single weighted image, thereby ignoring the contrast variations. [203][204][205][206][207][208][209][210] However, additional correction steps are required to reduce the effect of the small differences between the contrasts. Specifically, the trajectory or the order in which the k-space is sampled can be crucial to avoid artifacts.. 203 View-sharing can be combined with other reconstruction techniques such as regularized reconstruction, 66 subspace constrained reconstruction, 211 or PI reconstruction. 41,206,210 Learning-based reconstruction Techniques in the learning-based reconstruction subcategory use neural networks that learn to reconstruct artifact-free images. In this section, we focus on the works that exploit the redundant information in contrast images in the training process of the neural network. Learning-based image reconstruction is often seen as an artifact removal step where the input of the network is the aliased images that are reconstructed using the FTc techniques as described in Section 3.1 and the output is the artifact-free reconstructed contrast images. A neural network(Ž| ), which is parameterized by , is provided with aliased contrast images Ž as input, to provide artifact-free contrast imagesÎ as output: [212][213][214] The parameters of the network are defined in a training process, prior to image reconstruction, that can be described as:̂= where  denotes a suitable loss function, W t and Z t denote the training datasets indexed by t ∈ {1, 2, … , n t } with n t the number of training datasets, where Z t corresponds to the aliased training images as input of the network and W t can correspond to both fully sampled training images or fully sampled training k-space data, depending on the application. 215 As an example, a model-based loss function based on the sum of squared differences where training is performed using fully sampled k-space data 215 can be incorporated as follows: The neural network can be included in an iterative reconstruction. 188,[216][217][218] In this case, the optimization problem described by Equation (7) is modified such that the regularizer term includes a generative neural network as follows: The neural network can be trained on previously acquired datasets, 219,220 on synthetic datasets, 221 on a combination of synthetic and real datasets, 214 or on a dictionary generated by a biophysical signal model. 216,217 The training process can be performed on image patches instead of full field-of-view images, 222 which allows the design of a good-performing neural network even with a small amount of training data. Also, self-supervised learning-based techniques that do not require fully sampled data have been proposed. 213 The input data are often comprised of magnitude images. 181,215,222,223 However, in applications where the phase information is important for the subsequent parameter estimation step, complex neural networks have been developed that exploit the correlation between the real and imaginary parts of a complex image by processing the real and imaginary images with complex convolutions inspired by the multiplication of complex numbers. 219 Various architectures have been investigated for the neural network, including U-NET, 215,224 fully connected convolutional neural network, 223,225 multiscale residual network, 222 deep cascade of residual dense network, 221 and deep complex residual network. 219

Parameter estimation
The second step of the indirect reconstruction is parameter estimation. Papers that target to improve this step are categorized in this subsection. Apart from the conventional model-based technique introduced in Section 3, the parameter maps X can be estimated from the reconstructed contrast imagesÎ using dictionary matching and learning-based estimation.

Dictionary matching
In techniques in the dictionary matching sub-category, a model-based dictionary D ∈ C n c ×n d is generated that relates parameter values for specific scan settings to the signal evolution using the analytical model  . In a qMRI experiment with n q parameters to be estimated, the range of possible values for each parameter can be digitized with a particular number of steps and step size. If the parameters indexed by q ∈ {1, … , n q } have ranges that are divided into n d 1 , … , n d q steps, then the number of steps for all the parameters included in the dictionary is n d = ∏ q n d q . Note that n d is also the number of expected signal evolutions defined in a dictionary. D d is the signal evolution for one set of parameters and is referred to as a dictionary atom. The size of the dictionary scales exponentially with the number of qMRI parameters n q , which affects the computational cost of generation and usage of the dictionary. Hence, dictionaries are often generated with a small number of parameters, with a typical maximum of four. 226 Matching the qMRI parameters to the measured signal (pattern matching) can be performed by solving the optimization problem described by: where  is a function that computes the difference between the reconstructed signal and the simulated signal in the dictionary in each voxel. The set of parameters corresponding with thedth column of the dictionary are the estimated qMRI parameters in the jth voxel. For dictionary matching,  only needs to be evaluated once to create D for given acquisition settings. Various advanced simulation techniques have been used to implement  , such as Bloch simulations, 40 65,102,103,112,142,167,184,196,197,205,208,211,[237][238][239] and methods using advanced biophysical models of the tissue. 138,192,207,240 Auxiliary parameter estimation. The estimation of qMRI parameters can be improved if the effects of scanner-related parameters, such as B 0 and B 1 field inhomogeneities, are included in the dictionary generation. Such parameters are referred to as "auxiliary parameters" in this work. To include auxiliary parameters in the dictionary generation, the dictionary needs to be extended. Such extensions may include slice profile effects, 102,112,119,120,197,230 B 1 inhomogeneity, 36,66,165,197,[228][229][230][241][242][243] receiver phase, 230 B 0 inhomogeneity, 40,233,244 relaxation effects during the inversion and preparation pulses, 120 and partial volume effects. 227 Dictionary compression. The precomputed dictionary for dictionary matching often has a considerable size and can be compressed to facilitate its implementation and reduce the required number of comparisons. Such compression has been performed using different techniques and in different dimensions. For instance, singular value decomposition, 88,89,137,169,211,218,228,229,245 randomized singular value decomposition, 231,244 and group matching 246 are used to compress the dictionary in the contrast dimension. Similarly, the step size for generating different parameters in the dictionary can be increased, resulting in a coarser dictionary with lower computational costs. This approach is often combined with interpolation techniques, such as quadratic interpolation 244 or B-spline interpolation 230 of dictionary entries to reduce possible artifacts. Specifically, an automatic technique is proposed to determine the dictionary resolution corresponding to a specified interpolation error. 230 Pattern matching. Some studies focus on improving the efficiency of the dictionary matching step. For instance, the iterative brute-force searches are replaced with fast approximate nearest neighbor searches based on cover tree structures. 131 The gradient descent search algorithm, 245 and a learned Mahalanobis distance 233 are used for finding the closest match with dictionary entries. The robustness of results with respect to flow artifacts and intravoxel dephasing is improved by applying matching in multiple steps. 228 Combination with other techniques. Further improvements to techniques in the dictionary matching subcategory include modifying the dictionary to make it applicable in combination with other qMRI reconstruction and acquisition techniques. For instance, a combination of view-sharing reconstruction and dictionary matching requires D to be modified to account for the mixed contrast in the data. 205 Dictionary matching is also used in conjunction with SMS acquisition to cover a given volume faster.
Reconstruction from SMS data requires the signal of each slice to be unfolded. Next, the unfolded signal of each slice is matched to a precomputed dictionary. 40 However, if the acquisition settings vary across slices, the dictionary needs to be modified accordingly, and a separate dictionary for each slice needs to be generated. 232 Furthermore, dictionary matching can be followed by application of a convolutional neural network to reduce motion-induced artifacts. 134

Learning-based estimation
Learning-based estimation techniques use neural networks that learn an approximation of the inverse of the qMRI model function  described in Section 3.2 to reconstruct high-quality parameter maps from (often) aliased and noisy contrast images that are reconstructed in the image reconstruction step (cf. Section 3.1). A neural network (Ž| ) is used where the estimated parameter maps are the output of the neural network:X = (Ž| ). 38,42,85,132,178,209,224,[247][248][249][250][251][252][253][254][255][256][257][258][259] The network training procedure is similar to Equation (20) but with Z t as the set of (aliased) training images and W t denoting the artifact-free training parameter maps.
Various loss functions have been proposed to improve the performance of the neural network. For instance, the loss function is defined as the squared error of each estimate normalized by the corresponding Cramér-Rao bound (CRB) before averaging all the estimated parameters and all the samples in the training data. 170 Such a loss function that includes the CRB allows the network to account for variations in the noise propagation among the parameters. 170 Similarly, a total variation regularizer along with the mean squared error is used as a loss function, 258 which results in smooth reconstructed parameter maps.
Complex-valued contrast images have been used for training and testing the neural network, showing a high fidelity of parameter quantification. 133,260 Learning-based parameter estimation can be improved by joint estimation of the B 1 field map along with the relaxometry parameters. 223,261 Learning-based techniques relying on random-forest regression 45 and regression with kernels 262 have been proposed to estimate artifact-free quantitative parameter maps from the reconstructed contrast images which can have undersampling artifacts. Furthermore, deep learning frameworks have been applied to accelerate parameter estimation, particularly when many parameters need to be estimated. 263 It has been proposed to use attention-based neural networks to address the need for explainable architectures in neural networks used for parameter matching. 264

Direct reconstruction
Techniques in this category skip the image reconstruction step and directly estimate X from undersampled Y, usually through an iterative process. They can be categorized in two subcategories.

Direct model-based reconstruction
The techniques in the direct model-based reconstruction subcategory exploit the redundancies in the data by incorporating a physical signal model  in the reconstruction. [265][266][267][268][269] The basic formulation for direct model-based reconstruction is given by: Adjustment of the physical signal model. The signal model  can be adjusted to improve parameter estimation. The related studies are often focused on: • Developing signal models that more accurately describe the tissue characteristics and acquisition settings, such as the inversion recovery Look-Locker model for T 1 mapping, 270 the variable flip angle model for T 1 mapping, 271 the chemical shift model for water-fat separation imaging, 272 the echo-modulation curve model for T 2 mapping, 267 the Kalman filter model for T 2 mapping, 273 the generating function formalism for T 2 mapping, 274 and models based on Bloch equations for direct estimation of parameters from transient response signals. [275][276][277][278] • Developing novel models that parameterize the model jointly for auxiliary parameters. For instance: the joint estimation of off-resonance frequency in T * 2 mapping, 279 the joint estimation of T * 2 and field maps in water and fat parameter estimation, 272,280,281 the joint estimation of flip angles and proton density maps in T 2 mapping, 282 the joint estimation of steady-state signal, equilibrium signal, and effective relaxation rate in T 1 mapping, [283][284][285] the joint estimation of tracer kinetic model parameters in T 1 mapping, 286 and the joint estimation of phase parameters in multishot diffusion imaging. 287,288 Regularized direct reconstruction. To improve the performance of direct reconstruction, prior information about the qMRI parameters to be estimated can be incorporated by adding regularization terms, such as Tikhonov regularizer, 274 l 1 norm regularizer, 289 and weighted l 1 ball regularizer 268,290 to the cost function described by Equation (24). These priors can promote the sparsity of the parameter maps in different transform domains, such as total variation transform, 269,286,[291][292][293][294] the wavelet transforms, 281,[283][284][285][295][296][297][298] the total generalized variation transform, 299 the fractional variation constraint, 282 the surfacelet transform, 300,301 patch-based difference operator, 302 and the finite difference transform. 303 In these cases, the optimization problem (24) is modified to include regularizer terms yieldinĝ Furthermore, it has been proposed to jointly estimate the coil sensitivity maps along with the qMRI parameters in direct reconstruction, where a Sobolev norm is applied to these maps to enforce smoothness. 284,285,296 Additionally, to enforce phase smoothness in the joint reconstruction of phase and qMRI parameters, application of an l 2 norm regularization on the spatial gradient of the phase maps has been proposed. 294 Joint image reconstruction and parameter estimation. In some applications, parameter maps are estimated jointly with the images: 304,305 X,Î = arg min Estimating both contrast images and parameter maps allows the application of prior information in both the parameter and image domains. To promote local coherence of contrast images, it has been proposed to apply a low-rank prior on contrast images using a nuclear norm regularizer term in direct reconstruction. 286 Combination with other techniques. Techniques in the direct model-based reconstruction subcategory have been combined with other qMRI reconstruction techniques. High acceleration rates have been achieved by a combination of direct reconstruction of T 2 maps with PI techniques where the first echo in a multi-echo gradient echo scans is used for GRAPPA calibration. 306 The combination of direct reconstruction and PI allows for exploiting the prior knowledge from both local k-space dependencies and the physical signal model, which has been investigated in various organs in human body. [306][307][308] Furthermore, direct reconstruction has been combined with SMS acquisitions, by incorporating the SMS encoding matrix into the encoding operator . 285 Optimization algorithms. Minimizing the cost function of a direct model-based reconstruction problem is often challenging as it involves a nonlinear optimization problem. To solve this problem, exponential models can be approximated with their Taylor series expansion 281,302 or linearized in a small neighborhood using the Gauss-Newton approach. 299 Additionally, the nonlinear optimization problem can be treated as a tomography problem and solved using large-scale nonlinear optimization routines. 278,309 The computational time of techniques in direct reconstruction can be decreased using block-wise decomposition of the cost function in combination with a specific uniform undersampling pattern. 310 In case multiparametric models are used, advanced solvers such as Stochastic Gradient Langevin dynamics, 286 regularized trust region continuation techniques, 279 or alternating minimization 288,295 have been proposed to optimize the nonlinear cost-function jointly for different parameters. Furthermore, it has been proposed to automatically scale the unknown parameters in the joint optimization, which numerically balances the partial derivatives of the multiparametric model to improve the conditioning and hence ease the optimization. 311

Direct learning-based reconstruction
Techniques in the direct learning-based reconstruction subcategory aim to replace the direct model-based reconstruction with a neural network to estimate qMRI parameter maps directly from the acquired k-space data. The neural network  has undersampled k-space data as input, that is, Ž = Y and estimated qMRI parameter maps as output, that is,X = (Ž| ). The training process for optimizing network parameters is similar to Equation (20), with W t and Z t denoting the reference parameter maps and undersampled k-space data of a training dataset, respectively.
The computational costs of training a neural network that performs direct estimation of parameter maps from k-space data are high. Therefore, alternatively, neural networks have been proposed to improve parts of the techniques in the direct model-based reconstruction approach. [312][313][314][315][316] For instance, deep learning networks have been designed to regularize the estimation of parameter maps in Equation (24). [312][313][314] Similarly, neural networks have been proposed to estimate the initial parameter maps of the direct model-based reconstruction from undersampled k-space data. 314 Further, it has been proposed to predict the derivatives of a direct forward model using a pretrained neural network, decreasing the computation time for direct model-based reconstruction. 315 Finally, a self-supervised learning-based approach has been proposed to perform real-time parameter mapping for transient imaging. 316 In this approach, the Fourier Transform of the parameter maps is computed from the undersampled k-space data, after which parameter maps are obtained by Fourier Transform.

RESULTS: GRAPHICAL REPRESENTATION
The current section presents figures that visualize the distribution of the 292 reviewed papers over different categories. Figure 3A shows that 240 papers were categorized as indirect reconstruction, whereas the remaining 52 papers were categorized as direct reconstruction. Figure 3B shows a further subcategorization of the indirect reconstruction category, classifying the papers based on their image reconstruction step (top) and parameter estimation step (bottom). Note that the subcategories FTc

F I G U R E 3
Distribution of articles included in the review: (A) direct versus indirect categories; (B) Image reconstruction (top) and parameter estimation (bottom) within the indirect category.

F I G U R E 4
The distribution of articles included in this review over the different classes of target quantitative magnetic resonance imaging parameters. The letter J in front of the parameter indicates joint estimation of multiple quantitative parameters.
(top) and model fitting (bottom) correspond with the conventional image reconstruction and parameter estimation step of a classical two-step approach, respectively. Obviously, no paper selected in this review classifies for both these subcategories, since the inclusion criteria require a modification of at least one step of the conventional two-step approach. Figure 4 shows the distribution of the articles over the following classes of target qMRI parameters: 1. Relaxometry which includes T 1 , T 1 , T 1H 2 O , T 2 , T * 2 , proton density, spin density, magnetization transfer, myelin water fraction, quantitative susceptibility mapping, fat fraction, and water fraction. 2. Diffusion which includes apparent diffusion coefficients, mean diffusivity, fractional anisotropy, and parameters of more advanced diffusion models that are mentioned in the table available online. 3. Temperature based on proton resonance frequency shift.
4. Perfusion which includes parameters of arterial spin labeling and tracer-kinetic models. 5. Flow based on velocity mapping. 6. Others which consists of Hyperpolarized 13C, tissue sodium concentration, and Paramagnetic fluorine-19 parameters. Figure 4 shows that about 70% of the papers reviewed in this work focus on relaxometry. Figure 5 shows the number of papers reviewed in this work that consider a specific category of target parameters while presenting a novelty that can be categorized to one of the image reconstruction subcategories ( Figure 5A), parameter estimation subcategories ( Figure 5B) or direct techniques subcategories ( Figure 5B) that were described in Section 3 and 4 and summarized in Figure 2.
Similarly, Figure 6 shows the number of papers that target a specific application domain while presenting a novelty that can be categorized to one of the image reconstruction subcategories ( Figure 6A), parameter estimation subcategories ( Figure 6B), or direct techniques subcategories ( Figure 6B). It can be observed from Figure 6 that the main application domain is brain MRI, while lung MRI (not specified, but included in the category "others") is the least covered. Figure 7 demonstrates the distribution of articles over the publication year. It can be seen that 167 out of the 292 reviewed papers were published from 2019 to 2022. After 2012, the variety of targeted parameters ( Figure 7A) and the application domains ( Figure 7B) increased. Furthermore, it can be observed in Figure 7C,D that the learning-based techniques started to be used in parameter estimation in 2017 and in image reconstruction in 2019. Figure7E shows the number of articles that included a public link to the acquired datasets or the implemented code, in different years. The individual links are available in Table 3 and are summarized in this GitBook page.

F I G U R E 5
The number of articles that consider a specific class of target parameters, while presenting a novelty that can be categorized to one of the image reconstruction subcategories (A), parameter estimation subcategories (B), or direct reconstruction subcategories (B). The letter J in front of the parameter indicates joint estimation of multiple quantitative parameters.

F I G U R E 6
The number of articles that consider a specific application domain, while presenting a novelty that can be categorized to one of the image reconstruction subcategories (A), parameter estimation subcategories (B), or direct reconstruction subcategories (B).

DISCUSSION
This work describes and categorizes the different reconstruction techniques used for undersampled qMRI acquisitions in a unified mathematical framework, focusing on their methodological contribution while ignoring case-specific variables. Accelerated qMRI methods have enabled for example T 1 and T 2 mapping of the whole brain with 1 mm 3 resolution in just 3 min 126,152 and cardiac T 1 and T 2 mapping with a single breath-hold scan in about 16 s. 141 Accelerated qMRI is an active, fast-growing research field, a trend supported by the observation that around half of the papers included in this work were published in the last 3 years. The growing number of publications and the emergence of new methodologies, such as learning-based techniques, indicate the interest of the research community in further acceleration of qMRI as well as in more accurate and trustworthy techniques. Figure 7E, shows an increasing trend in the number of articles with publicly available code or data, suggesting improvements in transparent developments and collaborations.
The distribution of articles over the years suggests a gradually increasing application of all techniques. However, a jump can be observed in 2016 for FTc and multicontrast PI subcategories in image reconstruction and for dictionary matching subcategory in parameter estimation. This jump can be related to the emergence of transient response multi-parametric imaging techniques 236 that often reconstruct the contrast images using FTc and PI techniques and remove the undersampling artifacts in the parameter estimation step using dictionar y matching.
The distribution over target parameters that can be observed in Figure 7A shows an inclination toward investigating relaxometry. The reasons contributing to this could be the large variety and accessibility of relaxometry imaging sequences, and the wider availability of relaxation phantoms compared to other parameters. Other modalities, such as diffusion MRI, result in contrast images with lower signal-to-noise ratio values compared to relaxometry imaging, which creates a need for more data samples. 317 In perfusion imaging, in addition to the problem of low signal-to-noise ratio, there is a lack of a gold standard for validating the results, as designing a perfusion phantom is challenging. 318 These reasons can contribute to making relaxometry a better candidate for undersampled acquisitions.
The brain is the most investigated application domain for qMRI techniques. With the many properties that MRI can measure, it is uniquely positioned to study the complex composition and abnormalities in the brain. With the relative absence of motion and other disturbing factors, combined with strong interest from the neurological and neuroscience communities, it is logical that many techniques are first proposed for brain imaging. Next to the brain, cardiac applications are highly investigated with accelerated qMRI techniques. qMRI enhances the knowledge of cardiac function, volume, and tissue characterization. It allows for the spatial visualization of changes in the myocardium based on changes in myocardial relaxation times and flows, enabling the evaluation of diffuse changes within the myocardium. In this way, cardiac qMRI facilitates exploring the link between biology and the clinical manifestation of cardiac diseases and is an important imaging modality. However, cardiac imaging is more challenging than brain imaging due to substantial nonrigid motion during the scan, which is unavoidable during the long scan time of qMRI. Therefore, qMRI techniques facilitated by simultaneous motion correction are of great interest. 153,212 Lung imaging is only one of the least investigated application domains in the included studies. The reasons contributing to this are the lung's sparse soft tissue structures, low proton density, and larger susceptibility variations due to multiple interfaces between air and soft tissue, which hinder MRI signal generation. The lack of multiple contrasts in lung tissue limits the clinical application of qMRI for this organ. 319 Most of the techniques discussed in this work are applied to data acquired with MR sequence settings, such as echo time and flip angles, optimized for a fully sampled scan. Additionally optimizing the undersampling pattern along with MR sequence settings can further improve the results. 320 However, given the many possible combinations of acquisition settings and undersampling patterns, empirical optimization of in vivo precision is impractical. Hence, efforts for in-silico evaluation, such as predicting time efficiency, 12,266,[320][321][322][323][324][325] or accuracy, 118 are relevant. Recently, automated learning-based methodologies have been proposed to select an optimal sampling strategy independent of the model. 326 Similarly, owing to the development of modern data science tools such as auto differentiation ability, it is possible to develop neural networks to estimate parameter maps by minimizing a nonlinear loss function based on the Bloch equations. 327 The present study has some limitations. Although we report on the theoretical principle behind each technique, we do not report on implementation details as this was considered beyond the scope of this work. The implementation not only affects the computational cost and run time of each technique, but the choice of solver also affects the finding of the optimal solution. We neither aim to rank techniques because comparing performance in terms of reported acceleration factor, scan time, resolution, and signal-to-noise ratio is complicated by the significant differences in anatomical locations, hardware, and parameters with which the techniques have been presented. It should be noted that recent techniques for qMRI have been proposed that directly estimate qMRI maps from a single contrast image, using data-driven strategies. [328][329][330] These techniques are not included in this study.
qMRI is a promising field, but scan time reduction is only one of the main obstacles to its clinical application. In addition to the long scan time, accurate estimation of parameters is a fundamental challenge in qMRI, 331 which becomes even more prominent in accelerated qMRI. The reconstruction techniques often introduce some regularization and impose assumptions on parameter values that may be invalid and bias the estimated parameters. This creates a trade-off between acceleration and accuracy, which needs to be studied prior to introducing a technique. Moreover, it is crucial to account for variabilities introduced by hardware and acquisition settings. 332 To implement qMRI in clinical routines, further research on clinical interpretation of the biomarkers extracted from qMRI maps is necessary, since these biomarkers are conventionally defined on the weighted images. Synthetic weighted images produced from qMRI could bridge the direct use of quantitative maps and the current scheme where weighted images are used. 133 Clinical routines are often a combination of various sequences for different weighted images. A single qMRI scan from which synthetic weighted images are produced could reduce the repetitive acquisition of information shared among the weighted images, which can save scan time.
The unified mathematical framework presented in this work facilitates comparing accelerated qMRI techniques on theoretical grounds, such as the type of prior assumptions made. It highlights the essential aspects critical to the state-of-the-art qMRI reconstruction techniques, which may guide future studies. Furthermore, this review study provides an overview of the distribution of reviewed reconstruction techniques over the application domains, parameters of interest, and the years of publication, which facilitates exploring existing trends and gaps in current studies. We hope that this information helps researchers to propose new combinations of techniques or find new applications for promising techniques.

APPENDIX T A B L E A1
List of the symbols and mathematical notations.
The low-rank component of I when I is separated into its sparse and low-rank component

T A B L E A2
The search query used for each database.

T A B L E A3
List of the reviewed papers labeled with first authors' last names (Author) along with the year of the publication (year), reference number in this paper (Ref), focused qMRI parameter (Parameter), qMRI reconstruction categories that are proposed considering indirect and direct techniques (category) and the subcategories used for image reconstruction (Recon) and parameter estimation (Param) in indirect categories, application domain with respect to the type of the study (App) and the organ that is studied (Organ), and finally the accessibility of the study regarding the links to data or code that are provided (Data/Code).

Author
Year