{'fraction_non_alphanumeric': 0.03523570232112946, 'fraction_numerical': 0.014596793491265853, 'mean_word_length': 4.647466216216216, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust Quantitative Susceptibility Mapping via Approximate Message Passing

Purpose: It is challenging to recover magnetic susceptibility in the presence of phase errors, which may be caused by noise or strong local-susceptibility shifts in cases of brain hemorrhage and calcification. We propose a Bayesian formulation for quantitative susceptibility mapping (QSM) where a customized Gaussian-mixture distribution is used to model the long-tailed noise distribution. Theory: Complex exponential functions of the phase are used as nonlinear measurements. Wavelet coefficients of the susceptibility map are modeled by the Laplace distribution. Measurement noise is modeled by a two-component Gaussian-mixture distribution, where the second component is reserved to model the noise outliers. The susceptibility map and distribution parameters are jointly recovered using approximate message passing (AMP). Methods: The proposed AMP with built-in parameter estimation (AMP-PE) is compared with the state-of-the-art nonlinear L1-QSM and MEDI approaches that adopt the L1-norm and L2-norm data-fidelity terms respectively. They are tested on the simulated and in vivo datasets. Results: On the simulated Sim2Snr1 dataset, AMP-PE achieved the lowest NRMSE and SSIM, MEDI achieved the lowest HFEN. On the in vivo datasets, AMP-PE is more robust and better at preserving structural details and removing streaking artifacts in the hemorrhage cases than L1-QSM and MEDI. Conclusion: By leveraging a customized Gaussian-mixture noise prior, AMP-PE achieves better performance in challenging cases of brain hemorrhage and calcification. It is equipped with built-in parameter estimation, which avoids subjective bias from the usual visual-tuning step of in vivo reconstruction.


INTRODUCTION
The quantitative susceptibility mapping (QSM) technique recovers magnetic susceptibility from magnetic resonance (MR) phase images 1,2,3,4,5,6 . It is widely used to study iron deposition in the brain 7,8,9 , or pathologies such as hemorrhage 10,11 and calcification 12,13 . Since the values of raw phase images are within (− , ], we need to unwrap the phases first to eliminate the discontinuity caused by the transition at (or − ) 14 . In order to define the brain as the region of interest (ROI), a binary brain mask can be generated using the Brain Extraction Tool 15 . Magnetic field variations in the ROI are extracted from unwrapped phase images, and contain 1) the background field induced by global geometry, air-tissue interfaces and field inhomogeneities; 2) the local field induced by the brain. The background field is then removed to produce the local field map 16,17 . Recovery of the susceptibility map from the local field map is known as the dipole inversion, it is an ill-posed inverse problem due to the zeros in the dipole kernel along the magic angle 1,18 . In this case we shall rely on prior information about the susceptibility map to improve the image quality. The widely adopted smoothness prior on implies that the image gradients are sparse and mostly close to zero. This gives rise to the total-variation (TV) minimization approach that regularizes a data-fidelity term using the 1 -norm of image gradients 19,20,21 . Alternatively, the image wavelet coefficients are also sparse in general, and we can choose the 1 -norm of wavelet coefficients as the regularization term instead 22,23,24 . The advantage of the wavelet basis over the total-variation basis is that the wavelet transform is invertible, which offers us more freedom in developing suitable reconstruction algorithms. Additionally, anatomical information obtained from the magnitude image could serve as the morphology prior on to ensure the reconstruction is consistent with anatomical structures 25 .
Apart from the ill-posedness of the dipole inversion, phase unwrapping errors also make the problem challenging. The erroneous phase jumps in the phase image lead to severe streaking artifacts in QSM when a linear least-squares datafidelity term is used. In order to handle these phase jumps, Liu et al. mapped the local fields to the complex domain via the complex exponential function, and proposed a nonlinear least-squares data-fidelity term that was much more robust 14 . Both the linear and nonlinear least-squares datafidelity terms imply that the corresponding noise is modeled as additive white Gaussian noise (AWGN). However, in cases that involve brain hemorrhage or calcification, the resulting noise can be better modeled by long-tailed distributions such as the Laplace distribution and the Gaussian-mixture distribution. To this end, Milovic et al. used the 1 -norm to construct an L1 data-fidelity term with improved performance 26 .
Regularization approaches like TV-minimization use a parameter to balance the trade-off between the data-fidelity term and the regularization term. The QSM reconstruction typically has a variety of noise profiles depending on the subject's condition and the chosen processing pipeline. When we do not know beforehand if the noise profiles between the training and test sets would match, using a pre-tuned (fixed) parameter might not be ideal. People thus have turned to the L-curve method to find the parameter adaptively for each dataset 27 . However, the L-curve method is inherently heuristic, and there is no guarantee that the selected parameter would be optimal. In practice, visual fine-tuning is typically used to find working parameters or double-check pre-tuned parameters for in vivo reconstructions 26 . However, it also introduces subjective bias to the results inevitably.
In this paper we propose a probabilistic Bayesian approach to recover the susceptibility map and parameters jointly. We use the Laplace distribution to encode the sparse prior on the wavelet coefficients of susceptibility map, and a customized Gaussian-mixture distribution to model the longtailed noise distribution. The maximum-a-posteriori (MAP) estimations of the wavelet coefficients and distribution parameters are computed using the approximate message passing (AMP) 28,29 . In order to handle the phase unwrapping errors, we adopt the nonlinear measurement model where the complex exponential functions of phases are used as measurements. We then extend the standard linear AMP so that it could be used to solve the nonlinear dipole inversion. Experiments show that the proposed AMP with built-in parameter estimation (AMP-PE) successfully recovers susceptibility maps from both the healthy and hemorrhage scans with more robust and better performances.

