Structure Prior Effects in Bayesian Approaches of Quantitative Susceptibility Mapping

Quantitative susceptibility mapping (QSM) has shown its potential for anatomical and functional MRI, as it can quantify, for in vivo tissues, magnetic biomarkers and contrast agents which have differential susceptibilities to the surroundings substances. For reconstructing the QSM with a single orientation, various methods have been proposed to identify a unique solution for the susceptibility map. Bayesian QSM approach is the major type which uses various regularization terms, such as a piece-wise constant, a smooth, a sparse, or a morphological prior. Six QSM algorithms with or without structure prior are systematically discussed to address the structure prior effects. The methods are evaluated using simulations, phantom experiments with the given susceptibility, and human brain data. The accuracy and image quality of QSM were increased when using structure prior in the simulation and phantom compared to same regularization term without it, respectively. The image quality of QSM method using the structure prior is better comparing, respectively, to the method without it by either sharpening the image or reducing streaking artifacts in vivo. The structure priors improve the performance of the various QSMs using regularized minimization including L1, L2, and TV norm.


Introduction
Quantitative susceptibility mapping (QSM) generates tissue magnetic susceptibility property image from susceptibilitysensitive MRI data [1][2][3][4]. QSM can reduce blooming artifacts in susceptibility weighted imaging [5], is clinically useful for quantifying magnetic biomarkers that have susceptibilities different from their surroundings [6][7][8][9], and promises to probe oxygen metabolism [10][11][12] and inflammation [13]. The basic QSM physics includes estimating the local magnetic field (r) (measured main magnetic field B 0 ) from the MRI signal phase [14,15] and modeling the field as from dipole sources in tissue: = d ⊗ + n in image space (referred to as r-space) or Δ = DX + N in k-space, where d is a unit dipole in r-space (referred as dipole kernel), is the tissue susceptibility distribution in r-space, n is the observed noise in r-space, and Δ , X, D, and N are corresponding to , d, , and n in k-space, respectively [16][17][18]. The fundamental QSM algorithm is to solve the inverse problem from the field to the susceptibility source . Because the dipole kernel has a pair of opposing zero cone surfaces at the magic angle (54.7 ∘ ) with respect to the B 0 direction, the inversion from field to susceptibility is fundamentally ill-posed, and QSM requires additional information to select a unique susceptibility map from many possibilities [19]. For clinical data acquired using a single orientation, various regularization methods have been proposed for QSM [17]. The Bayesian formulation provides the main QSM approach to identify a susceptibility distribution of minimal streaking artifacts [20][21][22].
Bayesian QSM requires constructing a data fidelity term according to data noise property and a prior information term. Previously, we examine the importance of noise whitening in constructing the data fidelity term [17]. Here, we examine the importance of including structure information in the prior term. Specifically, we analyze the effects of the choice of prior in QSM based on three published methods in literature and several new algorithms with different regularization terms. All QSM methods are evaluated using simulations and phantom experiments where the susceptibility is known. Further, clinical applicability of QSM was evaluated with human brain data to investigate image quality.

Material and Methods
In Bayesian approach, the regularization prior is expressed as a cost function R that favors a solution of the desired property, and the degree to which it is favored is typically characterized by a regularization parameter . The maximum a posterior solution [17,23] is where E constitutes the data fidelity term. Since Gaussian noise in the complex MR signal domain should be accounted for in the data fidelity term and with proper noise weighting, noise effects in QSM can be reduced using Bayesian methods [17]. In the following section, we use E = ‖wz‖ 2 2 , with z = d ⊗ − and w as noise weighting.

Priors Used in Various
Algorithms. The following examples of regularization terms have been explored [17]: (iii) The total variation norm (TV) (iv) A wavelet domain such as a Daubechies wavelet (Φ) L1 norm (v) A combination of two sparsity terms such as total variation and wavelet [24] → Φ 1 + TV ( ) (vi) An L1 norm with structural consistency term m (MEDI) [25][26][27][28] = m 1 We evaluated the quantification accuracy and image quality of various QSM algorithms using numerical and experimental phantom data and image quality in clinical applications with in vivo patient data.

