{'fraction_non_alphanumeric': 0.03167206381802385, 'fraction_numerical': 0.008253692441355343, 'mean_word_length': 4.643637174058391, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Towards Quantitative Magnetisation Mapping

The starting point in quantitative susceptibility mapping (QSM) is a theoretical model that is used to map susceptibility distributions from magnetic field measurements. It requires regularisation techniques to avoid artefacts in the resulting image. The underlying problem is that the model was developed by starting with the so-called Lorentz sphere on a microscopic scale. After averaging over a macroscopic sample, it is assumed that the magnetic flux density vanishes in the center of the sample. For the macroscopic problem of a homogeneous sphere in a uniform field, we show that at the surface the normal component of the flux density is not continuous, which contradicts Maxwell's macroscopic theory. In this paper, we propose a model consistent with macroscopic magnetic field theory, in which we image magnetisation rather than susceptibility. This model is well-posed. Some simple but representative numerical examples show that it allows for high-resolution images.


Introduction
In Quantitative Susceptibility Mapping or QSM, the objective is to image the susceptibility profile within a part of the human body [2], [4], [13], [16]. Its development goes back more than 20 years ago [15] and QSM finds many applications in neuroradiology and neuroimaging, such as MR venography, oxygen saturation imaging, traumatic brain injury imaging, multiple sclerosis, and brain tumor imaging [11].
As is well known, the magnetic susceptibility is a constitutive parameter that relates the magnetisation M to the magnetic field H (to avoid confusion, we refer to H as the magnetic field and to B as the magnetic flux density). The spatial domain data model that is used in QSM to retrieve this constitutive parameter is obtained by starting from the classic field expression of a magnetic dipole [6] and taking the microscopic Lorentz correction into account [8]. The model can also be formulated in the spatial spectral domain (see [20], for example). In either case, the QSM problem of retrieving the susceptibility from magnetic field data is ill-posed and requires regularisation [9]. A mathematical analysis of the QSM model, the effects of particular regularisation strategies, and a discussion on the causes of streaking artefacts in reconstructed susceptibility maps are given in [1].
In the present paper, we focus on the fundamental QSM model and show that within this model the boundary condition for the magnetic flux density is violated. Specifically, we show that the normal component of the magnetic flux density actually jumps across a source-free interface, where the magnetisation exhibits a jump. Clearly, this is not in accordance with Maxwell's field theory.
Furthermore, in contrast to QSM, we propose to image the magnetisation of tissue rather than its susceptibility. We call this Quantitative Magnetisation Mapping (QMM) and we consider this mapping problem in the spatial domain. Similar to QSM, we assume that only scalar magnetic field data in the direction of the background field (usually the z-direction) is available. We show that in the QMM model the boundary conditions for the magnetic field H and magnetic flux density B are not violated. Moreover, we show that the QMM data operator is actually a shifted version of the QSM data operator and that the ill-posedness of the QSM model does not carry over to the QMM model.
Having a tissue magnetisation map available from QMM, the corresponding magnetic field can be determined. In case a model for the susceptibility function is available from microscopic or quantum mechanical considerations, the susceptibility function can be determined, because the magnetisation and the magnetic field are known. The QMM model does not depend on the particular microscopic or quantum mechanical susceptibility model that it used, since it reconstructs the macroscopic magnetisation from which the corresponding magnetic field can be computed. The susceptibility model is only required to relate the magnetisation to the magnetic field. Figure 1 gives a schematic overview of QSM and QMM, illustrating the difference between the two approaches. This paper is organised as follows. In section 2, we briefly consider the basic QSM model in the spatial domain and show a number of reconstruction results based on simulated and analytical data. These results verify that the QSM operator is indeed ill-posed and regularisation is required to obtain a unique and stable approximate solution. Furthermore, we show that within the QSM model, the boundary condition for the normal component of the magnetic flux density is not satisfied. Subsequently, in section 3, we present our QMM approach and demonstrate that within the QMM model the boundary conditions for the magnetic field H and magnetic flux density B are satisfied at an interface between two media with a different magnetisation. In section 4, numerical experiments demonstrating the performance of QMM are presented as well. In addition, we discuss how a susceptibility function of a particular material or tissue type can be reconstructed in QMM. Finally, in section 5, we compare the QSM and QMM operators and show how they are connected. The conclusions can be found in section 6.

Quantitative susceptibility mapping without regularisation
For each location with position vector r = (x, y, z), the mathematical problem in QSM is to solve the susceptibility χ(r) from the integral equation [20] ∆B Here, D ∈ R 3 is the domain where χ(r) = 0. Further, the integral − is a Cauchy principal value integral and ∆B is measured magnetic flux density data due to the presence of a nonzero susceptibility χ. Finally, d is the so-called dipole kernel, which is given by for r = r . If we substitute the right-hand side of (2) into (1) and interchange the order of integration and differentiation, we obtain Including the integration point r = r then leads to the integral equation where we have used that the principal value of the integral is equal to − 1 3 χ(r). Note that in the spectral domain with spectral vector k = (k x , k y , k z ), this equation becomes algebraic, cf. [20], where ∆B and χ are the spectral domain counterparts of ∆B and χ, which is clearly not invertible for all k ∈ R 3 satisfying k 2 z = 1 3 k · k. In the following section, our goal is to check the inversion of (4) in the spatial domain. Specifically, we first solve a forward problem by computing data ∆B from QSM equation (4) for a known susceptibility χ. Subsequently, we take this data as input and try to retrieve the susceptibility again using the same QSM equation. Although we then commit an "inverse crime" (see p. 154 of Colton and Kress [3]), this does allow us to check how well the inverted susceptibility reproduces the input of the forward problem. Finally, we also consider the case where we have analytic data for a homogeneous sphere and use this data to reconstruct the susceptibility of the sphere again using the QSM equation (4).

Iterative inversion based on circular convolutions
Equation (4) is a convolution type integral equation for the unknown magnetisation. Zwamborn and Van den Berg [21] have shown that such an integral equation can efficiently be solved using a weak form formulation together with the so-called iterative CGFFT method. After a spatial discretisation procedure, the discrete spatial convolution is computed as a circular convolution on an extended grid with zero padding. Although sampling in each Cartesian direction is then doubled, the discrete convolution is computed exactly. Here we follow a similar approach, but instead of using the conjugate gradient method, we use BiCGSTAB [19] as an iterative solver, since it has been shown that this method typically exhibits fast convergence for well-conditioned convolution type operators [18]. We realise, of course, that the QSM system is severely ill-conditioned, but for a fair comparison of QSM and QMM, we still use the BiCGSTAB method in QSM. In the discretized formulation, we use a regular grid of 128 × 128 × 128 nonoverlapping subdomains (voxels) with a sampling width of 0.002 m. On each subdomain we assume that the susceptibility χ is constant, and we replace the singular Green function 1/4π|r−r | by its spherical mean (weak form) [18]. In addition, we keep the ∂ 2 /∂z 2 operator in front of the integral over D and replace it by a finite-difference operator. The circular convolution is computed using 256 × 256 × 256 FFT samples and a fixed number of 200 iterations is used in the BiCGSTAB method, which for the current QSM problem is typically sufficient to obtain accurate field or susceptibility approximations provided the BiCGSTAB-method converges, of course.

Some simple test cases
First, we consider a discretised homogeneous sphere with radius a = 10 cm and susceptibility χ = 1, and compute synthetic data using (4). Subsequently, we attempt to retrieve the susceptibility from this data again using (4). In figure 2, we show the exact model images χ (left column), the reconstructed images χ map (middle column), and their difference χ−χ map (right column). In the first row we present the susceptibilities in the cross-sectional plane z = 0. Here, the discrepancies between the exact and reconstructed susceptibility are around 10%. In the second row we present the corresponding values at the cross-section x = 0. These differences are enlarged by the presence of the ∂ 2 /∂z 2 operator in the z-direction. In the caption of the figure, we also present the mean error over the 3D domain D, defined as mean error = χ−χ map D / χ D , where · D is the Euclidean norm on D. We note that increasing the number of BiCGSTAB iterations essentially does not lead to any further improvements in the reconstructions. We can conclude that even in an "inverse crime" situation, application of the standard QSM method in the spatial domain leads to reconstruction artefacts, which are also documented in the QSM literature (see e.g. [1]).
Second, we again consider a homogeneous sphere, but this time a defect is present inside the sphere. The defect is modeled as a smooth ellipsoid, where the susceptibility is given by but the reconstructed defect is clearly visible, although the shape is slightly distorted. As a final test case, we use analytical field data for the homogeneous sphere. In particular, for a sphere with radius a and a susceptibility χ = 1, equation (4) evaluates to (see (16) of [14], and Appendix A) 3z 2 − |r| 2 |r| 5 , for |r| > a and ∆B(r) = 0, for |r| < a. (6) Then, for these analytical data, the iterative inversion of the susceptibility from (4) diverges. The resultant image after 200 iterations is shown in figure 4. This confirms that the inversion in the spatial domain is ill-posed as well. In this regard, we mention that the second equation of (6) shows that for a sphere, any constant susceptibility is in the nullspace of the QSM operator for r ∈ D. Furthermore, Choi et al [1] stated that the solution of the QMS problem is a linear superposition of the actual susceptibility corresponding to the compatible portion of the data plus artefacts generated by the compatibility violators. They concluded that proper exclusion of data compatibility violations may be necessary. An obvious culprit in this case is that ∆B(r) in (6) jumps at the poles of the sphere (|r| = a and z = ±a in (6)), which is contrary to macroscopic Maxwell theory. In absence of magnetic surface sources, the normal component of the magnetic flux density at an interface should be continuous.
In this paper we propose a different inversion strategy. Instead of mapping the susceptibility, we propose to reconstruct the magnetisation within the domain of interest from a knowledge of magnetic flux density. In the next section, we start with classical Maxwell theory and formulate this reconstruction approach.

Quantitative magnetisation mapping
In a source-free smoothly varying medium, a static magnetic field satisfies the homogeneous Maxwell equations with magnetic field H and magnetic flux density B. In case the medium parameters are not smooth and jump across a source-free interface, these equations have to be supplemented by boundary conditions. In particular, the tangential components of H and the normal components of B must continuous upon crossing the interface. The relation between the magnetic flux density and the magnetic field is given by the relation where µ 0 is the permeability of vacuum. The magnetisation is related to the magnetic field, that is, M = M(H) and a constitutive relation describes how exactly M depends on the magnetic field H. This relation obviously depends on the magnetic properties of the medium, but we leave this dependence unspecified for the moment. We start by defining the background field as the field that is present when magnetisation effects are ignored. This field is denoted by {H b (r), B b (r)} and for this field we have the constitutive relation B b = µ 0 H b . Subsequently, we introduce the total field as the field that is present in case magnetisation effects are taken into account. At each location the total field is denoted by {H(r), B(r)}. In case M does not vanish on a bounded domain D, the difference between the total magnetic field and the background field is given by the integral equation, see e.g. [5], which holds for r ∈ R 3 . If the data [H(r) − H b (r)] is known for r ∈ D (that is, data inside the domain D is known), Friedman [5] has proven that in the Hilbert space of vector functions L 2 (D), a unique solution for the magnetisation M exists. We note that at this point the microscopic constitutive model is not important.
To check the continuity of the tangential components of the magnetic field at an interface where the magnetisation jumps, we apply the curl-operator to this integral equation. We observe that ∇ × H = ∇ × H b , since the curl of the gradient term vanishes on either side of the interface. Hence, at an interface with normal vector ν, the tangential components of the total magnetic field ν × H are continuous, since the tangential components of the background magnetic field ν × H b are continuous.
In QSM practice, magnetic flux densities are measured instead of magnetic field values. We therefore rewrite (10) in terms of the magnetic flux density. To this end, we multiply both sides of (10) by µ 0 and use relation (9) to arrive at the integral representation which holds for r ∈ R 3 . For known data [B − B b ] in D, (11) is an integral equation for the unknown magnetisation M in D. Since (9) is a linear relation and (10) has a unique solution, there exist a unique solution of (11).
To check the continuity of the normal components of the magnetic flux density, we use the relations and ∇ 2 (4π|r − r |) −1 = −δ(r − r ) with δ the Dirac function. With the help of these relations, the integral representations for the magnetic flux density can also be written as At an interface (with normal vector ν) where the magnetisation jumps, continuity of the normal component of B can now be verified by applying the div-operator to the representation of (13) for observation points located on either side of the interface. This gives ∇ · B = ∇ · B b , since the divergence of the curl vanishes at both sides of the interface. From this observation, we conclude that at an interface where the magnetisation M jumps, the normal component of the known magnetic flux density ν · B is continuous, since the normal component of the background field ν · B b is continuous. Finally, in QSM only the z-component of ∆B is known through measurements, while polarisation effects are neglected. In our formulation we take this into account as well, and rewrite the magnetisation as M = M i z . Ignoring polarisation effects, the vector integral representations of (10) and (11) then simplify to the scalar equations and Note that the definition of ∆B is slightly different from the one defined in the QSM model. In Appendix A we show that for the special case of a homogeneous sphere, the tangential component of the magnetic field and the normal component of magnetic flux density are indeed continuous across the surface of the sphere when the above data models are used. Our objective is now to reconstruct the magnetisation M from knowledge of ∆B and we refer to this inverse problem as quantitative magnetisation mapping or QMM.

Some simple test cases
In this section, we consider the same test cases as discussed for the QSM model, but instead of using QSM equation (4), we now consider the QMM equation (15). Again, we use the BiCGSTAB method to solve the latter equation.
First, we consider the discretised homogeneous sphere with radius a = 10 cm and magnetisation M = 1, and calculate numerically synthetic data using (15). Subsequently, we use this data and the same model (15) to retrieve the magnetisation. to the finite difference implementation of the ∂ 2 /∂z 2 operator. Increasing the number of BiCGSTAB iterations from 13 to 20, the mean error is reduced further to 0.03%. Second, we consider the homogeneous sphere with the spheroidal defect. The imaging results are presented in figure 6. Again, these results were obtained within 13 iterations and the mean error did not change, which confirms our expectation that the resolution achieved is the same over the whole observation domain.
Third, we use the analytical data for the homogeneous sphere. In contrast to the divergent results of figure 4, for an error criterion of 10 −4 , the BiCGSTAB method converges within 13 iterations. The results are shown in figure 7. Here, errors of half the maximum value of the model magnetisation are visible, but only at the interface of the sphere. The method shows the mismatch between the stair-case approximation and the exact spherical boundary. These differences are magnified by the presence of the ∂ 2 /∂z 2 operator. We further note that increasing the number of iterations reduces the residual error in the equation in question, but the mean error between M and M map does not significantly decrease. We also remark that, in this example, the jump in the magnetisation is large compared to magnetisation changes encountered in practice. For smaller magnetisation variations, the finite difference approximation of the ∂ 2 /∂z 2 operator is typically much better.
Our final example was not considered before and consists of a small sphere and a rectangular shape with smooth edges embedded in a larger sphere as illustrated in figure 8. The distance between the small sphere and the rectangular shape is equal to the sampling width of 2 mm. The reconstructed magnetisation is shown in the middle column of figure 8 and the difference between the exact and reconstructed magnetisation is shown in the third column. We observe that without noise and other modelling errors, QMM is able to provide high-resolution reconstructions. Reducing the sample width results in reconstructions with an even higher resolution, although at the expense of increasing computer memory and computation time.

Constitutive relation
As soon as the magnetisation M = M i z has been determined, the magnetic field H can be calculated from (10) and the constitutive parameter of the material (susceptibility) can be estimated by minimising some norm of the difference between M and H on the domain D. For an isotropic medium, the constitutive relation M = χH, with χ the susceptibility, can only be satisfied when H = H z i z . Fortunately, the magnetisation in QSM practice is much smaller than the background magnetic field, and one can make the approximation H ≈ H b z i z . Assuming that H b z = 0, the susceptibility can then be obtained as In principle, this is the effective susceptibility at the observation point r. When detailed information on a microscopic scale is desired, the constitutive relation of the material or tissue is required. This relation describes the response of the material to the local magnetic field. Strictly speaking, finding such a relation requires a quantum mechanical treatment and after averaging over a suitable domain, the continuum approximation of M(r) at the observation point r is obtained. In other words, a priori knowledge about tissue is required for inversion of a constitutive model. The question that remains is whether a priori knowledge of the microscopic model is required for medical diagnostics, and if so, which model should be used for a certain material or tissue type.

Comparison of the QSM and QMM operator equations
In QSM, the basic imaging equation can be written as where ϕ is the unknown susceptibility χ, f is the data ∆B of QSM, and A is the QSM operator given by (7). Finding ϕ from knowledge of f in D is an ill-posed problem and requires regularisation [7]. Examples of regularisation approaches are Tikhonov regularisation [17], total variation minimisation [12], or regularisation by analytical continuation [10].
In terms of operator A, the QMM equation of interest is given by where this time ϕ is the magnetisation M and f is the data ∆B of QMM. We observe that the QMM operator 2 3 I + A is actually a shifted version of the QSM operator with a shift given by the scaled identity operator 2 3 I. The ill-posed nature of QSM can be traced back to the underlying physical model in which Lorentz correction annihilates this shift and, loosely speaking, regularisation in QSM can be seen as an attempt to restore the shift. Since such an approach may introduce regularisation artefacts and is not necessary in QMM, we propose magnetisation imaging as an alternative to QSM.

Conclusions
It is well known that in QSM the problem of retrieving the susceptibility from knowledge of the (z-component of the) magnetic flux density inside the domain of interest is an ill-posed problem. In this paper we have shown that, in addition to this ill-posedness, the magnetic flux density within the QSM model also does not satisfy the appropriate boundary condition that should hold at an interface between two media with a different magnetisation.
To avoid some of the difficulties associated with QSM, we have proposed an imaging procedure that provides maps of the magnetisation instead of the susceptibility within a region of interest. We call this imaging procedure quantitative magnetisation mapping or QMM and we have shown that within the QMM model, the boundary conditions at the interface between two different magnetic media for the magnetic field H and magnetic flux density B are satisfied. Moreover, a comparison of QMM and QSM shows that the QMM operator is actually a shifted version of the QSM operator. More precisely, the operator of QSM is an integral operator of the first kind, while the QMM operator is an integral operator of the second kind [7]. Numerical experiments for elementary three-dimensional structures have also been presented and show that in QMM the magnetisation can be retrieved in just a few iterations of an iterative solver such as BiCGSTAB. Furthermore, our simple imaging examples exhibit very high resolution over the whole window of observation, which may have significant consequences for medical diagnostics.
With a magnetisation map M at our disposal, the corresponding magnetic field H can be determined from the integral representation that relates the magnetic field to the magnetisation. Subsequently, if a parametric model M (H) is known from either microscopic or quantum mechanical considerations, then the model parameters can be determined by matching the model to the reconstructed magnetisation. Obviously, the mapping procedure in QMM does not depend on a particular material or tissue model, since it directly reconstructs the magnetisation from measured field data.
Finally, we mention that in this paper we have focused on the basic QSM and QMM equations and we used simple 3D structures. We used either simulated or analytic data sets to demonstrate some of the fundamental properties of the QSM and QMM operators. In practice, however, we have to deal with measured data sets containing noise and other (background) field perturbations. The performance of QMM on such data sets will be investigated in future work.

Appendix A. Homogeneous sphere in a uniform field
To check the continuity of the tangential magnetic field and the normal magnetic flux density, we consider the special case of a homogeneous sphere that occupies the spherical domain B with radius a and center at the origin r = 0. We assume that, inside the sphere, the susceptibility χ is constant in the QSM model, and the magnetisation M is constant in the QMM model. For both models, we use the relation For the QSM model, we start with (4). Here, the magnetic flux density is obtained as ∆B = K(r) χ, for |r| > a, and ∆B = 0, for |r| < a. At these pole locations, ∆B is pointing in a direction normal to the sphere, but it is not continuous. The jumps amounts to 2χ/3. This contradicts the macroscopic theory of Maxwell's equations.
For the QMM model, we start with (14) and (15). The magnetic field is obtained as ∆H(r) = K(r) M, for |r| > a, and ∆H(r) = K(r) M, for |r| < a. which shows that the normal component of the magnetic flux density is continuous as well. Obviously, the QMM model is compatible with the macroscopic Maxwell equations.