THEORY
Let denote the produced local field after phase unwrapping and background field removal, and denote the corresponding phase at an echo time . We have where is the gyromagnetic ratio, 0 is the initial phase offset that depends on the coil-sensitivity, 0 is the main magnetic field, is the Fourier transform matrix, is the dipole kernel in the frequency domain, is the magnetic susceptibility, is the noise, and is the (combined) resulting operator applied on . The initial phase 0 can be estimated from multi-echo phase images by solving a nonlinear least-squares fitting problem 14 . In particular, the ill-posed dipole kernel is given by where is a diagonal weighting matrix, is the noise. The complex exponential function exp( ) of the phase image could greatly reduce the discontinuity caused by erroneous phase jumps. There are different ways to design suitable weighting matrices, here we can simply use the magnitude image as . We then aim to recover from the nonlinear measurements exp( ).

Bayesian Formulation
The sparse prior on the wavelet coefficients of the susceptibility is used to improve the image quality: where is the invertible wavelet transform. As shown in Fig.  1 (a), we use the sparsity-promoting Laplace distribution to model the distribution of wavelet coefficients. The coefficients in are assumed to be independent and identically distributed (i.i.d.): where > 0 is the unknown distribution parameter. As shown in Fig. 1 (b), we propose the following Gaussian-mixture distribution to model the long-tailed noise distribution where is the -th mixture weight, is the -th variance, the mixture means are zeros, and  (⋅) is the Gaussian density function.
The number of Gaussian mixtures needs to be chosen to avoid overfitting during reconstruction. From the experiments later in Section 4, we observed that two Gaussian mixtures produce the best performance. In particular, the variance 2 of the second Gaussian component should be large enough to cover the domain of , it is used to model the noise outliers that give rise to a long-tailed distribution. In practice, we can initialize the second variance 2 with a larger value than the first variance 1 .
Under the Bayesian formulation, we can recover by computing its MAP estimation where contains the measurements, and ( | ) is the posterior distribution of that can be computed via approximate message passing (AMP).

Approximate Message Passing with Parameter Estimation
AMP was originally designed for linear measurement systems 28 . In order to use it for the nonlinear dipole inversion, we can linearize the measurement model in (3) iteratively using the zero-th and first orders of Taylor series 14 . Letting ( ) denote the susceptibility at the -th iteration, we have where ( ) = ⋅ exp ( ) ⋅ ( ) − is a relative constant in the ( + 1)-th iteration that depends on ( ) . We can then recover ( +1) using AMP based on the linear approximation in (8).
As shown by the factor graph in Fig. 2 , the distribution parameters are treated as random variables and jointly recovered with the susceptibility 29 . The variable nodes are represented by "○" and contain random variables. The factor nodes are represented by "■" and encodes probability distributions of the random variables. The messages about the variable distributions are passed and discussed among the factor nodes until a consensus is reached.
For example, we shall use the following notations to denote the messages between the -th variable node and the -th factor node Φ in the -th echo: • Δ →Φ denotes the message from to Φ ,