QSM Algorithms in This
Paper. Six representative algorithms are summarized in Table 1 and in pseudocodes in the Appendix for reference: MGL2 (GL2 with structural consistency prior m), MTV (TV with m), and MEDI (GL1 with m). The equations of the GL2 and MGL2 methods were solved using the conjugate gradient method (CG). The MEDI, TV, MTV, and GL1 methods were solved using a lagged diffusivity fixed point method (LDFP). In methods involving a structure prior, m was estimated by setting the gradients of magnitude image greater than a certain threshold to 0 and to 1 otherwise [28]. The threshold was adjusted iteratively such that approximately 30% of the voxels within the interested brain region were 0, that is, considered to have gradients. This threshold was determined from a previous theoretical study where the threshold was varied and a global minimum was empirically found at 30% [19].
We applied a one-dimensional temporal phase unwrapping and a linear least squares fitting to estimate the rate of phase evolution to get field map. The noise weighting w was estimated as the SNR along with the field map estimation [17]. The PDF method [14] was used for removing the background field. The local field outside the brain parenchyma was cropped by a mask, which was manually segmented in the numerical and phantom experiments and was obtained using a Robust Brain Extraction (ROBEX) tool for in vivo brain data [29]. Voxels in the background region or within 3 mm to the background region were set to zero.

Numerical Simulation.
A 256 × 256 × 128 Zubal-type [30] numerical susceptibility phantom ( Figure 1) was built containing a spherical lesion of 5-pixel radius to simulate a low SNR region within the T2 * magnitude image of the brain parenchyma. The simulated susceptibility values were −0.05, 0.07, 0.09, 0.09, 0.19, 0.30, 0.90, and 0.00 ppm for white matter, thalamus, caudate nucleus, putamen, globus pallidus, veins, lesion, and other parenchyma, respectively. Complex Gaussian zero-mean noise with standard deviation ranging from 0.01 to 0.05 was added to the simulated complex MRI signal.

Phantom Experiments.
A Gd phantom was constructed consisting of six spherical balloons filled with solutions of various concentrations of Gd-DTPA (Magnevist; Berlex Laboratories, Wayne, NJ) and immersed in water within a cylindrical container with dimensions of 12.5 cm diameter and 30 cm height. The Gd concentrations ranged from 0.5% to 3.0% (using a 0.5% increment) with susceptibilities ranging from 0.81 to 4.89 ppm, corresponding to the Gd concentration in aortic arch measured in the first pass of a dynamic contrast enhanced MRI [18]. This phantom was scanned on a 1.5T MRI scanner (HDx, GE Healthcare, Waukesha, WI) with a body coil. A T2 * weighted multiecho gradient echo sequence was performed with the following parameters: FA = 30 ∘ ; TR = 30 ms; number of TEs = 3; first TE = 3.05 ms; uniform TE spacing (ΔTE) = 1.0 ms; BW = ±31.2 kHz; FOV = 12.8 cm; and acquired resolution = 2 × 2 × 2 mm 3 .

Implementation Details for
Algorithms. The regularization parameter ( ) was searched from 10 −5 to 10 1 (13 logarithmically spaced steps) for all QSM methods. The best parameter of every method for noisy numerical simulation and phantom was chosen according to the least error with respect to the true susceptibility [17,31] and to the prepared susceptibility in each balloon. Because the true susceptibility is not available in the real human brain, the best parameter of each method is chosen according to the balance of artifacts and contrast among brain components in one representative case by the 13-year experienced neuroradiologist (W. C.) with inspecting all the varying parameter's results reconstructed by the method. And all real cases use the chosen parameters.

