Y Kee

L Zhou

Z Liu

A Dimov

J Cho

J K Seo

Y Wang
email: yiwang@med.cornell.edu

Quantitative Susceptibility Mapping (QSM) Algorithms: Mathematical Rationale and Computational Implementations Youngwook Kee, Zhe Liu, Liangdong Zhou, Alexey Dimov, Junghun Cho, Ludovic de Rochefort, Jin Keun Seo, and Yi Wang*, Fellow, IEEE Abstract-Quantitative susceptibility mapping (QSM) solves the magnetic field-to-magnetization (tissue susceptibility) inverse problem under conditions of noisy and incomplete field data acquired using magnetic resonance imaging. Therefore, sophisticated algorithms are necessary to treat the ill-posed nature of the problem and are reviewed here. The forward problem is typically presented as an integral form, where the field is the convolution of the dipole kernel and tissue susceptibility distribution. This integral form can be equivalently written as a partial differential equation (PDE). Algorithmic challenges are to reduce streaking and shadow artifacts characterized by the fundamental solution of the PDE. Bayesian maximum a posteriori (MAP) estimation can be employed to solve the inverse problem, where morphological and relevant biomedical knowledge (specific to the imaging situation) are used as priors. As the cost functions in Bayesian QSM framework are typically convex, solutions can be robustly computed using a gradient-based optimization algorithm. Moreover, one can not only accelerate Bayesian QSM, but also increase its effectiveness at reducing shadows using prior knowledge based precondition. Improving the efficiency of QSM is under active development, and a rigorous analysis of preconditioning needs to be carried out for further investigation.

Index Terms-Bayesian image reconstruction, deconvolution, nonsmooth convex optimization, partial differential equations, preconditioning, quantitative susceptibility mapping (QSM).

I. INTRODUCTION

Q UANTITATIVE susceptibility mapping (QSM) [1, 2] is a noninvasive magnetic resonance imaging (MRI) method that enables quantitative investigation of the tissue magnetic susceptibility. QSM is based on the post-processing of the magnetic field map derived from the phase data of gradient echo (GRE) MRI that is readily available on most MRI scanners [START_REF] Bernstein | Handbook of MRI pulse sequences[END_REF]. QSM has been applied to a wide range of biomedical problems [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF][4][START_REF] Reichenbach | Quantitative susceptibility mapping: concepts and applications[END_REF][START_REF] Liu | Quantitative susceptibility mapping: contrast mechanisms and clinical applications[END_REF][START_REF] Eskreis-Winkler | The clinical utility of qsm: disease diagnosis, medical management, and surgical planning[END_REF][START_REF] Wang | Clinical quantitative susceptibility mapping (qsm): Biometal imaging and its emerging roles in patient care[END_REF][START_REF] Deistung | Overview of quantitative susceptibility mapping[END_REF] including: 1) demyelination, inflammation, and iron overload in multiple sclerosis [START_REF] Wisnieff | Quantitative susceptibility mapping (qsm) of white matter multiple sclerosis lesions: interpreting positive susceptibility and the presence of iron[END_REF][START_REF] Chen | Quantitative susceptibility mapping of multiple sclerosis lesions at various ages[END_REF][START_REF] Harrison | Lesion heterogeneity on high-field susceptibility mri is associated with multiple sclerosis severity[END_REF][START_REF] Li | Magnetic susceptibility contrast variations in multiple sclerosis lesions[END_REF][START_REF] Zhang | Quantitative susceptibility mapping and r2* measured changes during white matter lesion development in multiple sclerosis: myelin breakdown, myelin debris degradation and removal, and iron accumulation[END_REF][START_REF] Zhang | Longitudinal change in magnetic susceptibility of new enhanced multiple sclerosis (ms) lesions measured on serial quantitative susceptibility mapping (qsm)[END_REF][START_REF] Zhang | Magnetic susceptibility from quantitative susceptibility mapping can differentiate new enhancing from nonenhancing multiple sclerosis lesions without gadolinium injection[END_REF][START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF], 2) neurodegeneration and iron overload in Alzheimer's disease [START_REF] Acosta-Cabronero | In vivo quantitative susceptibility mapping (qsm) in alzheimer's disease[END_REF][START_REF] Hwang | Texture analyses of quantitative susceptibility maps to differentiate alzheimer's disease from cognitive normal and mild cognitive impairment[END_REF][START_REF] Moon | Patterns of brain iron accumulation in vascular dementia and alzheimers dementia using quantitative susceptibility mapping imaging[END_REF], Parkinson's disease [START_REF] Azuma | Lateral asymmetry and spatial difference of iron deposition in the substantia nigra of patients with parkinson disease measured with quantitative susceptibility mapping[END_REF][START_REF] Barbosa | Quantifying brain iron deposition in patients with parkinson's disease using quantitative susceptibility mapping, r2 and r2[END_REF][START_REF] Du | Quantitative susceptibility mapping of the midbrain in parkinson's disease[END_REF][START_REF] Guan | Regionally progressive accumulation of iron in parkinson's disease as measured by quantitative susceptibility mapping[END_REF][START_REF] Guan | Influence of regional iron on the motor impairments of parkinson's disease: A quantitative susceptibility mapping study[END_REF][START_REF] He | Region-specific disturbed iron distribution in early idiopathic parkinson's disease measured by quantitative susceptibility mapping[END_REF][START_REF] Ide | Internal structures of the globus pallidus in patients with parkinson's disease: evaluation with quantitative susceptibility mapping (qsm)[END_REF][START_REF] Langkammer | Quantitative susceptibility mapping in parkinson's disease[END_REF][START_REF] Lotfipour | High resolution magnetic susceptibility mapping of the substantia nigra in parkinson's disease[END_REF][START_REF] Murakami | Usefulness of quantitative susceptibility mapping for the diagnosis of parkinson disease[END_REF][START_REF] Chandran | Magnetic resonance imaging of the subthalamic nucleus for deep brain stimulation[END_REF][START_REF] Liu | Improved subthalamic nucleus depiction with quantitative susceptibility mapping[END_REF][START_REF] Schäfer | Direct visualization of the subthalamic nucleus and its iron distribution using high-resolution susceptibility mapping[END_REF] and Huntington's disease [START_REF] Van Bergen | Quantitative susceptibility mapping suggests altered brain iron in premanifest huntington disease[END_REF], 3) changes in metabolic oxygen consumption [START_REF] Haacke | Susceptibility mapping as a means to visualize veins and quantify oxygen saturation[END_REF][START_REF] Kudo | Oxygen extraction fraction measurement using quantitative susceptibility mapping: comparison with positron emission tomography[END_REF][START_REF] Xu | Flow compensated quantitative susceptibility mapping for venous oxygenation imaging[END_REF][START_REF] Zhang | Quantitative mapping of cerebral metabolic rate of oxygen (cmro2) using quantitative susceptibility mapping (qsm)[END_REF][START_REF] Zhang | Quantitative susceptibility mapping-based cerebral metabolic rate of oxygen mapping with minimum local variance[END_REF][START_REF] Zhang | Cerebral metabolic rate of oxygen (cmro2) mapping with hyperventilation challenge using quantitative susceptibility mapping (qsm)[END_REF][START_REF] Fan | Quantitative oxygenation venography from mri phase[END_REF], 4) hemorrhage including microhemorrhage and blood degradation [START_REF] Chang | Quantitative susceptibility mapping of intracerebral hemorrhages at various stages[END_REF][START_REF] Chen | Intracranial calcifications and hemorrhages: characterization with quantitative susceptibility mapping[END_REF][START_REF] Sun | Quantitative susceptibility mapping using a superposed dipole inversion method: application to intracranial hemorrhage[END_REF][START_REF] Klohs | Detection of cerebral microbleeds with quantitative susceptibility mapping in the arcabeta mouse model of cerebral amyloidosis[END_REF][START_REF] Liu | Cerebral microbleeds: burden assessment by using quantitative susceptibility mapping[END_REF][START_REF] Liu | Imaging cerebral microhemorrhages in military service members with chronic traumatic brain injury[END_REF], 5) bone mineralization [START_REF] Dimov | Bone quantitative susceptibility mapping using a chemical speciesspecific r2* signal model with ultrashort and conventional echo data[END_REF], and 6) drug delivery using magnetic nanocarriers [START_REF] Liu | Unambiguous identification of superparamagnetic iron oxide particles through quantitative susceptibility mapping of the nonlinear response to magnetic fields[END_REF][START_REF] Wong | Visualizing and quantifying acute inflammation using icam-1 specific nanoparticles and mri quantitative susceptibility mapping[END_REF]. QSM of dominant scalar susceptibility sources in tissue (including iron and calcium) has now become reasonably robust and accurate for routine use in research and practice [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF][START_REF] Eskreis-Winkler | The clinical utility of qsm: disease diagnosis, medical management, and surgical planning[END_REF].