Signal prior
The 1st echo The -th echo Noise prior The factor graph of the QSM task: "○" represents the variable node, and "■" represents the factor node. In particular, is the signal prior parameter, = { 1 , 2 , 1 , 2 } contains the noise prior parameters.
• Δ Φ → denotes the message from Φ to , where ∈ {1, ⋯ , }, ∈ {1, ⋯ , } and ∈ {1, ⋯ , }. Both Δ →Φ and Δ Φ → are functions of the variable , and they are expressed in the "log" domain in this paper. Detailed expressions of the messages in the factor graph are given in Appendix A.
Currently the convergence of the AMP algorithm has been established for random Gaussian measurement matrices 28 . It is still an open question on how to extend the convergence analysis to general measurement matrices. In practice, the damping and mean-removal operations are used to ensure the convergence of AMP for non-Gaussian matrices 30,31 . For the ill-posed measurement operator in (8) that contains the dipole kernel, we perform the damping operation on to stabilize AMP. Let ( ) denote the damped susceptibility in the -th iteration, and ( +1) the undamped susceptibility in the ( + 1)-th iteration. The MAP estimation of the -th entry in is We can compute the damped solution ( +1) as follows where ∈ (0, 1] is the damping rate on the susceptibility. For the QSM task, we need to choose a small damping rate that is around 0.01.
In particular, we estimate the parameter of the sparse signal prior by maximizing its posterior: where the posterior distribution ( | ) can be computed via AMP.
The estimation of the noise prior parameters requires more work. Due to the ill-posedness of the measurement operator, the maximum-a-posteriori (MAP) method over-estimates the weight 2 of the second Gaussian component that is reserved for the noise outliers. Here we propose the following two-step procedure to estimate the Gaussian mixture weights 1 , 2 :

Perform a preliminary reconstruction with only one
Gaussian component  ( |0, 0 ) to model the noise.
2. Estimate the weights { 1 , 2 } of the two-component Gaussian mixture using the preliminary reconstruction 0 . We first calculate the residual error using 0 : We next use to approximate the noise , and treat the residue entries that fall outside [−3 0 , 3 0 ] as outliers.
The weights 1 and 2 can then be estimated aŝ where 1(⋅) is the indicator function.
With the estimated weightŝ 1 ,̂ 2 fixed, we then estimate the Gaussian mixture variances 1 , 2 by maximizing their posteriors: When we compute the MAP estimations of the parameters , 1 , 2 , we also use the damping operations to overcome the ill-posedness of the dipole kernel and stabilize AMP. In the ( + 1)-th iteration, we havê wherê ( ) ,̂ ( ) are the damped parameters in the previous ( )-th iteration,̂ ( +1) and̂ ( +1) are the undamped MAP estimations computed via (11) and (16) in the ( + 1)-th iteration, ∈ (0, 1] is the damping rate on the estimated parameters. For the QSM task, we can choose to be around 0.1. The proposed AMP with built-in parameter estimation (AMP-PE) approach is summarized in Algorithm 1.

METHODS
We compare the proposed AMP-PE approach with the stateof-the-art nonlinear L1-QSM 26 and MEDI 14 approaches on the 2019 QSM challenge dataset and the in vivo 3D brain datasets.
Specifically, the nonlinear L1-QSM and MEDI approaches adopt the total-variation regularization and solve the following nonlinear dipole inversion problems: MEDI: If convergence is reached < , then break. 8: end for 9: Estimate the weightŝ 1 ,̂ 2 based on the residue error . 10: Use the preliminary reconstruction 0 as the initialization, and perform the final reconstruction with the two-component Gaussian mixture noise model. 11: for = 1, ⋯ , do 12: Linearize the measurement model via (8). 13: Recover the susceptibilitŷ ( ) via AMP with parameter estimation. 14: If convergence is reached ( ) and break. 15: end for 16: Return̂ .
where (⋅) is the gradient operator that computes image gradients,  is a 0-1 morphology mask and only selects the gradients that correspond to non-edges, and are the regularization parameters. The nonlinear L1-QSM uses the magnitude image as the weighting matrix , whereas the MEDI approach constructŝ based on the reliability of the measurements and iteratively updates it according to the residual. We can see that the nonlinear L1-QSM and MEDI approaches use the 1 -norm and the 2 -norm respectively to construct the data-fidelity terms, leading to different interpretations of the noise distribution as discussed in the introduction.