Data Analysis.
The normalized root mean square error (NRMSE) (normalized by root mean square of true susceptibility) of the whole volume and the linear regression of the measured versus known susceptibility was calculated for every method. This NRMSE is used to assess the accuracy of the numerical phantom reconstruction. In the Gd phantom, a linear regression between the measured and prepared susceptibilities of the various balloons was performed for each of the methods.
To assess difference between two QSM algorithms in numerical and experimental phantoms, the values of susceptibilities are measured and compared to their known values. Their differences between two methods were assessed using paired -test based on comparing linear regressions, and the method with the smaller error was reported as improvement when < 0.05.
To assess the quality of patient images reconstructed by different QSM methods, three neuroradiologists (W. C., C. P., and K. M.) reviewed images simultaneously in a random order blinded to the reconstruction methods. Overall image quality was scored with a 5-point scale based on radiological impression of smoothness and artifacts, where 5 is the highest quality and 1 is the lowest quality [17]. The radiologist (W. C.) assessed the image quality again 6 months later to assess intraobserver variability. Interobserver and intraobserver variabilities of image quality scores were assessed by using the intraclass correlation coefficient [32]. The following criteria for clinically relevant agreement were used to assess the calculated intraclass correlation coefficient: less than 0.40 was considered poor; 0.40-0.59, fair; 0.60-0.74, good; and greater than 0.74, excellent [33]. The significance of image score differences between reconstructed susceptibility maps assessed by the Wilcoxon rank sum test. Statistically significant with the higher image score was reported as an improvement when < 0.05.