Algorithmically, QSM solves a challenging ill-posed inverse problem. QSM was coined from a breakthrough in 2008 using the Bayesian maximum a posteriori (MAP) estimation theory [START_REF] De Rochefort | Quantitative susceptibility map reconstruction from mr phase data using bayesian regularization: validation and application to brain imaging[END_REF][START_REF] De Rochefort | Quantitative mr susceptibility mapping using piece-wise constant regularized inversion of the magnetic field[END_REF][START_REF] Kressler | Estimation of sparse magnetic susceptibility distributions from mri using non-linear regularization[END_REF][START_REF] Kressler | Nonlinear regularization for per voxel estimation of magnetic susceptibility distributions from mri field maps[END_REF][START_REF] Wang | Magnetic source mri: a new quantitative imaging of magnetic biomarkers[END_REF]]. Yet, it remains to be fully developed and there is plenty of room for further optimization. To establish this milestone and characterize associated challenges, this review paper defines the forward problem in QSM as a partial differential equation (PDE). Then, the difficulties of the field-to-susceptibility inversion problem are systematically analyzed based on general solutions for the PDE. This analysis emphasizes the limitations of early attempts [START_REF] Li | Quantifying arbitrary magnetic susceptibility distributions with mr[END_REF][START_REF] Morgan | Efficient solving for arbitrary susceptibility distributions using residual difference fields[END_REF][START_REF] Haacke | Imaging iron stores in the brain using magnetic resonance imaging[END_REF][START_REF] Shmueli | Magnetic susceptibility mapping of brain tissue in vivo using mri phase data[END_REF] compared to Bayesian approaches. Bayesian QSM methods optimally estimate the underlying tissue susceptibility distribution from noisy and incomplete magnetic field data (likelihood) using tissue structural information (a priori). To this end, variational methods are used to minimize a cost function associated with the negative posterior, where the optimal solution (MAP estimate) corresponds to a minimizer of the cost function. In Bayesian QSM, such a cost function is typically given as a convex objective so that one can use solvers developed in modern convex optimization [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF] for efficient computation of a solution. This principled approach avoids the drawbacks inherent in the non-Bayesian attempts and has laid the foundation for QSM, generating a great deal of interest in studying signal phase, the often-neglected other half of the complex MRI signal [START_REF] Li | Reducing the object orientation dependence of susceptibility effects in gradient echo mri through quantitative susceptibility mapping[END_REF][START_REF] Liu | Susceptibility tensor imaging[END_REF][START_REF] Schweser | Differentiation between diamagnetic and paramagnetic cerebral lesions based on magnetic susceptibility mapping[END_REF][START_REF] Wharton | Whole-brain susceptibility mapping at high field: a comparison of multiple-and single-orientation methods[END_REF][START_REF] Wharton | Susceptibility mapping in the human brain using threshold-based k-space division[END_REF][START_REF] Li | Quantitative susceptibility mapping of human brain reflects spatial variation in tissue composition[END_REF][START_REF] Liu | High-field (9.4 t) mri of brain dysmyelination by quantitative mapping of magnetic susceptibility[END_REF][START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF][START_REF] Liu | Morphology enabled dipole inversion (medi) from a single-angle acquisition: comparison with cosmos in human brain imaging[END_REF][START_REF] Bilgic | Mri estimates of brain iron concentration in normal aging using quantitative susceptibility mapping[END_REF][START_REF] Langkammer | Quantitative susceptibility mapping (qsm) as a means to measure brain iron? a post mortem validation study[END_REF][START_REF] Li | Mapping magnetic susceptibility anisotropies of white matter in vivo in the human brain at 7t[END_REF][START_REF] Liu | Morphology enabled dipole inversion for quantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map[END_REF][START_REF] Liu | Accuracy of the morphology enabled dipole inversion (medi) algorithm for quantitative susceptibility mapping in mri[END_REF][START_REF] Schweser | Quantitative susceptibility mapping for investigating subtle susceptibility variations in the human brain[END_REF][START_REF] Wu | Whole brain susceptibility mapping using compressed sensing[END_REF] There are a number of existing review papers on QSM, which can be categorized as follows: 1) General reviews [4,[START_REF] Reichenbach | Quantitative susceptibility mapping: concepts and applications[END_REF][START_REF] Deistung | Overview of quantitative susceptibility mapping[END_REF][START_REF] Schweser | Foundations of mri phase imaging and processing for quantitative susceptibility mapping (qsm)[END_REF], 2) phase processing methods [START_REF] Robinson | An illustrated comparison of processing methods for mr phase imaging and qsm: combining array coil signals and phase unwrapping[END_REF][START_REF] Schweser | An illustrated comparison of processing methods for phase mri and qsm: removal of background field contributions from sources outside the region of interest[END_REF], 3) biophysics [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF], and 4) clinical applications [START_REF] Eskreis-Winkler | The clinical utility of qsm: disease diagnosis, medical management, and surgical planning[END_REF][START_REF] Wang | Clinical quantitative susceptibility mapping (qsm): Biometal imaging and its emerging roles in patient care[END_REF]. Unlike these, this review focuses on mathematical and computational aspects of the deconvolution problem in QSM. Specifically, the source of artifacts in QSM images is identified as dipoleincompatible fields that deviates from the field generated by dipoles. Then, the artifacts are characterized as streaking and shadow by the fundamental solution of a wave-type operator. This characterization establishes a rationale for the design principles of a regularization term in the Bayesian QSM methods, implying that early attempts based on k-space truncation causes additional dipole-incompatibility in the field. These Bayesian QSM methods require numerical optimization whose implementation details and preconditioning techniques are provided here. Lastly, the review discusses remaining challenges in this field. This review will serve as a unique reference to QSM algorithm developers working on numerical optimization.

II. FORWARD PROBLEM-BIOPHYSICS MODEL

In an MRI scanner with a main static magnetic field B 0 = (0, 0, B 0 ), tissues gain a magnetization distribution M(r), where r = (x, y, z) ∈ R 3 is a position vector. This magnetization generates a magnetic field perturbation B(r) superpositioning on B 0 according to magnetostatic macroscopic Maxwell's equations [START_REF] Jackson | Electrodynamics[END_REF]

∇ • B(r) = 0 (1) ∇ × B(r) = µ 0 ∇ × M(r), (2)
where µ 0 is the vacuum permittivity. Using the Lorentz sphere correction model [START_REF] Lorentz | The Theory of Electrons and Its Applications to the Phenomena of Light and Radiant Heat: A Course of Lectures Delivered in Columbia University[END_REF][START_REF] Durrant | Magnetic susceptibility: further insights into macroscopic and microscopic fields and the sphere of lorentz[END_REF], the macroscopic field B(r) can be related to the microscopic field b (r) experienced by the spins of water protons in tissue, the MRI signal generators, as follows

b (r) = B(r) - 2 3 µ 0 M(r). (3)
MRI signal phase can be used for this field estimation, which includes nonlinear data fitting and field unwrapping as summarized in recent reviews [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF][START_REF] Schweser | Foundations of mri phase imaging and processing for quantitative susceptibility mapping (qsm)[END_REF]. Applying the curl operation to both sides of (2) and using ( 1) and ( 3), we obtain

-∆(b -B 0 ) = µ 0 ∇(∇ • M) - 1 3 ∆M , (4)
where ∆ is the Laplacian operator. The Fourier transform of (4), denoted as F, leads us to

|k| 2 F[b -B 0 ](k) = µ 0 k (k • F[M](k)) - 1 3 |k| 2 F[M](k) , (5)
where k = (k x , k y , k z ) ∈ R 3 is a spatial frequency vector. Isotropic biological tissue concerned here has weak magnetic susceptibility, i.e., |χ| 1, hence its magnetization can be expressed as M = (0, 0, M z ), where M z (B 0 /µ 0 )χ. In the following, we define the relative difference field (RDF) b, which is referred to as field, total field1 , or field inhomogeneity:

b(r) := b z (r) -B 0 B 0 ,
where b z denotes the z-component of b . We can rewrite [START_REF] Reichenbach | Quantitative susceptibility mapping: concepts and applications[END_REF] in terms of the z-component as follows

|k| 2 B(k) = k 2 x + k 2 y -2k 2 z 3 X(k), (6)
where

B(k) := F[b](k) and X(k) := F[χ](k).
For all k = 0, we can rewrite the expression (6) as

B(k) = 1 3 - k 2 z |k| 2 X(k) := D(k)X(k), (7)
where

D(k) = 1/3 -k 2 z /|k| 2 is the dipole kernel in k-space. Equivalently, b(r) = lim 0 |r-r |> d(r -r )χ(r ) dr , (8)
where d denotes the dipole kernel in r-space (image space) defined as

d(r) := F -1 [D](r) = 2z 2 -x 2 -y 2 4π|r| 5 .
This biophysics model (forward problem) describes that an RDF is generated by the dipole kernel convolved with the magnetic susceptibility distribution χ.

III. INVERSE PROBLEM-A WAVE PROPAGATION VIEW

Current acquisition and phase processing methods allow us to estimate an RDF from GRE phase data with high precision. The problem is to recover χ from the RDF b by fitting this with the forward model [START_REF] Eskreis-Winkler | The clinical utility of qsm: disease diagnosis, medical management, and surgical planning[END_REF] in k-space or (8) in image space. This leads to solving a deconvolution problem (inverse problem) whose nature is ill-posed because of a conic surface defined by

Γ 0 := {k ∈ R 3 | k 2 x + k 2 y -2k 2 z = 0}.
The surface Γ 0 is spanned by the magic angle ( 54.7 • ) with respect to the k z -direction as shown in Fig. 1. As can be seen in ( 6), the k-space points in Γ 0 cause division by zero when the deconvolution is performed in k-space.

This ill-posed inverse problem can be investigated in image space. Let us consider (4), the partial differential equation (PDE) associated with [START_REF] Liu | Quantitative susceptibility mapping: contrast mechanisms and clinical applications[END_REF], as follows:

P (∂)χ(r) = -∆b(r), (9)
where P (∂) is a wave-type operator (propagator) given as 9) is the inverse Fourier transform of (6), the inverse problem of recovering X(k) in (6) may be viewed as the kspace representation of the PDE [START_REF] Deistung | Overview of quantitative susceptibility mapping[END_REF]. This PDE is a wave-type equation in χ, where the z-axis is considered time, but is a Laplace equation in b. As noted in [START_REF] Evans | Partial Differential Equations[END_REF], one good strategy for investigating any PDE is to identify its explicit solutions; this can be achieved by deriving a fundamental solution. Since every constant coefficient linear PDE on R n has a fundamental solution [START_REF] Stein | Functional Analysis: Introduction to Further Topics in Analysis[END_REF], a solution of (9) denoted as χ can be expressed as