Sim2Snr1 Dataset from QSM Challenge 2.0
The QSM reconstruction challenge 2.0 (Seoul 2019) provides the simulated Sim2Snr1 dataset with a groundtruth susceptibility map for evaluation 32 . The dataset contains an intrahemispheric calcification that creates strong dephasing effect and is thus more challenging. The simulation parameters are as follows: the main magnetic field 0 = 7 T, repetition time (TR) = 50 ms; echo time (TE) TE1/TE2/TE3/TE4 = 4/12/20/28 ms; echo spacing = 8 ms; the flip angle = 15°; field of view (FoV) = 164×205×205 mm 3 and 1 mm 3 isotropic voxels. The Gaussian noise is added to the complex data, producing a peak SNR of 100. Since only the brain tissues are used to simulate the field perturbations, the background field removal is thus not needed. As recommended in 26 , the nonlinear L1-QSM approach uses the provided frequency map with a simulated TE of 10 ms to recover the susceptibility map, the L-curve method is first used to determine the regularization parameter with 150 iterations, and QSM is finally recovered using the computed parameter with 2500 iterations. As recommended by the MEDI-toolbox 14 , the MEDI approach uses the provided frequency map with a simulated TE that equals the echo spacing (i.e. 8 ms) to recover the susceptibility map. The AMP-PE approach shall use the unwrapped phase images to recover the susceptibility map.
According to the established scientific practice 33,34,35 , it is forbidden to tune the parameter on the test set with respect to the groundtruth. The regularization parameter can be either tuned on a training set or determined using the heuristic Lcurve method. For the nonlinear L1-QSM approach, we use the L-curve method to determine the parameter as suggested by the FANSI-toolbox 36,26 . For the MEDI approach, the Lcurve method was not able to produce a working parameter in this case. We thus used the tuned parameter = 1000 as suggested by the MEDI-toolbox 14 .
Unlike the nonlinear L1-QSM and MEDI approaches that use the total-variation regularization, AMP-PE uses the sparse prior on the image wavelet coefficients to improve the image quality. We shall use the Daubechies wavelet family to obtain the sparse representation of the susceptibility map 22 . The orthogonal "db1-db10" wavelet bases are usually used, and the complexity of a wavelet basis increases with respect to its order. In particular, the db1 wavelet is suitable for capturing the low-frequency information in the image, and the higher-order wavelet bases are more suitable for capturing the high-frequency information. Although we would like to keep structural details of the susceptibility map that are encoded in the high-frequency bands, the streaking artifacts also contain predominantly high-frequency information. We thus need to find a proper basis that balances the trade-off between low and high frequencies.
• For the simulated Sim2Snr1 dataset, a lot of the regions in the groundtruth susceptibility map are piecewise constant, and contain more low-frequency information. The db1 wavelet is thus a better choice in this case.
• For the subsequent in vivo datasets, the susceptibility map contains a lot of fine structure details that could not be fully captured by the db1 wavelet, and the more complicated db2 wavelet basis is chosen instead. Higher-order wavelet bases such as db3 could capture more details from the streaking artifacts, and is not recommended for the QSM task in general.