Numerical Simulation.
The structural prior m markedly reduced streaking artifacts originating from the simulated lesion for GL2 ( = 9.22×10 −7 ), TV ( = 0.02), and GL1 ( = 0.04) as seen in Figure 1(a), which was also reflected in the reduced NRMSE (Figure 1(b)) and in the increased regression slopes more closely approaching 1 as seen in Figure 1(c). The QSM reconstruction accuracy was improved from MGL2 to MTV ( = 3.49 × 10 −3 ) and from MTV to MEDI ( = 0.03).   Table 1 for every method, respectively. The spatial prior either sharpened the image or reduced streaking artifacts (Figure 2 black arrows). The overall image quality scores were higher when using structural consistency (m) prior compared to the same method without a structural consistency prior observed in the reconstructions of group 1. Compared to the method without m, the artifact from vessels (black arrow in Figure 2) is reduced and the appearance of small veins (dot boxes in Figure 2) is improved in the method with m, respectively. The difference is statistically significant ( = 6.37 × 10 −3 , 9.38 × 10 −6 , 2.52 × 10 −3 for GL2, TV, and GL1, resp.) between methods with and without m. The QSM reconstruction image quality was improved (Table 1 score  The spatial prior either sharpened the image or reduced streaking artifacts (Figure 3). The differences in overall image quality (Table 1 score for group 2) were statistically significant ( = 4.66×10 −2 , 1.89×10 −5 , 1.80×10 −4 for GL2, TV, and GL1, resp.) between reconstructions with and without structural consistency prior method when they are performed in group 2. The QSM reconstruction image quality improved from MGL2 to MTV ( = 9.47 × 10 −7 ), MGL2 to MEDI ( = 3.19 × 10 −7 ), and MTV to MEDI ( = 7.21 × 10 −4 ). The inter-and intraobserver results are shown in Table 2. These agreements ranged between good and excellent referred to the clinically relevant agreement.

Discussion
Our results of the various QSM methods indicate that QSM quality improves using a physical prior specific to the imaging situation. Investigations with a numerical Zubal lesion phantom, a Gd phantom, and 36 consecutive patients consistently . Similar to healthy subject images, the structure consistency m reduces overall streaking artifacts comparing QSM methods with and without it. The streaking artifacts were seen originating from vessels (black arrows) and the hemorrhage (hollow arrow heads) in GL2 and MGL in (c). This artifact was reduced to some extent in TV and MTV and even further in GL1 and MEDI. The averaging image quality scores of this patient were 1, 2, 2, 3.3, 3.3, and 4 for GL2, MGL2, TV, MTV, GL1, and MEDI, respectively. 7 corroborated this observation. These experimental results are concordant with the theoretical error analysis in QSM that the error in the reconstructed susceptibility comes from noise in the data and error in the prior [1,19]. The ill-posedness of the dipole inversion sets up a very challenging problem for QSM. The iterative solvers employed in more accurate QSM algorithms make it difficult to understand the contributions from various terms and parameters. The systematic comparison of various methods presented in this study tends to suggest the structure prior is useful to QSM methods using the regularized minimization. This observation is established in the comparisons in Table 1. The structural matching between the magnitude image and susceptibility map by matching their gradients tends to improve the QSM image quality and accuracy as seen in the paired comparison with and without using the magnitude gradient in Figures 1-3. It should be pointed that the new simulation with the different brain components' susceptibility has been run and the real cases were also particularly acquired with slightly different parameters in different time compared to our previous simulation and real cases in [17].
The value of the phantom experiment may be limited for assessing the performance of QSM algorithms, because the phantom was made of sphere of uniform [Gd] with perfectly identifiable edges. According to the above observation on QSM algorithms, structural information as defined by tissue contrast is a very important determinant for QSM performance. The tissue contrast of a human brain is much more subtle and complex than that of the phantom made of spheres. So what is learned from the phantom experiment is only the accuracies that various QSM methods can achieve under an ideal contrast situation. It should be pointed out that the susceptibility range of phantom experiment was much higher than that of in vivo human brain and the measurements were consistent with our simulation at similar susceptibility values. QSM algorithm comparison behaves similarly at a wide range of susceptibility values, which may be explainable with dimension analysis. When susceptibility is scaled and the regularization parameter is also scaled correspondingly (unchanged for L2 and scaled by the same factor for L1 and TV), the performances of the QSM algorithms remain the same.
It is also observed that GL1 is better than GL2 for structural matching (Figures 1-3), which is consistent with previous observations that the sparsity in the prior term can be achieved more effectively with the L1 norm than the L2 norm [28,34]. A recent publication also reached a similar conclusion that the L1 norm promoting the sparsity of spatial gradient is a better fit for QSM problem than the L2 norm [23,35]. GL1 and TV based methods both can be categorized into L1 norm QSM algorithm. It was also noted that GL1 is better than TV methods (Figures 1-3). This may be due to the fact that the cost of a particular edge in the objective function is the same regardless of its orientation in the TV norm, whereas the cost in the GL1 norm is dependent on the edge orientation. When there is no strong streaking artifacts as in Figure 2, the difference between GL1 and TV or between MEDI and MTV is small. When there is a large amount of streaking artifacts as in the hemorrhage case in Figure 3, the edge orientation sensitivity in GL1 and MEDI becomes obviously advantageous over the TV and MTV in suppressing the streaking artifacts.
The results by different solvers used in this paper, that is, LDFP and CG, show the consistency improvement of using structure prior. There are other solvers available including interior point methods or split Bregman method [36]. The reconstructed QSM image may be affected by the solver used in a given algorithm. However, in order to obtain reasonable answers to the concerning question on solver effects, the error propagation and accumulation in each solver may require detailed quantitative and analytic evaluation; other wellknown solvers may need to be implemented for all major QSM algorithms and included for comparison. However, this important work is beyond the scope of this paper.
Some folds in cortex in Figures 2 and 3 seem to be smoothed in MEDI; this may be interpreted as that the gradient consistency between magnitude image and susceptibility map used in current MEDI is still imperfect. Indeed, recent studies have suggested that the structure prior in current MEDI may be improved with more sophisticated identification and matching of structures for more accurate QSM, such as incorporating edges derived from the local field map [37,38]. It is also theoretically proven that a comprehensive detection of all the edges in the true susceptibility distribution will reduce the error in the reconstructed QSM [19]. However, the perfect prior as footstone of accurate regularization QSM method is still an ongoing research.

Conclusions
In summary, the structure prior with an effective match can improve the accuracy of QSM. Among compared methods, the MEDI method appears to be the most robust for quantitative susceptibility mapping. The more accurate priors or more physically meaningful priors should be studied in the future.