P (∂) := - 1 3 ∆ + ∂ 2 ∂z 2 . Because (
χ = g * (-∆b) = - R 3 g(r -r )∆ r b(r ) dr , (10)
where r = (x , y , z ) ∈ R 3 is a position vector and

g(r) = 3 4π √ z 2 -2(x 2 +y 2 ) if r ∈ Υ; 0 otherwise (11)
is the fundamental solution of the differential operator

P (∂) with F[g](k) = 3 4π 2 (k 2
x +k 2 y -2k 2 z ) =: 1 P (k) as displayed in Fig. 1 and Υ defined as

Υ := {r ∈ R 3 | z 2 -2(x 2 + y 2 ) > 0}. g(r)
1/P (k) D(k) Observe that the conic surface of the fundamental solution is perpendicular to that of D(k).

The detailed derivation of (11) can be found in Appendix 0.3 in [84]. Note that Υ defines a double-cone in R 3 whose conic surface is spanned by the angle 35.3 • , the complimentary magic angle. The fundamental solution g is shown in Fig. 1.

We can also obtain from [START_REF] Eskreis-Winkler | The clinical utility of qsm: disease diagnosis, medical management, and surgical planning[END_REF], the Fourier transform of (9), a solution in k-space as follows:

X(k) = B(k)/D(k), k ∈ R 3 \ {0}. (12)
How do we interpret these two explicit forms-( 10) and ( 12)-of a solution to the field-to-source inversion problem, particularly when there is some perturbation, e.g., noise, in the RDF b?

A. Dipole-Compatible Field Data

Suppose that RDF data is "purely" generated by dipoles (denoted here as b c (r), B c (k),), i.e.,

B c (k) = P (k)X c (k) 4π 2 |k| 2 = D(k)X c (k), k ∈ R 3 \ {0}, (13)
and we assign some constant to B c (0). Then, the field B c is said to be dipole compatible. For any B c , one can perfectly recover X c for k ∈ R 3 \ Γ 0 by the relation

X c (k) = B c (k) D(k) = 4π 2 |k| 2 B c (k) P (k) . ( 14
)
For k ∈ Γ 0 \ {0}, X c is recovered by the following limit [START_REF] Choi | Inverse problem in quantitative susceptibility mapping[END_REF]:

X c (k) = lim h 0 4π 2 |k + he z | 2 B c (k + he z ) P (k + he z ) (15)
where e z = (0, 0, 1). It should be noted that neither lim k→0 X c (k) nor X c (0) exist. This can easily be shown by considering two different paths (such as

k x = k y = 0, k z → 0, then 1/D → -3/2 or k x = k z = 0, k y → 0, then 1/D → 3).
Therefore, X c is only recovered up to a constant shift, i.e., X c (0), by using the direct reconstruction formulas [START_REF] Zhang | Quantitative susceptibility mapping and r2* measured changes during white matter lesion development in multiple sclerosis: myelin breakdown, myelin debris degradation and removal, and iron accumulation[END_REF] and [START_REF] Zhang | Longitudinal change in magnetic susceptibility of new enhanced multiple sclerosis (ms) lesions measured on serial quantitative susceptibility mapping (qsm)[END_REF]. Note that X c (0) is a symbolic notation for some constant shift in image space. The performance of ( 14) and ( 15) is shown in Fig. 2 where we assume that X c (0) = 0. Remark 1: That X c (k) is undefined at the center of kspace can be also checked from the PDE (9): We see that Hence, susceptibility values are often reported with respect to an estimate of reference tissue susceptibility. This is also reflected in (8) that lim 0 |r-r |> Cd(r ) dr = 0.

B. Dipole-Incompatible Field Data-Streaking and Shadow

Let

Γ B 0 be the set {k ∈ Γ 0 | B(k) = 0} ⊆ Γ 0 \ {0}. Then, regarding the dipole-compatible part, we see Γ Bc 0 = ∅ since B c (k) = D(k)X c (k) where D(k) = 0 in Γ 0 \ {0}. Let B i be a measured field such that Γ Bi 0 = ∅;
then we call such a field dipole-incompatible as the field cannot be expressed by the dipole model [START_REF] Li | Magnetic susceptibility contrast variations in multiple sclerosis lesions[END_REF]. Consider a dipole-compatible field with additive noise, B i (k) = D(k)X(k) + N (k), where N (k) is (Gaussian) noise in k-space. Then, we express the underlying susceptibility χ i in terms of the fundamental solution g in ( 10) and [START_REF] Chen | Quantitative susceptibility mapping of multiple sclerosis lesions at various ages[END_REF] as follows:

χ i (r) = (g * (-∆b i ))(r) = R 3 -3∆ r b i (r ) 4π E(r, r ) dr , (16)
where E(r, r

) := (z -z ) 2 -2((x -x ) 2 + (y -y ) 2 ) > 0
defines the double-cone Υ above and below the incompatible source -∆ r b i (r ). Strong propagation of -∆ r b i (r ) is present along the double-cone surface where E(r, r ) vanishes, which manifests as streaking in sagittal and coronal views and ringing in the axial view.

To see this more concretely, suppose that an incompatible source is expressed as a linear combination of the Dirac delta functions, i.e., -∆b i (r) = N j=1 c j δ(r -r j ), where c j ∈ R and r j = (x j , y j , z j ) ∈ R 3 is a position vector. Therefore, χ i is expressed as a linear combination of Green's functions, where each of the functions defines the double-cone Υ located at r 1 , . . . , r N as follows:

χ i (r) = g(r) * N j=1 c j δ(r -r j ) = N j=1 c j g(r -r j ).
Notice that the singular support of the double-cone Υ defined as the set of r ∈ R 3 such that z 2 -2(x 2 + y 2 ) = 0 forms a conic surface spanned by the complimentary magic angle ( 35.3 • ) along the z-axis. This conic surface is associated with Γ 0 defined in k-space, a zero cone surface spanned by the Remark 2: Observe that E(r, r ) in the denominator of (16) that causes propagating artifacts does not come into play when the field is dipole-compatible. As can be seen in the direct reconstruction formulas ( 14) and (15), 1/P (k) = F[g](k) is effectively canceled out when the field is dipole-compatible. This observation reassures the existence of streaking-free solutions proved in [84,[START_REF] Choi | Inverse problem in quantitative susceptibility mapping[END_REF].

b i in axial b i in sagittal b i in coronal χ i in axial χ i in sagittal χ i in coronal χ c + χ i in axial χ c + χ i in sagittal ∆b i in sagittal
Not only is streaking explained by [START_REF] Zhang | Magnetic susceptibility from quantitative susceptibility mapping can differentiate new enhancing from nonenhancing multiple sclerosis lesions without gadolinium injection[END_REF], but also the expression allows χ i to have other types of "apparent artifacts". Suppose that ∆b i is smoothly distributed across the image volume; then, χ i is expected to exhibit slowly varying spatial artifacts that manifest as dark shadows as demonstrated in Fig. 4. Note that if the sign of ∆b i is flipped, bright clouds will be seen instead. Although these artifacts appear different with low-frequency components as opposed to streaking with highfrequency, it is crucial to keep in mind that the expression ( 16) remains the same in both cases. In other words, dark shadows are simply a superposition of the double-cones Υ computed by [START_REF] Zhang | Magnetic susceptibility from quantitative susceptibility mapping can differentiate new enhancing from nonenhancing multiple sclerosis lesions without gadolinium injection[END_REF]. How does a solution expressed as [START_REF] Wisnieff | Quantitative susceptibility mapping (qsm) of white matter multiple sclerosis lesions: interpreting positive susceptibility and the presence of iron[END_REF] behave? If a given RDF is dipole-compatible, perfect recovery is possible up to a constant shift using ( 14) and [START_REF] Zhang | Longitudinal change in magnetic susceptibility of new enhanced multiple sclerosis (ms) lesions measured on serial quantitative susceptibility mapping (qsm)[END_REF]. If the RDF is dipoleincompatible, streaking and shadow artifacts will be present according to [START_REF] Zhang | Magnetic susceptibility from quantitative susceptibility mapping can differentiate new enhancing from nonenhancing multiple sclerosis lesions without gadolinium injection[END_REF] and the shape of dipole-incompatible source distributions.

b i in axial b i in sagittal b i in coronal χ i in axial χ i in sagittal χ i in coronal χ c + χ i in axial χ c + χ i in sagittal ∆b i in sagittal

C. Minimizing Streaking and Shadow: Early Attempts based on k-space truncation versus Bayesian MAP Estimation

Dipole-incompatible field data that deviates from the dipole model ( 13) includes noise, anisotropic sources, chemical shift, or discretization error. Projecting or filtering out the incompatible part directly from the GRE data is very difficult, because 1) noise is omnipresent in practice, 2) anisotropic sources cannot be measured reliably in a single-angle acquisition, and 3) discretization is indispensable for numerical computation. Consequently, the formulas ( 14) and ( 15) cannot be used in practice, as their use on the MRI estimated field data causes severe artifacts from the dipole-incompatible part.

Early attempts to solve the inverse problem are based on kspace truncation to modify (13) near the zero-cone Γ 0 where 1/D(k) → ∞ for k near Γ 0 . These methods include truncated k-space division (TKD) [START_REF] Shmueli | Magnetic susceptibility mapping of brain tissue in vivo using mri phase data[END_REF], truncated singular value decomposition (TSVD) or Tikhonov-regularized minimal norm (MN) [START_REF] Kressler | Nonlinear regularization for per voxel estimation of magnetic susceptibility distributions from mri field maps[END_REF][START_REF] Wharton | Whole-brain susceptibility mapping at high field: a comparison of multiple-and single-orientation methods[END_REF], iterative susceptibility weighted imaging and susceptibility mapping (iSWIM) [START_REF] Tang | Improving susceptibility mapping using a threshold-based k-space/image domain iterative reconstruction approach[END_REF], and analytic continuation (AC) [START_REF] Natterer | Image reconstruction in quantitative susceptibility mapping[END_REF]. Truncation in k-space changes the dipole kernel near the zero-cone Γ 0 , which is equivalent to adding dipole-incompatible field components. Consequently, k-space truncation causes additional streaking according to [START_REF] Zhang | Magnetic susceptibility from quantitative susceptibility mapping can differentiate new enhancing from nonenhancing multiple sclerosis lesions without gadolinium injection[END_REF], and is demonstrated in Fig. 5.

Remark 3: Observe that in Fig. 5, strong streaking appears at the interface between two different homogeneous regions (jump set). A rigorous analysis on this type of streaking can be found in [START_REF] Choi | Inverse problem in quantitative susceptibility mapping[END_REF][START_REF] Zhou | Dipole incompatibility related artifacts in quantitative susceptibility mapping[END_REF] where the notion of wave front set [START_REF] Hörmander | The analysis of Linear Partial Differential Operators I, distribution theory and Fourier Analysis[END_REF] plays a crucial role in identifying its appearance.

Bayesian inference approaches overcome the limitations of k-space truncation by directly penalizing streaking and shadow in image space, which can be solved by a variational method. In this framework, one seeks an optimal solution that minimizes

E(χ) := E data (b|χ) + λE reg (χ),
where the data fidelity term E data measures the disagreement between χ and b, and the regularization term E reg that encodes additional information measures a certain property of χ, e.g., smoothness, sparsity, and consistency with known morphologic and biologic information. The regularization parameter λ > 0 balances the two terms, which can be empirically chosen by the L-curve heuristic [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the l-curve[END_REF][START_REF]Discrete inverse problems: insight and algorithms[END_REF]. Note that this variational method is intrinsically connected to Bayesian MAP estimation [START_REF] Chan | Image processing and analysis: variational, PDE, wavelet, and stochastic methods[END_REF], i.e., p(b|χ) = e -E data (b|χ) and p(χ) = e -Ereg(χ) . If the objective E(χ) is convex, its globally optimal solution (MAP estimate) can be computed by solving the corresponding Euler-Lagrange equation-whose solutions are functions for which E(χ) is stationary. Note that writing out an explicit solution for the Euler-Lagrange equation is typically impossible, making numerical optimization required to compute a numerical solution. In general, the data fidelity term is derived from the forward problem; in most cases it is of the form ||d * χ -b|| 2 w , a weighted 2-norm. However, there are a wide range of variations in the regularization term E reg to reduce streaking and shadows, generating various Bayesian QSM algorithms as outlined in the next section.

IV. BAYESIAN QSM METHODS

In this section, we give an overview of Bayesian QSM methods based on an archetypal model. These methods are capable of reducing both streaking and shadow when proper priors are employed-which is under active development.

A. Reducing Streaking-Morphology Enabled Dipole Inversion (MEDI)

The morphology-enabled dipole inversion (MEDI) method [START_REF] Liu | Morphology enabled dipole inversion for quantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map[END_REF] is an archetypal model that has been widely used for QSM. Let χ : (Ω ⊂ R 3 ) → R, then MEDI minimizes

η Γ 1 η Γ 2 η Γ 3 χ 1 χ 2 χ 3 χ 1 -χ c χ 2 -χ c χ 3 -χ c
E MEDI (χ) := Ω |w(d * χ -b)| 2 + λ Ω |M∇χ| 1 . (17)
Here, w :

L 2 (Ω) → L 2 (Ω) is the SNR (signal-to-noise ratio)
weighting and M is a matrix acting on ∇χ. Then, the integrand in the regularization term is defined as

|M(r)∇χ(r)| 1 = i=x,y,z |(1 -δ i (r))∂ i χ(r))|.
Here, δ i (r) is an edge indicator for the i-direction at r. For example, δ x (r 0 ) = 1 if there is an edge along the xaxis at r 0 ; otherwise δ x (r 0 ) = 0. These edge indicators are derived from the magnitude image associated with the unknown susceptibility map under the assumption that edges in the two images coexist. Therefore, the weighted anisotropic total variation Ω |M∇χ| 1 penalizes only the regions where tissue structure is not expected [START_REF] Wang | Structure prior effects in bayesian approaches of quantitative susceptibility mapping[END_REF].