In vivo 3D Brain Dataset
We acquired in vivo 3D brain data on a 3T MRI scanner (Prisma model, Siemens Healthcare, Erlangen, Germany), with written consent obtained from the healthy subjects before imaging. We also obtained clinical MR images from patients with hemorrhage as a retrospective chart review study, under the approval from the Institutional Review Board of Emory University. The patient scans were performed using a 3T MRI scanner (Tim Trio model, Siemens Healthcare, Erlangen, Germany). We reconstructed susceptiblity maps from both the healthy subject scans and patient scans from clinical cases with brain hemorrhage. The data were acquired using the GRE sequences that have the acquisition parameters: • Five healthy subject scans (H1-H5). We have the main magnetic field 0 = 3 T, the flip angle = 15°, the number of echoes = 4, the first echo time = 7.32 ms, echo spacing = 8.68 ms, TR = 38 ms, slice thickness = 0.7 mm, in-plane resolution = 0.6875 mm × 0.6875 mm, bandwidth per pixel = 260 Hz, and acquisition matrix size = 320 × 320 × 208.
In order to showcase the robustness of using complex exponential measurements exp( ), we also recovered susceptibility maps using the conventional linear measurements as follows: where the weighting matrix is chosen to be the magnitude image, > 0 is the regularization parameter.
The raw phase images were first unwrapped using Laplacian-based phase unwrapping 37 , and the background field was then removed using projection onto dipole fields (PDF) 16 . In order for the linear recovery, nonlinear L1-QSM, and MEDI approaches to achieve the best performances, the background-removed phase images were first converted to the corresponding local field map (in Hz) of each echo, the average local field map was then mapped to the phase image (in radian) at a chosen TE and used to recover the susceptibility map. The linear recovery approach used the first echo time TE1 as the chosen TE. As recommended by the FANSI toolbox, a TE of 20 ms was used for the in vivo reconstruction by the nonlinear L1-QSM. The MEDI approach used the echo spacing as the recommended TE. For the AMP-PE approach, the processed multi-echo phase images were used to recover the susceptibility as before, the db2 wavelet basis was adopted to capture enough high-frequency details in the in vivo reconstruction.

Sim2Snr1 Dataset from QSM challenge 2.0
Using the code provided by the QSM challenge 2.0, we computed the normalized root mean square error (NRMSE), ROI-based error metrics for tissue, blood and deep gray matter (DGM), calcification error metrics, structural similarity index measure (SSIM) and high frequency error norm (HFEN). The evaluation metric scores are shown in Table 1 , and the reconstructed QSMs are shown in Fig. 3 . In terms of global metrics, AMP-PE (db1) achieves the lowest NRMSE and the highest SSIM, MEDI achieves the lowest HFEM. For the rest region-specific metrics, each approach has its own strong suit: the nonlinear L1-QSM achieves the lowest calcification streaking and the lowest detrend NRMSE on blood and DGM, MEDI achieves the lowest detrend NRMSE on the brain tissue, AMP-PE (db1) achieves the lowest calcification error. The MEDI approach adopts the least-squares data-fidelity term that implies the Gaussian noise prior. On the other hand, the nonlinear L1-QSM approach adopts the 1 -norm data fidelity term that implies the Laplace noise prior. The AMP-PE approach adopts the Gaussian-mixture noise prior. In the calcification case, the latter two noise priors fit the longtailed noise distribution better than a single Gaussian-noise prior. Hence the nonlinear L1-QSM and AMP-PE approaches achieve lower NRMSE than the MEDI approach. The advantage of MEDI lies in its ability to incorporate structural information by penalizing non-edge gradients through the morphology mask . The edges in QSM recovered by MEDI are thus better preserved, which leads to a significantly lower HFEN. As mentioned earlier, it is forbidden to tune the parameters with respect to the groundtruth according to the established scientific practice 33,34,35 , the heuristic L-curve method is thus used to compute the parameter = 5.0119 −4 for L1-QSM. However, the L-curve method could not work with MEDI since the weighting matrix in (20) is updated iteratively with respect to the residue, which introduces additional variability to the data-fidelity term that changes the shape of the L-curve. We thus used the default parameter = 1000 for MEDI, which turns out to be working on both the simulated Sim2Snr1 and in vivo datasets.

In vivo 3D Brain Dataset
Recovered susceptibility maps from the healthy scan "H1" and hemorrhage scans "P1,P2" are shown in Fig. 4 -6 . Additional QSMs from "H2-H5" and "P3" are show in Fig. S1-S5 in the Supporting Information. For the nonlinear L1-QSM approach, although the L-curve method worked very well on the simulated Sim2Snr1 dataset, the parameters determined by the L-curve method over-regularized the reconstructions on the in vivo datasets, leading to the loss of finer details in recovered maps. As suggested in 36 , we shall resort to visual fine-tuning to find working parameters for the nonlinear L1-QSM when the computed parameters fail. However, visual fine-tuning would introduce subjective bias to recovered susceptibility maps inevitably. For the MEDI approach, the default parameter = 1000 was used. In order to ensure MEDI achieves its best performance, we need to double check the results through visual fine-tuning as well. The recovered maps produced under the default parameter in MEDI were generally consistent with the visual fine-tuning results. For the linear recovery approach, the parameter = 0.1 was determined through visual fine-tuning. On the other hand, the proposed AMP-PE did not require visual fine-tuning since  it estimated the parameters from the data automatically and adaptively.
For the healthy scans "H1, H2" shown in Fig. 4 and Fig. S1 in the Supporting Information, the linear recovery and AMP-PE approaches were better than the nonlinear L1-QSM and MEDI at recovering the structural details. For the healthy scans "H3-H5" shown in Fig. S2-S4 of the Supporting Information, the four approaches performed almost equally well.
The advantage of the three nonlinear approaches is revealed better on the more challenging hemorrhage scans. Compared to the linear recovery approach, the nonlinear L1-QSM, MEDI and AMP-PE approaches are more robust and did a better job at removing the streaking artifacts around the hemorrhage spot, thanks to complex exponential measurements they used. For the patient scan "P1" shown in Fig.  5 , AMP-PE revealed more structural details than the nonlinear L1-QSM and removed the streaking artifacts around the hemorrhage spot better than MEDI, indicating that the Gaussian-mixture noise prior is a better fit in this case. For the patient scan "P2" shown in Fig. 6 , MEDI and AMP-PE performed almost equally well. However, the nonlinear L1-QSM failed to remove the streaking artifacts caused by hemorrhage spots. For the patient scan "P3" shown in Fig. S5 in the Supporting Information, the nonlinear L1-QSM, MEDI and AMP-PE all performed equally well. In general, AMP-PE  is more robust than the nonlinear L1-QSM and MEDI across the healthy and hemorrhage scans.

DISCUSSION
In the simulated Sim2Snr1 dataset from the QSM challenge, only the brain tissue within the ROI is used to generate the field variations and there is no need for background field removal. It is thus beneficial to enforce the binary brain mask  on by setting the susceptibility outside the ROI to 0 during optimization. In practice, however, the background field removal would not be perfect for the in vivo reconstruction. Usually there is still some background field left in the ROI. If we enforce the brain mask on as in the simulated case, the residual background field would lead to incorrect susceptibility variations in the brain tissue. The left plot in Fig. 7 (a) shows the recovered susceptibility map from AMP-PE when we enforce the brain mask on during optimization. The right plot in Fig. 7 (a) shows the recovered map when we allow the susceptibility outside the ROI to be nonzero during optimization, and we only use the brain mask on the completed reconstruction to extract the brain in the last step. By allowing the susceptibility outside the ROI to account for the residual background field, we can recover a much better susceptibility map in Fig. 7 (a). The L-curve method has been used previously to compute the regularization parameter. Due to its heuristic nature, it is not robust enough against varying signal and noise conditions. The L-curve was originally proposed for the leastsquares data-fidelity term with Tikhonov regularization (i.e. the squared 2 -norm) and achieved relatively robust performance 27 . However, when it was used for the 1 -norm datafidelity term with the 1 -norm regularization in the nonlinear L1-QSM, the computed parameter might not be close to the optimal parameter, since the statistical property of the 1 -norm is different from that of the squared 2 -norm. To address this issue, the L-curve method was combined with visual finetuning in 26 for in vivo reconstructions. This requires human intervention and introduces subjective bias to the results. In this paper we propose the Bayesian approach for QSM reconstruction that allows us to estimate the parameters adaptively and automatically. The estimated parameters maximize their posterior distributions, which avoids the subjective bias introduced by visual fine-tuning.
In order for AMP-PE to use the sparse signal prior to recover the susceptibility map, the sparsifying transform on the image needs to be invertible. The total-variation transform is not invertible and thus could not be used. Here we use the wavelet transform to get the sparse representation of image in the form of wavelet coefficients. Fig. 7 (b) shows the recovered susceptibility maps using AMP-PE with different wavelet bases. The db1 basis could not capture enough high-frequency information from the image, and the recovered map shows pixelation effect. As the order of wavelet basis gets higher, more high-frequency information can be captured. However, the db6 basis captures the high-frequency information from both the image and streaking artifacts (see the coronal views in Fig. 7 (b)). We need to balance the trade-off between the captured low-and high-frequency information. In practice, it is recommended that the db2 basis be chosen for the in vivo QSM task. On the other hand, we should note that the db2 basis might not be the optimal choice for other MRI tasks. For example, the 1 and * 2 mappings do not involve the removal of streaking artifacts, and the db6 basis was used in 38 to obtain the optimal reconstruction of * 2 map. As shown in Fig. 1 (b), we adopt the two-component Gaussian-mixture distribution to model the noise in QSM. The second component is used to take care of the noise outliers that are caused by phase unwrapping errors or brain hemorrhage and calcification. Fig. 8 (a) shows the comparison of recovered susceptibility maps in the case of phase unwrapping errors where we use the (one-component) Gaussian distribution and the (two-component) Gaussian-mixture to model the noise. We can see that the AMP-PE with a Gaussian-mixture noise prior is better at removing the streaking artifact than the AMP-PE with a single-Gaussian noise prior. As discussed earlier in Section 2.2, AMP-PE overestimates the weight 2 of the second Gaussian component due to the ill-posed dipole kernel. In order to address this problem, we proposed a two-step procedure where 1 , 2 are first estimated based on a preliminary reconstruction and then fixed in the final reconstruction. The left plot in Fig. 8 (b) shows the recovered map using the standard procedure where 1 , 2 are freely updated; and the right plot in Fig. 8 (b) shows the recovered map using the proposed two-step procedure where the streaking artifacts are removed better and the images are cleaner.
The convergence of AMP can be established for random Gaussian measurement matrices. The damping and meanremoval operations are used to ensure AMP's convergence for non-Gaussian measurement matrices. Since the means of rows of the transform operator in (1) are approximately zero, the mean of the resulting measurement operator in (8) is also approximately zero, we thus do not need to perform the mean-removal operation in this case. However, the distribution ( ) of the entries in is not Gaussian. Take the simple case where the size of is 10 3 × 10 3 for example, the normal probability plot in Fig. 9 (a) shows that ( ) is both left-and right-skewed with discontinuity in the distribution, the histogram in Fig. 9 (b) also corroborates the observation that ( ) is quite different from the Gaussian distribution. As discussed in Section 2.2, we use the damping operations on the estimated susceptibility and parameters to achieve the convergence of AMP. In addition, since the QSM reconstruction is a nonconvex problem, the initialization also plays an important role in stabilizing the AMP algorithm. We found that the simple least-squares solution could be computed efficiently and used to initialize AMP.