Remark 4: The 2-norm is also applicable for the integrand in the regularization term, which is known as weighted isotropic total variation (TV). As opposed to anisotropic TV, isotropic TV is rotation invariant. However, the difference in accuracy between the two norms in the weighted TV for QSM is very subtle; this invites further investigation.

Incorporating tissue edge information into TV is a notable feature of MEDI [START_REF] Wang | Structure prior effects in bayesian approaches of quantitative susceptibility mapping[END_REF]. Recall that streaking and tissue boundaries have similar penalties in terms of TV, but streaking only appears along the singular support of the double-cone Υ, making it distinct from the tissue edges. Furthermore, SNR weighting in the data fidelity term plays an important role in robust fitting. These have been shown in [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF][START_REF] Zhou | Dipole incompatibility related artifacts in quantitative susceptibility mapping[END_REF][START_REF] Wang | Noise effects in various quantitative susceptibility mapping methods[END_REF]. Note that the objective in [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF] does not truncate data in k-space so that more data is exploited.

B. Variations and Flexibility

Numerous QSM methods have been proposed in the context of Bayesian image reconstruction with MEDI. These methods can be understood as variations of MEDI-their common goal is to reduce streaking artifacts-and include:

• Nonlinear MEDI [START_REF] Liu | Nonlinear formulation of the magnetic field to source relationship for robust quantitative susceptibility mapping[END_REF] • Compressed sensing compensated (CSC) inversion [START_REF] Wu | Whole brain susceptibility mapping using compressed sensing[END_REF] • Homogeneity enabled incremental dipole inversion (HEIDI) [START_REF] Schweser | Quantitative susceptibility mapping for investigating subtle susceptibility variations in the human brain[END_REF] • Total generalized variation (TGV) based inversion [START_REF] Langkammer | Fast quantitative susceptibility mapping using 3d epi and total generalized variation[END_REF][START_REF] Bredies | Total generalized variation[END_REF] • Iterative LSQR (iLSQR) [START_REF] Li | A method for estimating and removing streaking artifacts in quantitative susceptibility mapping[END_REF] • Quantifying susceptibility by inversion of a perturbation (QSIP) model [START_REF] Poynton | Quantitative susceptibility mapping by inversion of a perturbation field model: correlation with brain iron in normal aging[END_REF] • Vessel orientation constrained QSM [START_REF] Cetin | Vessel orientation constrained quantitative susceptibility mapping (qsm) reconstruction[END_REF] • Structural featured based collaborative reconstruction (SFCR) [START_REF] Bao | Quantitative susceptibility mapping using structural feature based collaborative reconstruction (sfcr) in the human brain[END_REF] • MEDI based on anisotropic weighting [START_REF] Kee | Coherence enhancement in quantitative susceptibility mapping by means of anisotropic weighting in morphology enabled dipole inversion[END_REF] These methods can be combined in a plug-and-play fashion to tackle targeted applications. For instance, the nonlinear data fidelity term proposed in [START_REF] Liu | Nonlinear formulation of the magnetic field to source relationship for robust quantitative susceptibility mapping[END_REF] can be combined with TGV [START_REF] Langkammer | Fast quantitative susceptibility mapping using 3d epi and total generalized variation[END_REF][START_REF] Bredies | Total generalized variation[END_REF] to reduce noise and effects of phase unwrapping failures as Fig. 7. Comparison between magnitude image (left), QSM reconstructed using MEDI (middle) and QSM0 (right) for in vivo data. Ventricles appear more homogeneous in QSM0 with reduced shadows compared to MEDI. Here, SPURS [START_REF] Dong | Simultaneous phase unwrapping and removal of chemical shift (spurs) using graph cuts: application in quantitative susceptibility mapping[END_REF] and PDF [START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF] were used for phase unwrapping and background field removal, respectively. well as the undesirable staircasing effect [START_REF] Chan | Recent developments in total variation image restoration[END_REF], i.e., artificial piecewise constant regions, caused by TV regularization. Two regularization terms can be also combined, e.g., anisotropic weighting [START_REF] Kee | Coherence enhancement in quantitative susceptibility mapping by means of anisotropic weighting in morphology enabled dipole inversion[END_REF] and TGV to promote the parallel orientation between the magnitude and susceptibility edges, and to reduce the staircasing effect. The flexibility of the Bayesian MAP estimation framework is a unique and appealing feature in that new priors will continue to be proposed.

C. Reducing Shadow and zero reference-QSM0

MEDI and its variations listed above can effectively reduce strong streaking by TV-based regularization. However, they are inherently incapable of penalizing slowly fluctuating artifacts as demonstrated in [START_REF] Zhou | Dipole incompatibility related artifacts in quantitative susceptibility mapping[END_REF] because the gradient of such fluctuation is much smaller than that of streaking. A naïve approach is to make use of a high pass filter, however this also filters out desired susceptibility variations below the cutoff frequency. Inspired by [START_REF] De Rochefort | Quantitative susceptibility map reconstruction from mr phase data using bayesian regularization: validation and application to brain imaging[END_REF], cerebrospinal fluid (CSF) regularization to enforce a natural zero susceptibility reference needed for absolute quantification has been shown to reduce shadow artifacts [START_REF] Liu | Qsm0-qsm with automatic uniform csf zero reference[END_REF][START_REF] Liu | Qsm0: Quantitative susceptibility mapping with automatic uniform cerebrospinal fluid zero reference[END_REF], with the following cost function

E QSM0 (χ) := Ω |w(d * χ -b)| 2 + λ 1 Ω |M∇χ| 1 + λ 2 CSF |χ -χ CSF | 2 ,
where λ 1 and λ 2 are regularization parameters and χ CSF is the arithmetic mean of χ inside the CSF region. The rationale for the second regularization term is as follows: CSF is chemically almost pure water so that its susceptibility distribution is expected to be homogeneous (i.e., the variance of χ in the CSF region is very small), and the CSF region can be automatically segmented. There are often apparent susceptibility variations in the CSF region, because anisotropic white matter near the region contributes to a dipole-incompatible field. Hence, imposing CSF homogeneity in the QSM reconstruction can be an effective way of reducing shadow as demonstrated in Figs. 6 and7.

Remark 5: As its name QSM0 reveals, E QSM0 was proposed in the context of CSF zero reference [START_REF] Liu | Qsm0-qsm with automatic uniform csf zero reference[END_REF][START_REF] Liu | Qsm0: Quantitative susceptibility mapping with automatic uniform cerebrospinal fluid zero reference[END_REF]. Additionally, it has the benefit of reducing shadow that has low spatial frequency variations over the CSF region. In general, shadow artifacts may be minimized by penalizing such low spatial frequency variations where tissues are expected to be uniform.

V. COMPUTATIONAL METHODS-SOLVERS

Bayesian QSM methods need to be solved by iterative methods because writing out an explicit solution for the Euler-Lagrange equation associated with the objective is difficult or impossible. Here, we focus on minimizing E MEDI in [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF]. The major algorithmic challenge in this optimization comes from the presence of TV which is nonsmooth (albeit convex). We present three well-established algorithms for MEDI which can be also applied to the variations listed in Section IV-B: Quasi-Newton, split-Bregman, and primal-dual algorithms.

We start with a standard discrete setting for minimizing [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF]. As discussed in [107], an improper discretization scheme for the differential operator in [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF] leads to inaccuracy as well as visual artifacts in the final susceptibility map. Keeping this in mind, consider a 3D regular Cartesian grid of size N x × N y × N z as follows:

{(ih x , jh y , kh z ) : 1 ≤ i ≤ N x , 1 ≤ j ≤ N y , 1 ≤ k ≤ N z } ,
where h x , h y , and h z denote the voxel size and (i, j, k) denotes the voxel location. Then, χ is discretized in H = R Nx×Ny×Nz where a scalar product is defined as

χ, ψ H = i,j,k χ i,j,k ψ i,j,k , χ, ψ ∈ H.
For χ ∈ H, the discrete gradient is defined by (∇χ) i,j,k = ((∇χ) x i,j,k , (∇χ) y i,j,k , (∇χ) z i,j,k ), where

(∇χ) x i,j,k = (χ i+1,j,k -χ i,j,k )/h x if i < N x , 0 if i = N x ; (∇χ) y i,j,k = (χ i,j+1,k -χ i,j,k )/h y if j < N y , 0 if j = N y ; (∇χ) z i,j,k = (χ i,j,k+1 -χ i,j,k )/h z if k < N z , 0 if k = N z .
Notice that the forward difference scheme with the Neumann boundary condition is used for the discrete gradient [START_REF] Maerz | Image quality improvement using short range finite difference in qsm reconstruction[END_REF][START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF]. Also, the gradient vector ∇χ is a vector in K = H × H × H where a scalar product is given as

ξ, ζ K = i,j,k ξ x i,j,k ζ x i,j,k + ξ y i,j,k ζ y i,j,k + ξ z i,j,k ζ z i,j,k , ξ, ζ ∈ K.
Now, we define a linear map M : K → K such that

ξ i,j,k → (Mξ) i,j,k = (M x i,j,k ξ x i,j,k , M y i,j,k ξ y i,j,k , M z i,j,k ξ z i,j,k ).
The objective in ( 17) is discretized accordingly, yielding the following optimization problem in a finite dimensional vector space:

χ ∈ argmin χ ||w(d * χ -b)|| 2 2 + λ||M∇χ|| 1,1 , (18)
where w : H → H is SNR weighting, d is the discretized dipole kernel, and b ∈ H is the RDF. We define || • || 1,1 as

||ξ|| 1,1 = i,j,k |ξ x i,j,k | + |ξ y i,j,k | + |ξ z i,j,k |.
After vectorization, the formulation ( 18) can be rewritten as:

χ ∈ argmin χ ||W (F H DF χ -b)|| 2 2 + λ||M Gχ|| 1 ( 19
)
where W is a diagonal matrix associated with w, F is the discrete Fourier transform (DFT) matrix, D is a diagonal matrix whose diagonal elements correspond to the Fourier transform of d, and the superscript H denotes the Hermitian transpose in the data fidelity term. The finite difference matrix G : R NxNyNz → R 3NxNyNz is constructed from the forward difference scheme with the Neumann boundary condition, and M is a diagonal matrix associated with M.

A. Quasi-Newton Method (Gauss-Newton in QSM)

Since the optimization problem ( 19) is convex, a necessary and sufficient condition for optimality of ( 19) is to find

χ such that 0 ∈ ∂ χ E( χ), where E(χ) := ||W (F H DF χ -b)|| 2 2 + λ||M Gχ|| 1 . Here, ∂ χ E(χ) is given as ∂ χ E(χ) = 2(W F H DF ) H (W F H DF )χ -(W F H DF ) H b + λ(M G) H M Gχ diag(M Gχ) 2 . ( 20
)
Then, the lagged diffusivity fixed-point iteration (LDFPI) method [START_REF] Vogel | Iterative methods for total variation denoising[END_REF][START_REF] Chan | On the convergence of the lagged diffusivity fixed point method in total variation image restoration[END_REF] is used to linearize the nonlinear equation ∂ χ E(χ) = 0 by substituting 1/ diag(M Gχ) 2 with 1/ diag(M Gχ k ) 2 + , where χ k denotes the k-th iterate and is a strictly positive relaxation parameter which avoids division by zero. Note that ∂ χ ||M Gχ|| 1 leads to a nonlinear diffusion equation, and this linearization procedure makes the diffusivity coefficient lag one step (iteration) behind. LDFPI then computes χ k+1 by solving the following linear system

L(χ k )χ k+1 = b, (21)
where

L(χ k ) := 2(W F H DF ) H (W F H DF )χ + λ(M G) H (M G)/ diag(M Gχ k ) 2 + , and b := (W F H DF ) H b. Once χ k+1 is computed, χ k ← χ k+1
and solve [START_REF] Azuma | Lateral asymmetry and spatial difference of iron deposition in the substantia nigra of patients with parkinson disease measured with quantitative susceptibility mapping[END_REF] until ||χ k+1 -χ k || becomes sufficiently small (fixed-point iteration). The rate of convergence largely depends on how fast LDFPI solves a series of the linear systems [START_REF] Azuma | Lateral asymmetry and spatial difference of iron deposition in the substantia nigra of patients with parkinson disease measured with quantitative susceptibility mapping[END_REF]. Krylov subspace methods such as conjugate gradient (CG) can be used with preconditioning techniques [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] From ( 20) and ( 21), observe that b = L(χ k )χ k -∂ χ E(χ k ). Using this relation, we can rewrite [START_REF] Azuma | Lateral asymmetry and spatial difference of iron deposition in the substantia nigra of patients with parkinson disease measured with quantitative susceptibility mapping[END_REF] as follows: where δχ = χ k+1 -χ k . Then, solving [START_REF] Barbosa | Quantifying brain iron deposition in patients with parkinson's disease using quantitative susceptibility mapping, r2 and r2[END_REF] with respect to δχ is a quasi-Newton method where L(χ k ) is an estimate of the Hessian matrix of E(χ) at χ k . To see this, consider the second order Taylor expansion of E(χ) around χ k as follows:

L(χ k )δχ = -∂ χ E(χ k ), (22)
E(χ k + δχ) E(χ k ) + ∂ χ E(χ k ) δχ + 1 2 δχ [HE(χ k )]δχ,
where HE(χ k ) is the Hessian matrix. Newton's method sets the derivative of this expansion with respect to δχ equal to zero:

HE(χ k )δχ = -∂ χ E(χ k ),
where

HE(χ k ) (23) = L(χ k ) + λ(M G) H (M Gχ k )∂ χ 1 diag(M Gχ) 2 + χ=χ k . L(χ k
) is an estimate of the Hessian HE(χ k ) by dropping the last term from [START_REF] Du | Quantitative susceptibility mapping of the midbrain in parkinson's disease[END_REF]. The quasi-Newton approach, therefore, iteratively approximates the original nonlinear minimization problem ( 19) by a series of quadratic minimization problems. An advantage of the quasi-Newton form ( 22) over the fixedpoint form [START_REF] Azuma | Lateral asymmetry and spatial difference of iron deposition in the substantia nigra of patients with parkinson disease measured with quantitative susceptibility mapping[END_REF] is that [START_REF] Barbosa | Quantifying brain iron deposition in patients with parkinson's disease using quantitative susceptibility mapping, r2 and r2[END_REF] tends to be less sensitive to roundoff error [START_REF] Vogel | Computational methods for inverse problems[END_REF].

The accuracy of the quasi-Newton form [START_REF] Barbosa | Quantifying brain iron deposition in patients with parkinson's disease using quantitative susceptibility mapping, r2 and r2[END_REF] depends on the value of the numerical relaxation parameter , which, by default, is set to 10 -6 for QSM. Accuracy can be improved by decreasing but this comes at the cost of significantly slowing down the rate of convergence [107]. This can be overcome by using the so-called continuation method. That is, at each outer iteration, we reset by the rule ← γ • where γ is a decreasing factor. We initialize = 0.1 and set γ to values ranging from 0 to 1. Fig. 8 shows relative errors when γ is 0.01, 0.05, 0.1, 0.3, 0.5, and 0.8.

Remark 6: This quasi-Newton method has been frequently referred to as Gauss-Newton in the QSM literature [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF]. This may be because Newton's method has been often referred to as the Newton-Raphson method for finding roots in physics.

B. Split-Bregman Method

The split-Bregman method [START_REF] Goldstein | The split bregman method for l1regularized problems[END_REF] makes use of variable splitting and the Bregman iteration. Unlike the quasi-Newton method, this solves (19) without numerical approximation for TV regularization. The method of variable splitting reformulates [START_REF] Hwang | Texture analyses of quantitative susceptibility maps to differentiate alzheimer's disease from cognitive normal and mild cognitive impairment[END_REF] as

( χ, v) (Constrained) = argmin (χ,v) W (F H DF χ -b) 2 2 + λ v 1 s.t. v = M Gχ.
A quadratic penalty is incorporated and we have