CONCLUSION
We propose a Bayesian formulation to recover the susceptibility map from nonlinear complex exponential measurements that are robust to phase errors. We adopt a customized two-component Gaussian mixture noise prior to model the long-tailed noise distribution that appears in cases of brain hemorrhage or calcification. Combined with the sparsitypromoting Laplace signal prior, we recover the susceptibility using the proposed AMP with automatic and adaptive parameter estimation (AMP-PE). In vivo experiments showed that AMP-PE is more robust and better than the nonlinear L1-QSM and MEDI at preserving structural details and removing Gaussian noise prior

FIGURE 9
The distribution of the entries in the transform operator .

SUPPORTING INFORMATION
The following supporting information is available as part of the online article: Figure S1. Healthy scan (H2): recovered susceptibility maps using the linear recovery approach, nonlinear L1-QSM, MEDI, and AMP-PE with the db2 wavelet basis. Figure S2. Healthy scan (H3): recovered susceptibility maps using the linear recovery approach, nonlinear L1-QSM, MEDI, and AMP-PE with the db2 wavelet basis. Figure S3. Healthy scan (H4): recovered susceptibility maps using the linear recovery approach, nonlinear L1-QSM, MEDI, and AMP-PE with the db2 wavelet basis. Figure S4. Healthy scan (H5): recovered susceptibility maps using the linear recovery approach, nonlinear L1-QSM, MEDI, and AMP-PE with the db2 wavelet basis. Figure S5. Patient scan (P3): recovered susceptibility maps using the linear recovery approach, nonlinear L1-QSM, MEDI, and AMP-PE with the db2 wavelet basis.

A MESSAGES EXCHANGED BETWEEN THE VARIABLE AND FACTOR NODES
In the following we derive the messages exchanged on the factor graph in Fig. 2 : • We have the following messages passed from Φ to in the ( + 1)-th iteration.
• We further have the following messages passed from to Φ in the ( + 1)-th iteration. (A10)