( χ, v) (Unconstrained) = argmin (χ,v) W (F H DF χ -b) 2 2 + λ v 1 + ρ 2 v -M Gχ 2 2 .
This can be solved by alternating continuation, i.e., taking ρ → ∞ one alternatively minimizes the objective with respect to χ and v. However, the problem becomes ill-conditioned as ρ → ∞ [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Instead of taking ρ → ∞, the Bregman iteration [START_REF] Yin | Bregman iterative algorithms for \ell 1-minimization with applications to compressed sensing[END_REF] is used with a fixed value of ρ. This procedure makes its numerical behavior much more stable; the iterative sequence is as follows:

χ k+1 (24)
= argmin χ W (F H DF χ -b) 2 2 + ρ 2 v k -M Gχ -z k 2 2 = argmin χ W F H DF ρ/2 M G χ - W b ρ/2(v k -z k ) 2 2 , v k+1 (25)
= argmin v λ||v|| 1 + ρ 2 ||v -M Gχ k+1 -z k || 2 2 = prox ρ/λ,||•||1 (M Gχ k+1 + z k ), z k+1 (26)
= z k + (M Gχ k+1 -v k+1 ),
where Krylov subspace methods, e.g., CG or LSQR, can be used for [START_REF] Guan | Regionally progressive accumulation of iron in parkinson's disease as measured by quantitative susceptibility mapping[END_REF]. The proximal operator prox ρ/λ,||•||1 (•) in ( 25) is given as soft-thresholding [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. The last line of the algorithm ( 26) is known as the Bregman update and its derivation can be found in [START_REF] Yin | Bregman iterative algorithms for \ell 1-minimization with applications to compressed sensing[END_REF]. In QSM, the split-Bregman method was applied to [START_REF] Hwang | Texture analyses of quantitative susceptibility maps to differentiate alzheimer's disease from cognitive normal and mild cognitive impairment[END_REF] without SNR weighting in [START_REF] Bilgic | Fast quantitative susceptibility mapping with l1-regularization and automatic parameter selection[END_REF] to accelerate the speed of deconvolution. The method was subsequently considered with SNR weighting in [START_REF] Chung | Computational methods for image reconstruction[END_REF].

Remark 7: Note that the split-Bregman method with linear constraints is equivalent to the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Hence, its convergence can be analyzed by ADMM.

C. Primal-Dual Algorithm

Similar to the distributional derivative for a function whose differentiation is not pointwise defined [START_REF] Stein | Functional Analysis: Introduction to Further Topics in Analysis[END_REF][START_REF] Blanchard | Mathematical Methods in Physics: Distributions, Hilbert Space Operators, Variational Methods, and Applications in Quantum Physics[END_REF], we can derive a generalized (dual) definition of the weighted (anisotropic) TV [107, 120] as follows:

Ω |M∇χ| 1 = sup ξ:||ξ||∞≤1 Ω M∇χ, ξ , (27)
where the dual variable ξ : Ω → R 3 is a smooth map and || • || ∞ is the pointwise maximum norm. Replacing the TV term in E MEDI with its dual formulation [START_REF] Ide | Internal structures of the globus pallidus in patients with parkinson's disease: evaluation with quantitative susceptibility mapping (qsm)[END_REF], we notice that minimizing E MEDI in χ is equivalent to finding a saddle-point ( χ, ξ) such that

min χ max ξ:||ξ||∞≤1 Ω |w(d * χ -b)| 2 + λ Ω M∇χ, ξ .
Then, one can use gradient descent in χ and gradient ascent in ξ with the relation M∇χ, ξ = χ, (M∇) H ξ . The data term can be further dualized to avoid solving a linear system that appears as a nested routine as seen in the quasi-Newton and split-Bregman methods. For technical and implementation details, we refer the reader to [107].

D. Final Remarks on Solvers

As opposed to Newton's method in optimization that makes use of the Hessian matrix of the objective function, first-order methods (quasi-Newton, split-Bregman, and primal-dual) have a linear rate of convergence for the complete class of nonsmooth convex problems [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. Since only subtle differences in speed and accuracy are observed, any method can be used, for instance NESTA [START_REF] Becker | Nesta: A fast and accurate firstorder method for sparse recovery[END_REF] can be also used. That being said, it is unclear if second-order methods (Newton-type methods) would outperform the first-order methods because they require computation of the Hessian matrix at each step; this is not only expensive, but also noise sensitive [START_REF] Jaggi | Sparse convex optimization methods for machine learning[END_REF].

New regularization terms may appear nonconvex, therefore the first-order methods may not be useful in finding a nearoptimal solution. In this case, the Hessian matrix can help escape saddle points towards a globally optimal solution [START_REF] Dauphin | Identifying and attacking the saddle point problem in high-dimensional non-convex optimization[END_REF]. Algorithmic advances will be then necessary to catch up with realistic nonconvex energies.

VI. COMPUTATIONAL METHODS-PRECONDITIONING

Numerical computation of a solution in Bayesian QSM can be accelerated using preconditioning techniques that incorporate prior knowledge of the targeted susceptibility solution. As known in matrix inversion [START_REF] Golub | Matrix computations[END_REF][START_REF] Trefethen | Numerical linear algebra[END_REF], a well-chosen left/right preconditioner can accelerate the rate of convergence of an iterative solver. In Bayesian MAP estimation, preconditioning may play a role in guiding the path of a minimizing sequence in such a way that it emphasizes certain properties of a final solution [START_REF] Calvetti | Bayes meets krylov: preconditioning cgls for underdetermined systems[END_REF][START_REF] Arridge | A priorconditioned lsqr algorithm for linear ill-posed problems with edge-preserving regularization[END_REF]. Here, we present an overview of recent advances in background field removal methods in the context of preconditioning.

A. Background Field Removal

The field generated by the sources inside a tissue volume of interest (VOI) such as the brain, known as tissue (local) field, turns out to be much smaller than the field generated by the sources outside the VOI, which include strong susceptibility gradients at air-tissue interfaces and in the skull, known as background field [START_REF] Zhou | Background field removal by solving the laplacian boundary value problem[END_REF]. A typical QSM strategy is to first perform data fitting and removal of this background field on the total field estimated from MRI phase, which is referred to as background field removal [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF][START_REF] De Rochefort | Quantitative susceptibility map reconstruction from mr phase data using bayesian regularization: validation and application to brain imaging[END_REF][START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF][START_REF] Schweser | An illustrated comparison of processing methods for phase mri and qsm: removal of background field contributions from sources outside the region of interest[END_REF]. The tissue field after background field removal is then inverted to generate a tissue susceptibility map. Because background field removal does not require precise susceptibility source localization, it can be executed rapidly. With the removal of the background field, the susceptibility range to be searched for mapping tissue susceptibility can be reduced by ∼ 10 fold (but this breaks down in case of hemorrhage [START_REF] Liu | Nonlinear formulation of the magnetic field to source relationship for robust quantitative susceptibility mapping[END_REF][START_REF] Liu | Preconditioned total field inversion (tfi) method for quantitative susceptibility mapping[END_REF], leading to a rapid convergence in solving [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF]. Therefore, current QSM algorithms typically consist of two steps [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF]: background field removal and tissue field inversion.

As discussed in [START_REF] Wang | Quantitative susceptibility mapping (qsm): decoding mri data for a tissue magnetic biomarker[END_REF], the lack of MRI signal in the background is one of the QSM challenges to determine the background field inside the VOI accurately. All current methods to remove the background field in the tissue VOI are based on Maxwell's equations, which state that the background field in the tissue VOI obeys the homogeneous Laplacian equation. The lack of boundary condition makes it an ill-posed problem to solve this Laplacian equation [START_REF] Zhou | Background field removal by solving the laplacian boundary value problem[END_REF]. Numerous formulations of prior knowledge have been proposed for estimating background field, including:

• High pass filtering (HPF) [START_REF] Haacke | Susceptibility weighted imaging (swi)[END_REF][START_REF] Langham | Retrospective correction for induced magnetic field inhomogeneity in measurements of large-vessel hemoglobin oxygen saturation by mr susceptometry[END_REF][START_REF] Wang | Artery and vein separation using susceptibility-dependent phase in contrast-enhanced mra[END_REF] • Projection onto dipole fields (PDF) [START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF] • Sophisticated harmonic artifact reduction for phase data (SHARP) [START_REF] Schweser | Differentiation between diamagnetic and paramagnetic cerebral lesions based on magnetic susceptibility mapping[END_REF] • Laplacian boundary value (LBV) method [START_REF] Zhou | Background field removal by solving the laplacian boundary value problem[END_REF] • Iterative spherical mean value method (iSMV) [START_REF] Wen | An iterative spherical mean value method for background field removal in mri[END_REF] • Regularization-enabled SHARP (RESHARP) [START_REF] Sun | Background field removal using spherical mean value filtering and tikhonov regularization[END_REF] For a complete review on background field removal methods, we refer the reader to [START_REF] Schweser | An illustrated comparison of processing methods for phase mri and qsm: removal of background field contributions from sources outside the region of interest[END_REF]. It should be noted that these prior knowledge formulations contain error, and the background field removal assumption of reducing the susceptibility search range fails when susceptibility values are large, such as in the case of hemorrhage [START_REF] Liu | Preconditioned total field inversion (tfi) method for quantitative susceptibility mapping[END_REF]. It was recently shown that a tissue susceptibility map could be recovered without background field removal but this required a long processing time [START_REF] Liu | Preconditioned total field inversion (tfi) method for quantitative susceptibility mapping[END_REF], because a tissue susceptibility map was still computed from the total field by minimizing E MEDI in [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF]. This work elucidates the rate of convergence of an iterative solver, demonstrating that the background field removal step plays a role in preconditioning in a classical sense. It also demonstrates that the error in a background field removal method can be avoided by focusing on minimizing E MEDI from the total field.

The remaining concern for minimizing E MEDI is systematically incorporating preconditioning in such a way that its solution is computed as fast as the current two-step procedure (background field removal + tissue field inversion) in Bayesian QSM. In the following sections, we present two notable examples.

B. Preconditioned MEDI

Preconditioning can be systematically incorporated into a Bayesian QSM method [START_REF] Liu | Preconditioned total field inversion (tfi) method for quantitative susceptibility mapping[END_REF] by modifying the cost function Notice that preconditioning dramatically reduces the hypointense artifact around the hemorrhage site. For both MEDI and TFI, SPURS [START_REF] Dong | Simultaneous phase unwrapping and removal of chemical shift (spurs) using graph cuts: application in quantitative susceptibility mapping[END_REF] was used for phase unwrapping, and PDF [START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF] was additionally used for background field removal for MEDI.

E MEDI in (17) with a preconditioner P :

E PMEDI (P y) := Ω |w(d * P y -b)| 2 + λ Ω |M∇P y| 1 .
which is minimized with respect to y to generate a susceptibility solution χ = P y. Here, b is the total field, and the preconditioner P is chosen from prior knowledge about the distribution of brain susceptibility. In the case of QSM reconstruction of a patient with intracerebral hemorrhage, P is chosen as a binary matrix with higher weights in the hemorrhage than other regions to reflect the strong susceptibility contrast between the hemorrhage and the surrounding brain tissue [START_REF] Liu | Preconditioned total field inversion (tfi) method for quantitative susceptibility mapping[END_REF]. Both simulation and in vivo results (Fig. 9) show that the preconditioner helps suppress shadow artifacts adjacent to the hemorrhage site by improving the convergence behavior of CG involved in the inner loop of the quasi-Newton method presented in Section V-A.

C. Differential Model

Along with the integral approach in [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF] where the data term is derived from (8), differential approaches proposed in [136- (150 CG iterations). MEDI + SMV produces more homogeneous map at the cost of brain erosion. Here, SPURS [START_REF] Dong | Simultaneous phase unwrapping and removal of chemical shift (spurs) using graph cuts: application in quantitative susceptibility mapping[END_REF] and PDF [START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF] were used for phase unwrapping and background field removal, respectively. 139] are also commonly used Bayesian QSM methods which can be understood in the context of preconditioning. These approaches make use of the PDE (4) or [START_REF] Deistung | Overview of quantitative susceptibility mapping[END_REF] as the data fidelity term of E MEDI in [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF]. As a result, a background field removal step is no longer needed, making these approaches as fast as the two-step (e.g., PDF + MEDI) procedure. Additionally, it has been shown that some of the approaches known as singlestep QSM [START_REF] Bilgic | Single-step qsm with fast reconstruction[END_REF]139] have less shadow artifacts. Recall that the integral form (8) is derived from the differential form (4) with the assumption that the Fourier transform of each term in (4) exists and k = 0. Therefore, both integral and differential representations appear to be equivalent. Yet, a rigorous analysis needs to be carried out because the differential approach seems more capable of reducing shadow artifacts than the integral approach.

From a computational point of view, this behavior may be explained by the 4π|k| 2 left preconditioner that comes from the Laplacian operator. In other words, as opposed to the data term of the integral approach in (17), the following differential term

Ω | w∆(d * χ -b)| 2 = Ω | wF -1 [4π|k| 2 (DX -B)]| 2
may emphasize the high frequency components of a candidate susceptibility distribution more than its low frequency components containing slowly varying dark shadows. Note that unlike the single-step QSM model [START_REF] Bilgic | Single-step qsm with fast reconstruction[END_REF]139], here we have considered the SNR weighting w which is modified from w in the integral approach (17) as follows: Let L be the Laplacian operator in a matrix form, and let n be the random noise vector n = d * χ -b. Since we define w as the inverse of the square root of the covariance matrix of the noise vector n, a proper way of determining the new SNR weighting w is to set it to the inverse of the square root of the covariance matrix associated with Ln.

The differential approach can be implemented by using the spherical mean value (SMV) kernel and incorporated into E MEDI ; this is denoted by MEDI + SMV. Here, we use the diagonal components of the inverse of the square root of the covariance matrix associated with Ln for w. Fig. 10 shows both simulation and in vivo results, where the tissue field contains unresolved background fields and is used as an input for both MEDI and MEDI + SMV. This SMV kernel considerably improves the rate of convergence of CG. Visually, MEDI + SMV produces a more homogeneous map than MEDI. However, any SMV implementation causes erosion of the VOI region since the field outside the VOI cannot be measured. This causes loss of anatomical information at the boundaries of the brain, especially when the voxel size is large.

Remark 8: These Bayesian preconditioning techniques have only been empirically explored. A rigorous analysis is required as preconditioning seems to give promising results in QSM.

VII. DISCUSSION AND CONCLUSION

QSM based on the Bayesian MAP estimation framework can accurately and robustly map strong isotropic sources, such as biometals including iron in ferritin and deoxy-heme in hemoglobin [START_REF] Wang | Clinical quantitative susceptibility mapping (qsm): Biometal imaging and its emerging roles in patient care[END_REF][START_REF] Langkammer | Quantitative susceptibility mapping (qsm) as a means to measure brain iron? a post mortem validation study[END_REF], calcium in bones and other pathologic calcifications [START_REF] Chen | Intracranial calcifications and hemorrhages: characterization with quantitative susceptibility mapping[END_REF][START_REF] Dimov | Bone quantitative susceptibility mapping using a chemical speciesspecific r2* signal model with ultrashort and conventional echo data[END_REF], and highly paramagnetic gadolinium and iron in contrast agents [START_REF] Liu | Unambiguous identification of superparamagnetic iron oxide particles through quantitative susceptibility mapping of the nonlinear response to magnetic fields[END_REF][START_REF] Wong | Visualizing and quantifying acute inflammation using icam-1 specific nanoparticles and mri quantitative susceptibility mapping[END_REF]. The accuracy and robustness of Bayesian MAP estimation are well established, such as in image processing known to outperform classical multistep algorithms for tasks including denoising [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], deblurring [START_REF] Chan | Total variation blind deconvolution[END_REF], segmentation [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF], and optical flow estimation [START_REF] Horn | Determining optical flow[END_REF]. Effective QSM applications in clinical problems that are often complex require further technical innovation. Although we have used the convex cost in [START_REF] Eskreis-Winkler | Multiple sclerosis lesion geometry in quantitative susceptibility mapping (qsm) and phase imaging[END_REF] for elucidating the Bayesian QSM framework, more sophisticated priors capable of precisely modeling the underlying susceptibility distribution are desired so that streaking and shadow can be better penalized. Fast and robust algorithms that can handle nonconvex priors are often necessary as such priors may well be no longer convex, such as when fat presents in tissue [START_REF] Dimov | Bone quantitative susceptibility mapping using a chemical speciesspecific r2* signal model with ultrashort and conventional echo data[END_REF][START_REF] Dong | Simultaneous phase unwrapping and removal of chemical shift (spurs) using graph cuts: application in quantitative susceptibility mapping[END_REF][START_REF] Dimov | Joint estimation of chemical shift and quantitative susceptibility mapping (chemical qsm)[END_REF][START_REF] Sharma | Quantitative susceptibility mapping in the abdomen as an imaging biomarker of hepatic iron overload[END_REF] or when SNR is poor [START_REF] Liu | Nonlinear formulation of the magnetic field to source relationship for robust quantitative susceptibility mapping[END_REF]. As seen in Section VI, Bayesian methods benefit from further investigation in preconditioning techniques. These innovations may be made based on a deeper understanding of the biophysics behind tissue magnetic susceptibility and solvers in numerical optimization.

A particularly complex problem is the molecular microstructure and its anisotropic susceptibility of white matter (WM) in the brain [START_REF] Wisnieff | The influence of molecular order and microstructure on the r2* and the magnetic susceptibility tensor[END_REF]. This WM susceptibility anisotropy has been attributed to anisotropic lipid molecules that are radially oriented in the myelin sheathes of white matter fibers [START_REF] Li | Mapping magnetic susceptibility anisotropies of white matter in vivo in the human brain at 7t[END_REF][START_REF] Li | Magnetic susceptibility anisotropy of human brain in vivo and its molecular underpinnings[END_REF][START_REF] Wisnieff | Magnetic susceptibility anisotropy: cylindrical symmetry from macroscopically ordered anisotropic molecules and accuracy of mri measurements using few orientations[END_REF][START_REF] Lounila | Effects of orientational order and particle size on the nmr line positions of lipoproteins[END_REF]. Accordingly, a hollow cylinder fiber model may account for effects from both microstructure and anisotropy of WM [START_REF] Wharton | Fiber orientation-dependent white matter contrast in gradient echo mri[END_REF], which is used to study the orientation of WM fibers [START_REF]Gradient echo based fiber orientation mapping using r2* and frequency difference measurements[END_REF][START_REF]Effects of white matter microstructure on phase and susceptibility maps[END_REF][START_REF] Aggarwal | B0-orientation dependent magnetic susceptibility-induced white matter contrast in the human brainstem at 11.7 t[END_REF]. WM microstructure can be inspected from the orientation dependent local frequency shift [START_REF] He | Biophysical mechanisms of phase contrast in gradient echo mri[END_REF][START_REF] Luo | Magnetic susceptibility induced white matter mr signal frequency shiftsexperimental comparison between lorentzian sphere and generalized lorentzian approaches[END_REF]. The molecular anisotropy of WM in a non-local frequency perturbation [START_REF] Lee | Sensitivity of mri resonance frequency to the orientation of brain tissue microstructure[END_REF] may be estimated using various inversion algorithms [START_REF] Liu | Susceptibility tensor imaging[END_REF][START_REF] Li | Mapping magnetic susceptibility anisotropies of white matter in vivo in the human brain at 7t[END_REF][START_REF] Li | Magnetic susceptibility anisotropy of human brain in vivo and its molecular underpinnings[END_REF][START_REF] Wisnieff | Magnetic susceptibility anisotropy: cylindrical symmetry from macroscopically ordered anisotropic molecules and accuracy of mri measurements using few orientations[END_REF]. A recent and more complex model combines the generalized Lorentzian approach [START_REF] He | Biophysical mechanisms of phase contrast in gradient echo mri[END_REF] with the lipoprotein structure of myelin sheath [START_REF] Sukstanskii | On the role of neuronal magnetic susceptibility and structure symmetry on gradient echo mr signal formation[END_REF] and water attempting to explain experimental phase data [START_REF] Yablonskiy | Generalized lorentzian tensor approach (glta) as a biophysical background for quantitative susceptibility mapping[END_REF]. Although QSM can be used to study WM fibers using multiorientation scans, a reliable measurement of anisotropy in WM and microstructure in a voxel from a single-angle acquisition has not been achieved. In fact, clinical MR scanners do not allow "perfect" multi-orientation sampling using a combination of uniformly distributed field directions because of the limited range of orientations with human subjects [START_REF] Wisnieff | Magnetic susceptibility anisotropy: cylindrical symmetry from macroscopically ordered anisotropic molecules and accuracy of mri measurements using few orientations[END_REF].

In summary, we have tried to provide a systematic overview of the current status and recent advances in QSM methods, and remaining challenges. The notion of a fundamental solution in the theory of linear PDEs makes it possible to write out an explicit form of the solution for the magnetic field to susceptibility source inversion problem. We have demonstrated that streaking and shadows are two types of artifacts characterized by the PDE solution for field data deviating from or incompatible with the continuous dipole field pattern. As dipole-incompatibility in field data is unavoidable in practice, it is critical to devise effective methods capable of reducing streaking and shadow artifacts. Early attempts that make use of truncation in k-space unnecessarily cause dipoleincompatibility in the field, giving rise to additional streaking and shadows. On the other hand, Bayesian QSM methods provide desired results with significant improvements over the truncation-based methods, allowing for the design of problemspecific and physics-based priors. The QSM community continues to devise new priors that will improve the current state of the art. As these new priors may be nonconvex, they will be more capable of modeling the complicated biophysical nature behind QSM; therefore, the community looks forward to algorithmic advances with preconditioning techniques.

Fig. 1 .

 1 Fig. 1. The fundamental solution g(r) is displayed on the left with the window level [0, 5] in sagittal view. Notice that the conic surface is spanned by ∼ 35.3 • . Its Fourier transform F [g](k) = 1/P (k) = 1/(4π|k| 2 D(k)) is displayed in the middle with the window level [-100, 100] in sagittal view. The dipole kernel D(k) is displayed on the right with the window level [-0.01, 0.01] in sagittal view. Here, the conic surface is spanned by ∼ 54.7 • .Observe that the conic surface of the fundamental solution is perpendicular to that of D(k).

Reference χ c b c Recovered χ c Fig. 2 .

 c2 Fig. 2. The dipole compatible field bc displayed in the middle with the window level [-0.3, 1] was generated from the dipole model (13) using the reference χc shown on the left with [-0.5 ppm, 1 ppm]. We recovered χc shown on the right with [-0.5 ppm, 1 ppm] using (14) and[START_REF] Zhang | Longitudinal change in magnetic susceptibility of new enhanced multiple sclerosis (ms) lesions measured on serial quantitative susceptibility mapping (qsm)[END_REF]. Here, the largest sphere is 0.7 ppm and the background (outside the sphere) is 0 ppm. The four spheres inside the largest one, have susceptibility values of 1, 0.85, 0.45, and 0.2 ppm, going counterclockwise starting from the bottom.

Fig. 3 .

 3 Fig. 3. The top row shows a dipole-incompatible field b i (point sources) with the window level [-0.3, 1]; the susceptibility values of white and black dots are +1 and -1, respectively. The rightmost image in the bottom row is the Laplacian of b i in sagittal with the window level [-0.5, 1]. The middle row shows the recovered susceptibility χ i with the window level [-0.5 ppm, 1 ppm]. Artifacts manifest as ringing in the axial view and streaking in sagittal and coronal views. The leftmost and middle panels in the bottom row show the recovered map with the window level [-0.5 ppm, 1 ppm] fromthe RDF bc+b i where bc is given as in Fig.2. Susceptibility was reconstructed using TKD[START_REF] Shmueli | Magnetic susceptibility mapping of brain tissue in vivo using mri phase data[END_REF] with 0.005 as the value of truncation.

Fig. 4 .

 4 Fig. 4. The top row shows a dipole-incompatible field b i (smooth distribution) with the window level [-0.3, 1]; the values of susceptibility distribution that generates the field b i are all negative. The rightmost image in the bottom row is the Laplacian of b i with the window level [-0.5, 1]. The middle row shows the recovered susceptibility χ i with the window level [-0.5 ppm, 1 ppm]. Artifacts manifest as dark shadows. The leftmost and middle panels in the bottom row show the recovered map with the window level [-0.5 ppm, 1 ppm] from the RDF bc + b i where bc is given as in Fig. 2. Susceptibility was reconstructed using TKD [58] with 0.005 as the value of truncation.

Fig. 5 .

 5 Fig. 5. The top row shows the indicator function of Γ denoted by η Γ , where Γ is defined as {k ∈ R 3 : |D(k)| < } with different values of , where 1 = 0.1, 2 = 0.01, and 3 = 0.001. That is, η Γ (k) = 1 for k ∈ Γ and η Γ (k) = 0 for k ∈ R 3 \ Γ , i.e., white = 1 and black = 0. For the dipole-compatible field described in Fig. 2, the second row shows the corresponding TKD solutions defined as χ = F -1 (Bc(k)/D(k)(1η Γ (k))) = χc -F -1 (Bc(k)/D(k)) * F -1 (η Γ (k))) with the window level [-0.5 ppm, 1 ppm]. The bottom row shows the difference image between χc in Fig. 2 and χ with the window level [-0.5 ppm, 1 ppm]. Notice that even for a dipole-compatible field, truncation causes streaking artifacts around the high-frequency regions (interfaces between regions-jump set).

Fig. 6 .

 6 Fig. 6. Simulation result. Top row: Reference QSM (left) which does not contain white matter anisotropy, MEDI (middle) and QSM0 (right) generated from a field which contains white matter anisotropy. Bottom row: Highlighted QSM for CSF. Notice the reduced shadow artifact indicated by the red arrows in QSM0 compared to MEDI. Meanwhile, the homogeneity of CSF is improved using QSM0.

Fig. 8 .

 8 Fig.8. The continuation method helps speed up the quasi-Newton formulation of LDFPI (qNewton) if the decreasing factor γ is properly chosen (see the purple, green, and red curves below the black). Continuation converges rather slowly when γ ≥ 0.3; see the yellow, blue, and cyan curves above the black.

Fig. 9 .

 9 Fig. 9. MEDI versus preconditioned MEDI (TFI). (a) Plot of error between estimated and true brain susceptibility with respect to the CG iteration number in numerical simulation (b) QSM reconstructed by MEDI and TFI in numerical simulation and in vivo.Notice that preconditioning dramatically reduces the hypointense artifact around the hemorrhage site. For both MEDI and TFI, SPURS[START_REF] Dong | Simultaneous phase unwrapping and removal of chemical shift (spurs) using graph cuts: application in quantitative susceptibility mapping[END_REF] was used for phase unwrapping, and PDF[START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF] was additionally used for background field removal for MEDI.

Fig. 10 .

 10 Fig. 10. Comparison between MEDI and MEDI + SMV. (a) Plot of error between estimated and true brain susceptibility with respect to the CG iteration number in numerical simulation. (b) QSM reconstructed by MEDI and MEDI + SMV in numerical simulation (100 CG iterations) and in vivo(150 CG iterations). MEDI + SMV produces more homogeneous map at the cost of brain erosion. Here, SPURS[START_REF] Dong | Simultaneous phase unwrapping and removal of chemical shift (spurs) using graph cuts: application in quantitative susceptibility mapping[END_REF] and PDF[START_REF] Liu | A novel background field removal method for mri using projection onto dipole fields[END_REF] were used for phase unwrapping and background field removal, respectively.

Note that the acronym RDF has been often referred to as tissue (local) field (:= total fieldbackground field) in the literature. In this paper, however, RDF means total field by its definition.

ACKNOWLEDGMENT

This work was supported in part by grants from National Institute of Health (R01 NS095562, R01 NS090464, R01 NS072370, S10 OD021782). Jin Keun Seo was supported by NRF grant 2015R1A5A1009350.


Andreas Herzig
email: herzig@irit.fr

Jérôme Lang
email: lang@irit.fr

Dominique Longin
email: longin@irit.fr

I thought you didn't know! On belief revision in dynamic doxastic logic

Introduction

Suppose agent i believes that agent j does not believe that p. What is i's belief state like after j informs i that p?

Scenarios of this kind can be expressed nicely in dynamic doxastic logics. There, to every action α is associated a dynamic operator [α], to every agent i is associated a belief operator B i , and to every set of agents {i, j, . . .} is associated a common belief operator B {i,j,...} . Speech acts can be viewed as particular actions, and we can write j's action of informing i that p as inform(j, i, p).

To simplify our exposition let us make the hypothesis that every agent is sincere: he cannot assert things he doesn't believe himself. This can be expressed by the axiom

¬B j ϕ → [inform(j, i, ϕ)]⊥
Then standard KD45 principles entail that B i ¬B j ϕ → B i [inform(j, i, ϕ)]⊥, i.e. in our scenario i believes j cannot tell him that p.

We moreover make the hypothesis that sincerity is the only precondition of inform-actions. We therefore have the axiom

B j ϕ → inform(j, i, ϕ)
Such scenarios, where an agent learns that some action α he wrongly believed to be inexecutable has nevertheless occurred, pose problems for the logics of actions that have been proposed in the literature to account for the evolution of knowledge and belief. Some of them [START_REF] Moore | A formal theory of knowledge and action[END_REF][START_REF] Scherl | The frame problem and knowledge producing actions[END_REF][START_REF] Shapiro | Specifying communicative multi-agent systems[END_REF][START_REF] Thielscher | Representing the knowledge of a robot[END_REF][START_REF] Shapiro | Iterated belief change in the situation calculus[END_REF][START_REF] Herzig | A logic for planning under partial observability[END_REF] deliberately restrict their attention to true belief (alias knowledge) in order to avoid such scenarios. Indeed, for the other approaches in [START_REF] Scherl | The frame problem and knowledge producing actions[END_REF][START_REF] Gerbrandy | Reasoning about information change[END_REF][START_REF] Gerbrandy | Bisimulations on Planet Kripke[END_REF][START_REF] Baltag | The logic of public announcements, common knowledge, and private suspicions[END_REF][START_REF] Baltag | A logic of epistemic actions[END_REF][START_REF] Bleeker | The epistemics of encryption[END_REF]

a principle of permutation B i [α]ϕ → [α]B i ϕ is valid.
Semantically, let B i (w) be the set of worlds agent i considers possible at the actual world w, and suppose R α (w) is the set of possible outcomes of action α when executed in w. Then in the above approaches the worlds i considers possible after α are obtained by applying α to each of the possible worlds in B i (w) ("mentally executing α"), and then restricting the resulting set (B i • R α )(w) in some way.

In the case of actions that do not modify the world but only the agents' belief states ('purely doxastic actions'), the permutation principle roughly speaking expresses that the possible worlds making up i's belief state after α, are a subset of i's possible worlds before α. In terms of AGM belief change operations, every such α is thus an expansion.

Obviously, the permutation principle must be weakened to

¬B i [α]⊥ ∧ B i [α]ϕ → [α]B i ϕ Semantically, we should have that (R α • B i )(w) ⊆ (B i • R α )(w) only when (B i • R α )(w) = ∅.
Sure, this prevents i from running into inconsistencies in our scenario. But this does tells us nothing about i's belief state after inform(j, i, p). In this paper we introduce a permutation principle for this case.

Action laws and enabling actions

If i believes α is inexecutable, this means that there is no v ∈ B i (w) where the preconditions of α's executability laws hold. By an executability law we mean a principle associated to α that has the form A → α . Let pre(α) be the disjunction of all executability preconditions of α. Then all the executability laws for α can be represented by the single pre(α) → α . Now when i learns that α has nevertheless occurred, he cannot just mentally execute α to form his new belief state: he should first change his beliefs about pre(α). We suppose that i does this by adjusting each of his possible worlds so as to make pre(α) true, by mentally executing a particular action whose effect is pre(α).

Formally, we associate to every atomic action α an action enable α , and we say that enable α makes α executable.

We postulate that enable α can always be executed. Hence

pre(enable α ) =
Semantically, for every action there must exist at least one possible world where it is executable. This forces us to exclude actions that are never executable, such as the action inform(j, i, ⊥) of unsincerely asserting an inconsistency. 1As usually done in AI, we suppose that to every action α there is associated a set of effect laws describing the change α brings about. Such laws take the form A → [α]C. For inform-actions we have e.g. the effect law [inform(j, i, ϕ)]B {i,j} B j ϕ expressing that after j informs i that ϕ it is a common belief of i and j that j believes ϕ. We also suppose ¬pre(α) → [α]⊥ to be among α's effect laws.

What are the effect laws for enable α ? As enable α makes α executable, we must have:

[enable α ]pre(α)

We have supposed that sincerity is the only precondition of inform-actions. Therefore pre(inform(j, i, ϕ)) = B j ϕ. Hence we have the effect law [enable inform (j,i,ϕ) ]B j ϕ

Constructing the new belief state

Semantically, if w is the actual world then the worlds i considers possible after a surprising occurrence of α can now be obtained from the worlds in B i (w) by first applying enable α and then applying α.

That is, if (R α • B i )(w) = ∅ then we have (R α • B i )(w) ⊆ (B i • R enable α • R α )(w)
Axiomatically, the permutation principle is weakened to

B i [α]⊥ ∧ B i [enable α ][α]ϕ → [α]B i ϕ
Just as it has been done for the strong permutation principle in the AI field of reasoning about actions [START_REF] Scherl | The frame problem and knowledge producing actions[END_REF][START_REF] Shapiro | Specifying communicative multi-agent systems[END_REF][START_REF] Herzig | A logic for planning under partial observability[END_REF], such a principle can then be combined with existing solutions to the frame problem. (For our permutation principles this has been done in [START_REF] Herzig | Sensing and revision in a modal logic of belief and action[END_REF].) The result is a constructive characterization of the agents' belief states after an action.

Using our permutation principles and the sincerity axiom we can prove with our action laws for inform-and enable-actions that [inform(j, i, ϕ)]B {i,j} B j ϕ is valid.

Discussion and conclusion

Successor state axioms

Inspired by work of Moore [START_REF] Moore | A formal theory of knowledge and action[END_REF], Scherl and Levesque have introduced a situation calculus based framework for resaoning about action and knowledge [START_REF] Scherl | The frame problem and knowledge producing actions[END_REF]. They adopt the strong permutation axiom, and even strengthen it to an equivalence that they call a successor state axiom. They show how such a principle allows for regression, which is a powerful reasoning technique in the situation calculus. In [START_REF] Shapiro | Iterated belief change in the situation calculus[END_REF], Shapiro et col. adapt Scherl and Levesque's successor state axiom to belief, integrating a revision-like operation that is based on plausibility orderings. They define B i ϕ as truth of ϕ in the most plausible among the possible worlds. If a doxastic action eliminates these most plausible possible worlds, then previously less plausible worlds become the most plausible ones. The plausibility ordering should be kept fixed.

While being intuitively appealing, their solution has several drawbacks. (1) As the authors note, it is restricted to deterministic actions. ( 2) "The specification of [the plausibility ordering] over the initial situation is the responsibility of the axiomatizer of the domain." [START_REF] Shapiro | Iterated belief change in the situation calculus[END_REF] This is particularly demanding because (3) in order to guarantee that after α the set of possible worlds is nonempty, the authors require the set of possible worlds to contain enough worlds initially, restricting thus the agent's 'doxastic freedom'. [START_REF] Bleeker | The epistemics of encryption[END_REF] The approach is unsatisfactory when applied to communication. Consider the following example: agent k is competent at p, and j is not. Agent i is completely ignorant initially, and all possible worlds are equally plausible for i. Then (under adequate hypotheses of cooperation) we can expect that when j asserts p, then i adopts p, i.e. [inform(j, i, p)]B i p. Moreover, as all worlds were equally plausible, p holds in every world possible for i. Therefore when subsequently k asserts ¬p, i will unavoidably move to an empty set of possible worlds.

Is this revision?

In our scenario, what i must do is to revise his beliefs about the executability of inform(j, i, p). The normative framework for belief revision being the AGM theory [1], Segerberg's DDL [START_REF] Segerberg | Belief revision from the point of view of doxastic logic[END_REF][START_REF] Segerberg | Two traditions in the logic of belief: bringing them together[END_REF] integrates AGM revision into a doxastic logic. Linder et al. have a similar framework [START_REF] Bernd Van Linder | Actions that make you change your mind[END_REF].

In the case of an expansion action, the agent's belief state is completely determined in DDL (by an axiom that correspond to the above permutation axiom). In the case of revision actions, apart from the AGM persistence postulate there is no principle relating an agent's belief state before and after α. In a sense, we have integrated into DDL a particular operation of revision, which establishes such a link.

Which of the AGM postulates do we satisfy? With a similar encoding as that of Shapiro et col. it can be shown that we satisfy the basic postulates (K*1) -(K*4), and (K*6). (The names of the postulates are as in [START_REF] Shapiro | Iterated belief change in the situation calculus[END_REF]). If we define update actions as in [START_REF] Shapiro | Iterated belief change in the situation calculus[END_REF] we satisfy the update postulates (K 1), (K 2), (K 4), and (K 5) just as there.

This is also the reason why we consider atomic actions only. Indeed, let toss be the action of tossing a coin, let its executability precondition be that one must hold a coin (pre(toss) = holdCoin ), and let one of its effects be that one no longer holds the coin (¬holdCoin ). Then it is impossible to execute toss twice. This means that toss; toss cannot be enabled, illustrating that enable cannot be applied to complex actions.