Model-Based Reconstruction for Magnetic Particle Imaging in 2D and 3D

We contribute to the mathematical modeling and analysis of magnetic particle imaging which is a promising new in-vivo imaging modality. Concerning modeling, we develop a structured decomposition of the imaging process and extract its core part which we reveal to be common to all previous contributions in this context. The central contribution of this paper is the development of reconstruction formulae for MPI in 2D and 3D. Until now, in the multivariate setup, only time consuming measurement approaches are available, whereas reconstruction formulae are only available in 1D. The 2D and the 3D (describing the real world) reconstruction formulae which we derive here are significantly different from the 1D situation -- in particular there is no Dirac property in dimensions greater than one when the particle sizes approach zero. As a further result of our analysis, we conclude that the reconstruction problem in MPI is severely ill-posed. Finally, we obtain a model-based reconstruction algorithm.

1. Introduction. Magnetic particle imaging (MPI) is an emerging imaging modality which was developed by Gleich and Weizenecker in 2005 [6]. Its goal is to determine the spatial distribution of magnetic nanoparticles by measuring the non-linear magnetization response of the particles to an applied magnetic field. For example, in angiographic applications, the particles are injected in humans' systems of blood vessels in order to image this system by determining the concentration of the particles.
The particles are not observed directly, but a certain current induced by them is measured. More precisely, the non-linear response of superparamagnetic nanoparticles to a temporally changing applied magnetic field is exploited: Changing the magnetic field results in a change of magnetization locally near a point, say x; the local change of magnetization induces a current in a set of receive coils where the current reflects the concentration of the particles near x.
MPI is a very promising imaging modality for biomedical diagnostics because of two features: firstly, MPI offers a high dynamic spatial and temporal resolution [17]. Secondly, in contrast to many other tomographic methods such as PET or SPECT (cf. [23,32]), MPI does not employ any ionizing radiation. In particular, patients are not exposed to radioactive substances in the imaging process. Furthermore, only the magnetic particles are imaged which avoids artefacts stemming from nearby tissue as in angiographic CT. Potential applications of MPI are any tracer based diagnostics, such as blood flow imaging and cancer detection [18,Ch. 7] as well as quantitative stem cell imaging [35]. A basic reference on MPI is the just mentioned book [18], as well as the thesis [15]. Concerning early work, we exemplarily refer to [7,34,27]. As further references we would also like to mention [16,19,20,30,11,29,22,5]. These references also include references to related imaging setups and to related questions concerning particle design as well.
In a proper imaging setup, the MPI forward operator is well modeled as a linear operator which has to be inverted (in the sense of inverse problems) to reconstruct the signal from the measured voltage. The present approaches to represent this forward operator can be categorized into measurementbased approaches, e.g., [33,28,24] and model-based approaches, e.g., [21,12]. In the measurement based approach, a basis (or, more general, a dictionary) is considered and, for each member, the system response is measured physically. This is rather expensive; in particular, when reconstructing with higher resolution. This is actually a time limiting factor in practice. Knowing these measurements, the reconstruction consists of the (regularized) solution of the corresponding system of linear equations with the rows of the corresponding matrix being the measured system responses. In order to reduce the time consumption, a model-based approach is desirable.
Presently, there are basically three approaches to derive a model for the MPI forward operator which are intimately related. They have been developed over the past five years; for a detailed overview, we refer to [12]. All these models employ Faraday's law of induction with respect to a volumetric coil; see [18,15] for a derivation. Combined with the Langevin model of the magnetization response a description of the received signal is obtained. The applied magnetic field is in general a combination of a static field and a time-dependent drive field. The latter determines the movement of the so-called field free point that scans the test object.
Let us be more precise. In [27], Rahmer et al. consider a 1D scenario and study the Fourier series of the signal derived from a delta distribution. They conclude that, in the case of ideal particles, the system function can be written in terms of Chebyshev polynomials. In the 2D/3D scenario such a representation is not available [28]. In contrast to [27], Goodwill and Conolly [9,10] and Schomberg [31] follow a geometric approach. In [9,10], the received signal is studied in terms of the trajectory of the field free point whereas the approach of [31] studies the integral of the received signal instead of the signal itself. The first authors observe that the field free point in the spatial domain (or x-space in their terminology) is uniquely determined by the magnetic drive field and vice versa.
Their main result is that the signal generation process can be described by a certain kind of convolution operator based on a matrix-valued point-spread function or kernel (cf. the "Generalized MPI signal equation" in [10]) involving the field-free point and the tangent to its trajectory.
Schomberg [31] obtains a formulation of convolution type in which the kernel describes the magnetization response of the system. In contrast to [9,10], the kernel is vector-valued and a different quantity is described. For the 1D case, it is shown by Grüttner et al. [12] that the x-space formulation is mathematically equivalent to the frequency space formulation of the model in [27].
Albeit some previous approaches allow for arbitrary trajectories, to our knowledge, no previous approach is detached from them. Furthermore, in contrast to 1D, there are no reconstruction formulae and, thus, no algorithmic approaches based on such formulae are available. The previous model-based reconstruction algorithms in 2D and 3D mainly consist of calculating the entries of the corresponding matrices instead of measuring them.
Certainly, reconstruction formulae in 2D and 3D are very desireable. We believe that such formulae have not been derived yet for the following reason: in contrast to 1D, the problem cannot be modeled as a standard time-invariant system which would result in a standard deconvolution problem. Instead, a matrix-valued kernel enters the scene which modifies the problem in such a way, that, in the recent years, no solution was figured out -although desired. As it turns out in this paper, a simple but striking observation allows us to develop an approach which solves the problem. Indeed, applying this striking observation casts the problem into a simpler problem one can deal with easier.
Contributions. In this paper we make three contributions: (i) a unified spatial model in terms of an MPI core operator; (ii) reconstruction formulae for multivariate MPI in 2D and 3D; (iii) reconstruction algorithms for multivariate MPI in 2D and 3D.
Concerning (i), we revisit the modeling of the MPI time signal generation process and decompose it. We identify an MPI core operator (2.15) acting on phase space. The phase space point of view has the advantage of being totally free of the concept of any trajectories. Trajectories only enter the scene as a means of getting data for the operator equation involving the core operator. This point of view is essential for understanding the spatial action of an MPI scanner on a spatial particle distribution (the image of interest) and is key for the derivation of our reconstruction formulae. The MPI core operator is formulated in a common form for all three cases 1D, 2D and 3D.  [17]). The rectangle in the center of the scanner is the field of view (FOV). Two permanent magnets and the outer coil pair aligned with the y-direction generate the selection field H S . Both outer coil pairs are used to generate the drive field H D which moves the field free point (FFP) in the FOV for example on Lissajous trajectories. The inner coil pairs record the signal.
The most important contribution of this paper is the derivation of reconstruction formulae (Theorems 3.2 and 3.5) for multivariate MPI in 2D and 3D (ii). This is a crucial step to bring a modelbased reconstruction approach for dimensions higher than one to the MPI practinoners. Until now, such formulae only existed in 1D (cf. [27,28]), while for 2D and 3D only extremely time consuming measurement approaches have been employed in these practically relevant situations. These formulas are a result of our analysis of the MPI core operator. Our analysis shows that the 1D situation is significantly different from the 2D and 3D situation; we also emphasize that the 2D and the 3D situation themselves are also significantly different. This difference is owed to a close relationship between the MPI core operator and the fundamental solution of Laplace's equation which we highlight in Theorem 3.1. Our analysis also reveals the severe ill-posedness of the MPI reconstruction problem (Theorem 3.8) which explains why regularization is needed.
Taking this into account, we build reconstruction algorithms (Algorithm 1) for multivariate MPI which are based on our reconstruction formulae (iii). In this way, we provide a model-based reconstruction approach for the 2D and 3D case. To our knowledge, this is the first approach based on a reconstruction formula for the multivariate case.
Outline of the paper. In Section 2, we provide a thorough mathematical description of the MPI signal encoding. In particular, we provide a trajectory free description revealing the core of the MPI signal generation process. In Section 3, we analyze this MPI core operator. Here, we consider first a high resolution idealization limit in Section 3.1. Using the derived results, we consider the nonideal case in Section 3.2. In Section 4, we apply the derived reconstruction formulas for the MPI core operator to establish a reconstruction algorithm. In Section 4.1, we start with the discretization; then, we propose a reconstruction algorithm in Section 4.2; finally, we illustrate our algorithm with a numerical example in Section 4.3.

Mathematical
Model of the MPI Signal Encoding. We start by providing the basic ideas and principles of MPI in Section 2.1. We discuss briefly the one dimensional situation of particle detection in Section 2.2 and then turn to the multivariate situation which is the central topic of this paper. In Section 2.3 we derive the model describing the MPI signal generation. Then we reveal the core of the MPI signal generation process in Section 2.4. In Section 2.5 we relate to the existing models.
2.1. Basic Ideas and Principles. The idea of magnetic particle imaging is to leverage the nonlinear response of superparamagnetic nanoparticles to a temporally changing applied magnetic field H(x, t). Figure 2.1 shows a schematic of a 2D MPI scanner. The combination of the permanent magnets and the outer two coil pairs generates the magnetic field H(x, t) while the inner coil pairs record the signal. (A similar set-up for a 3D MPI scanner is possible by using three outer and three inner coil pairs.) The recording coils receive a voltage u(t) which is described by Farraday's law of induction as the temporal change of the magnetic flux Φ(t), i.e.
The flux Φ(t) of (2.1) is caused by the applied field H(x, t) as well as the magnetization response M(x, t). The parameter µ 0 denotes the magnetic permeability and R(x) ∈ R 3×3 the sensitivity pattern of the three recording coil pairs. The second important principle is the Langevin theory of paramagnetism [3,13]: it gives a description of the magnetization response M(x, t) for superparamagnetic nanoparticles in terms of the Langevin function L: Here m is the magnetic moment of a single particle, ρ is the concentration of the particles. Moreover, the magnetization field M is aligned with the applied field H and vanishes if H is zero. The goal of magnetic particle imaging is to reconstruct the spatial particle concentration ρ from the received time signal u(t). Since only the magnetization response depends on ρ the signal u(t) is reduced to the wanted signal s(t), by subtraction of a reference voltage [31] obtained from a scan of an empty field of view (FOV), i.e., with M = 0.

The MPI Principle of Particle Detection in 1D.
The effect of the non-linear response of superparamagnetic nanoparticles is explained best in a simplified one-dimensional scenario, cf. Figure 2.2. Here, we assume a point sample is located at position x, i.e., ρ(y) = δ(x − y). Moreover, we have a single pair of selection and drive coils such that H(x, t) = H(x, t)e 1 and a single pair of recording coils such that R(x) = R(x)e T 1 . With this configuration the signal is where we used the fact that the Langevin function is a sigmoid (yellow curve in 2.2). We assume now that H(x, t) as function of time t oscillates with small amplitude about a reference valueĤ. If the reference valueĤ is about zero we will receive a non-trivial signal because of the non-linearity of the Langevin function near zero. Otherwise ifĤ is large in magnitude we will receive the zero signal since the Langevin function is almost constant away from zero. Thus, control over the applied field allows for the detection of the point sample.

The MPI Signal Generation in Higher Dimensions.
Let us now consider the general multivariate situation. In contrast to the univariate setting where one scans along a straight line, it becomes important for the multivariate setting that the signal is obtained by scanning along a trajectory which covers an area or a volume.
Scanning. By design, the applied field H is the superposition of a static selection field H S and a dynamic drive field H D : (2.5) The static part H S is generated by a direct current in the selection field coils. The dynamic part is generated via an alternating current in the drive field coils and is described by Here I(t) denotes the electric current in the coils and P(x) is the sensitivity profile of the drive field coils according to the Biot-Savart law [3,13]. The locations where H ≈ 0 are important for the detection of particles (cf. Section 2.2) which takes us directly to the definition of the field free point (FFP) denoted by r(t). The FFP is the geometric locus where H vanishes: From the second equality we learn that the current I(t) drives the FFP to scan the object. Because the FFP trajectory r(t) describes the locations where particles can be detected, the hull of r(t) (under the assumption that r(t) fills a portion of space very well) defines the actual field of view. Signal Encoding. Schomberg [31] points out that by solving the right-hand side of (2.7) for I(t), I(t) = −P(r(t)) −1 H S (r(t)), the applied field H can be described in terms of the FFP rather than the current: By design, the sensitivity pattern P is almost homogeneous in the FOV which simplifies the expression to The static field H S is generated by permanent magnets and/or a selection coil pair in Maxwell configuration [3,13]. It is therefore linear in x with a constant gradient G, i.e., (2.10) The factor g denotes the nominal gradient of the static field. Finally, the linearity of H S implies that the applied field is of the form [31,10] H(x, t) = G(x − r(t)), and also that I(t) = −P −1 Gr(t). (2.11) Again by design, the sensitivity pattern R of the recording coils is also almost homogeneous in the FOV. By substituting (2.11) into (2.2) and (2.3), the MPI signal is given by 2.4. Decomposition of the Signal Generation Process and Derivation of MPI Core Operator. Until now, the signal is described in terms of the trajectory of the field free point. Now, we see that actually less information is required; it turns out that only phase space information is needed meaning that the signal is uniquely determined by the location and the tangent to the trajectory. This reveals the core part of the MPI signal generation process, and makes the model independent of the particular trajectory under consideration.
In order to extract the essential part in (2.12), we employ the spatial flux function Φ defined by Using Φ, we can rewrite (2.12) as (2.14) We derive the MPI core operator MPI[ρ] as the part which operates on ρ, i.e., We see that MPI[ρ] is a spatial convolution of ρ with a matrix-valued kernel. The convolution evaluated at the FFP r(t) produces a matrix which is a applied to the velocity vector v =ṙ(t) of the FFP. The central advantage of (2.15) is that it is completely independent of a particular trajectory, it describes MPI in space and gets rid of all decorative elements. Plugging in a trajectory r(t) into (2.15) the signal is obtained by A basic first step in this direction was already done by Goodwill et al. [10] in their x-space formulation where they relate the signal description with spatial convolution.
Scaling. Now, we want to determine a single dimensionless parameter h which characterizes the resolution of the MPI scanning setup. Because the MPI core operator (2.15) depends on many different parameters, namely the nominal gradient g, the size L of the field of view, and a saturation parameter H sat of the particles (which in turn is a function of the temperature T , the particle diameter d, and the saturation magnetization M sat ), we combine them in a single parameter h. The smaller the parameter h the better is the resolution. We will see here that in a typical setup h is about h ≈ 10 −2 . In Section 3 which is concerned with the analysis of the MPI core operator we will study the idealization limit h → 0 and its relation to the case of h > 0. The support of the particle concentration function ρ is usually contained in the field of view Ω and of the same size, which determines the length scale. As the field of view is a cube, we write it as a linear deformation of the standard cube [−1, 1] 3 withĜ as of (2.10) and the length scale L of the physical field of view. We use this now to transform variables in (2.12) by That isx andr are elements of a non-dimensional field of viewΩ = [−1, 1] 3 . Moreover, we write the receive coil sensitivity as R = R 2R with a dimensionless matrixR which gives only the receive pattern. Equation (2.12) transforms then to where c is a dimensional constant, while h is now a dimensionless shape parameter. By refering to the transformed spatial flux function we have The spatial resolution of an MPI scanner depends on how well the sigmoidal Langevin-function approximates the sign-function on the relevant domain. By the non-dimensionalization above we obtain, with s := |r −x|, s ∈ [0, 2 √ 3], the diameter ofΩ as the relevant domain of L(s/h). The size of h determines the shape of L(s/h) on that domain. With smaller values for h the Langevinfunction gets closer to the sign-function. We have the following relationship between h and the involved physical parameters:  giving a quite good approximation of the sign-function, see Figure 2.3.
Influence of the Electrical Current on the Scan Trajectory. The scan trajectory is controlled via the electrical current I(t). In dimensional form we have (2.25) By decomposing the current as well as the drive coil sensitivity into magnitude and direction, i.e. I(t) = IÎ(t) and P = P 2P we obtain a dimensionless form of (2.25): That means in order to scan a FOV of length scale L we would need a magnitude of for the electrical current (assuming thatP does not shear too much). Given that, then, up to a linear transformation, the dimensionless scan trajectorŷ inherits the signal shape ofÎ(t).
Lower Dimensional Particle Distributions. So far we have considered the signal generation and the MPI core operator in 3D. But for lower dimensional particle distributions we can employ lower dimensional models, i.e. 1D or 2D. Here, we provide a common canonical form valid for all three cases.
We note that the physics of MPI is intrinsically 3D. Hence, even if the particle distribution is lower dimensional, we initially have to start from a 3D setup. Let us first assume that we are given a one-dimensional particle distributionρ(x) = ρ(x 1 ) δ(x 2 ) δ(x 3 ). Here, one would use a single pair of recording coils, i.e.R = e T 1 , and a single pair of drive coils, i.e.r =r 1 e 1 . In that case we have that s 1 (t) = d dr 1Φ 1 (r 1 )ṙ 1 (t), wherē Accordingly we proceed in the case of two-dimensional particle distributionsρ(x) = ρ(x 1 ,x 2 ) δ(x 3 ). In summary, we get the canonical form with a full rank matrixR ∈ R n×n andr ∈ [−1, 1] n , which is valid for any dimension n ∈ {1, 2, 3}.
2.5. Relation to the Previously Proposed Models. Rahmer et al. [27] consider a 1D scenario and study the Fourier series of the signal derived from a delta distribution. In the idealized 1D case, they derive a representation of the convolution kernel in terms of Chebyshev polynomials. In the 2D and 3D scenario they study the Fourier transform of the signal obtained from inserting the Lissajou trajectory into the MPI core operator. But a similar representation of the convolution kernel as in 1D is not available, see [28].
In contrast to [27], Goodwill and Conolly [9,10] study the received signal in terms of the trajectory of the field free point. In [9,10] no Fourier transforms are considered and a general trajectory is inserted into the MPI core operator. For the 1D case, it was shown in [12] that the x-space formulation is mathematically equivalent to the frequency space formulation of [27].
The approach of Schomberg [31] studies the integral of the received signal instead of the signal itself. The result is a formulation of convolution type involving the spatial flux Φ of the MPI core operator.
Summing up, all previous approaches are based on studying the concatenation of the MPI core operator (2.15) or the related flux (2.13) with some individually different operators.
3. Analysis of the MPI Core Operator. As seen in the previous section, the central task in analyzing MPI is to study the MPI core operator A h [ρ] which operates on the particle density ρ . The topic of this section is to start this investigation. The MPI core operator is the derivative of the spatial flux function Φ. The starting point for our analysis is (2.30), where we from now on omit hats and bars in the notation. Moreover, to have a handy notation, we consider the signal multiplied byR −1 and call it s again.
Precisely, in this section, we consider Then the mathematical form of the MPI core operator is The operator A h acts on the density function ρ and produces the matrix-valued function A h [ρ] which is a function of space r.
The time-dependent signal s(t) is obtained by inserting the scan trajectory in the spatial argument r = r(t) and multiplying by the velocity vector v =ṙ(t), i.e., We see from (3.2) that the information on ρ is encoded in the matrix-valued function A h [ρ](r), whereas the vectors v in turn depend only on the employed scan trajectories which are independent of ρ. This observation motivates the focus on the core operator A h and the relationship between A h [ρ] and ρ which we investigate here. By scanning the object with different trajectories we can in principle gain full knowledge about A h [ρ](r); practically, denser trajectories which visit each point several times allow for better approximations.
We start out in Section 3.1 to study the MPI core operator A h in the limiting case h → 0. This reveals the central part of the core operator by freeing it from its blurring part; in particular, we observe dependence on the dimension and see that the previously conjectured Dirac property only holds in 1D but is no longer valid in a multivariate setting. Then in Section 3.2, we consider the non-idealized situation h > 0 taking the blur into account.
3.1. The Idealized MPI Core Operator -the High Resolution Limit h → 0 . We learned in Section 2 that the resolution parameter h is of size 10 −2 . Hence, the question of what happens in the ideal high resolution limit as h → 0 is also interesting and natural from an applied point of view. From a theoretical point of view, it turns out that there are several relations between the idealized and the non-idealized MPI operator: e.g., in the physically most important 3D case, the trace of A h for a particular h is just the concatenation of the trace in the idealization limit with a blurring operator. We will see this in Theorem 3.7 later on.
Equation (3.2) tells us that A h is a convolution operator which operates linearly on ρ. In 1D, the kernel is scalar yielding a standard convolution operator and a standard deconvolution problem. In contrast the multivariate setting comes with a matrix-valued kernel giving a non-standard convolution operator. But, as can be seen in (3.2), the components of the matrix-valued kernel act seperately on the scalar-valued ρ. The striking point here is that the trace of A h [ρ] -as we will see -contains already all the information about ρ. This has not been realized before but is the crucial step because application of the matrix-trace to A h [ρ] gives a scalar deconvolution problem. These facts allow for the derivation of reconstruction formulae. Moreover, we notice an interesting connnection to the fundamental solution of the Laplace equation. In the multivariate setting, n > 1, the limiting convolution kernel is

5)
and (3.4) holds, for any ρ in the space BV 0 (Ω) of functions of bounded variation with zero boundary, where Ω is the field of view.
In the univariate setting, n = 1, we have a Dirac-kernel

6)
and (3.4) holds for continuous functions ρ with zero boundary. Depending on the spatial dimension n ∈ {1, 2, 3}, we have three different relations of the convolution kernel to the fundamental solution Φ(r, x) of the Laplace equation, −∆ r Φ = δ(r − x). In 3D, we have In 2D, we have In 1D, we have Before proving this statement, we record some immediate consequences. The first consequence is the following reconstruction formula. (3.10) That means, we take the pointwise trace and then deconvolve w.r.t. the dimension-dependent κ given by (3.7),(3.8), or (3.9). Proof. Given the matrix-valued data F = A 0 [ρ], we apply the matrix-trace. As a consequence of Theorem 3.1 we obtain the deconvolution problem trace F = κ * ρ. (3.11)

By inversion we obtain
A next consequence is the following statement on ill-posedness in the sense of inverse problems [4,14]. Proof. The degree of ill-posedness of the fundamental solution of the Laplace operator is two by the properties of the Laplace operator [14]. Then, the statement is a consequence of (3.7) and (3.8) for the respective dimensions.
Having recorded these consequences, we turn to the proof of Theorem 3.1. Proof. We consider the MPI core operator A h [ρ] and notice that we can swap the gradient w.r.t r for a gradient w.r.t. x (giving a minus sign) to obtain Then, the matrix-trace applied to A h [ρ] turns the gradient into the divergence operator For the next steps we assume ρ ∈ C 1 c (Ω), which we relax later. By assumption, ρ has compact support in the FOV, and we apply integration by parts to get As the integrand is dominated by |∇ x ρ(x)|, we obtain by Lebesgue's theorem on dominated convergence that Finally, we apply integration by parts again and undo the swap of derivatives which results in We now distinguish between the multivariate and the univariate case. The multivariate case, n > 1. The limiting convolution kernel can be expanded to This may be seen by the following two steps: the Jacobian of the function r/|r| is where I is the identity matrix. The divergence div r (r/|r|) = (n − 1)/|r| is then given by the trace of the Jacobian. The dependence on the spatial dimension n is because of trace I = n. For n > 1 the result (3.16) holds also for functions ρ ∈ BV 0 (Ω) because these can be approximated by functions from C ∞ c (Ω) in the strict topology on BV 0 (Ω); cf.
[1]. The univariate case, n = 1. Here we have This is because for n = 1 the divergence expression comes down to The last equality is to be understood in the sense of distributions.) For n = 1 the result (3.16) holds also for continuous functions ρ ∈ C 0 (Ω) because these can be approximated by functions from C ∞ c (Ω) in the norm topology on C 0 (Ω). Finally, we relate the convolution kernels to the fundamental solution of the Laplace equation Depending on the spatial dimension n ∈ {1, 2, 3} we have three different scenarios: a) In 3D, n = 3, the fundamental solution is given by Φ(r, x) = 1 4π |r−x| , while the convolution kernel is In 2D, n = 2, the fundamental solution is given by Φ(r, x) = − 1 2π log(|r − x|) and its gradient w.r.t r is The convolution kernel is c) In 1D, n = 1, the fundamental solution is given by Φ(r, x) = −|r − x|/2. The convolution kernel is a Dirac and thus related to the second derivative of the fundamental solution: This was the last assertion of the theorem and completes the proof. The Dirac property in 1D has already been reported by Gleich and Weizenecker in [6], and was until now assumed to hold in higher dimensions, too. The news is that the convolution kernels with parameter h in higher dimensions are not Dirac sequences. In contrast, we end up with non-trivial convolution kernels which are lower-order derivatives of the fundamental solution of the Laplacian. This makes sophisticated deconvolution techniques necessary even in the ideal high resolution limit, as we have seen in Theorem 3.3.

The MPI Core
Operator -the Non-Ideal Case h > 0. In concrete applications, it is, for sure, important to reconstruct in the non-idealized case h > 0. Loosely speaking, it will turn out that a non-zero resolution parameter h > 0 results in blured data.
We start out by deriving a statement analogous to the first part of Theorem 3.1 which shows, that, also in the non-idealized case, the trace α h [ρ] of the MPI core operator can be written as a convolution operator with a scalar-valued convolution kernel κ h . In particular, we see that this convolution kernel is an analytic radially symmetric function. The analyticity makes the reconstruction problem even severely ill-posed from an inverse problem point of view.
Theorem 3.4. The trace α h [ρ] of the MPI core operator A h [ρ] is given by

26)
where the convolution kernel κ h takes the form where L is the Langevin function. In particular, f is analytic with singularities only on the imaginary axis. Near zero, f has the power series expansion with a convergence radius of π. (B l denotes the l-th Bernoulli number.) Because of this the kernel κ h is itself analytic and (3.26) holds for any ρ ∈ BV 0 (Ω) independently of the spatial dimenion n ∈ {1, 2, 3}. We postpone the proof of this statement to the end of this section and record some immediate consequences. The first consequence is the following reconstruction formula.
Theorem 3.5. (Reconstruction Formula for the Non-Idealized Case.) Suppose that noise-free data F(r) = A h [ρ](r) at each point r is given. Then, That means, we take the pointwise trace and then deconvolve w.r.t. κ h given by (3.27). The inverse of the convolution operator κ h is denoted by κ −1 h . Proof. Given the matrix-valued data F = A h [ρ], we apply the matrix-trace. As a consequence of Theorem 3.4 we obtain the deconvolution problem trace F = κ h * ρ. (3.30) . A next consequence is the following statement giving a first relationship between the convolution kernel κ h and the idealization limit κ. It reveals again the dimension dependence from a different viewpoint.
Corollary 3.6. (Relation to the Idealization Limit I.) The kernels κ, κ h representing the traces of the idealized and the non-idealized MPI operator are related by where L is the Langevin function.
In the univariate case, n = 1, the second summand of (3.31) vanishes, while the first summand is a Dirac-sequence: In the multivariate case, n > 1, the factor (n − 1)/|y| of the second summand is already the idealized convolution kernel κ(h). The first summand is a zero-sequence, while the second summand tends to κ(y) = (n − 1)/|y|.
The limits here are understood in the sense of distributions.
Proof. This is an immediate consequence of Theorems 3.4 and 3.1: Theorem 3.4 yields the trace of the MPI core operator in terms of a convolution α h [ρ] = κ h * ρ with Theorem 3.1 gives the limit In the univariate case, n = 1, the second summand of (3.34) vanishes. The limiting kernel is κ(r − x) = 2 δ(r − x), thus the first summand is a Dirac-sequence.
In the multivariate case the factor (n − 1)/|r − x| is the limiting convolution kernel while the factor L (|r − x|/h) tends to one almost everywhere, thus the first summand is a zero-sequence.
A particularly nice relation between the traces α[ρ] and α h [ρ] of the idealized and the nonidealized MPI operators is the already announced connection in terms of convolution. It tells us that in 3D the non-idealized α h [ρ] is a smoothed version of the idealized α[ρ].
Theorem 3.7. (Relation to the Idealization Limit II.) In the three-dimensional case, n = 3, we have the following relationship between the non-idealized α h [ρ] and the idealized α[ρ]: In 3D, we can replace the Dirac-δ with the negative Laplacian of the fundamental solution. In addition we move the Laplace operator over to κ h to obtain By Theorem 3.1, we have κ = 8π Φ. Now we express Φ in terms of the idealized kernel κ which results in Finally, we use the result of Theorem 3.1 stating that κ * ρ = α[ρ], which yields equation (3.36). Furthermore we obtain the following statement on severe ill-posedness in the sense of inverse problems [4,14].
Theorem 3.8. (Severe Ill-Posedness.) The non-idealized MPI problem of Theorem 3.4 is severely ill-posed in the following sense: there are no two spaces H s , H t in the Sobolev scale, such that the trace α h of the MPI operator induces an isomorphism α h : H s → H t between these two spaces.
Proof. This is an immediate consequence of the analyticity of the kernel κ h and the fact that an isomorphism has a closed image. From Theorem 3.7 we learn that α h [ρ] is analytic for any ρ; hence the image of α h is necessarily not closed in H t .
Finally, we supply the proof of Theorem 3.4.

Proof. By applying the matrix-trace to
The convolution kernel in the non-idealized case is thus given by By employing the product rule for the divergence, the kernel can be expanded to We express the kernel via the scalar function f as From the Laurent series of the coth we learn that the Langevin function L has the power series where B 2k are the Bernoulli numbers. It is analytic at zero with radius of convergence π. Consequently, the function f is also analytic near zero; it can be written as All further singularities of coth(z) (besides the one at 0 which was involved above) are nonzero and lie on the imaginary axis. This implies that the singularities of f are nonzero and lie on the imaginary axis as well. Using y = r − x as short hand we see that the convolution kernel κ h (y) is analytic near y = 0 by And away from y = 0 it is concatenation of analytic functions.

A Model-Based Reconstruction Algorithm for Multivariate MPI.
In this section, we apply the derived reconstruction formulae for the MPI core operator to establish a reconstruction algorithm for MPI in the multivariate (2D and 3D) cases. We start with the discretization in Section 4.1. In particular, we recast the discrete time data into spatial data and obtain a reconstruction problem for the discretized MPI core operator. Then we present our reconstruction algorithm in Section 4.2. Finally, we illustrate our algorithm with a numerical example in Section 4.3.

Discretization and Reduction to a Reconstruction
Problem for the Discretized MPI Core Operator. As described in Section 2, the measured data is the time-dependent data s(t) which is associated with a scan trajectory r(t). The scan trajectory inherits the shape of the electrical current I(t) given in (2.28) by r(t) = P I(t).
(4.1) (Here we refer to dimensionless quantities, but omit hats.) One advantage of the present approach is that it is independent of the particular trajectory type employed. Furthermore, in our approach, even different trajectories might be combined. The important point is that the data is given as living on a discrete sampling of the trajectory r(t), i.e., we observe values s(t k ) at location r(t k ), with the trajectory having tangent v(t k ) at time t k , for finitely many measurements indexed by k. These data are a discrete sampling of the MPI core operator MPI(r k , v k ) = A h [ρ](r k )v k . Our first task is to employ these data to reconstruct a discrete version of the trace α h [ρ]. To reconstruct these trace data in a discrete setup, we consider discrete functions defined on a grid of N 1 × N 2 cells. (We explain the 2D setting in the following; the 3D setting is obtained by straight forward modifications.) On each cell a function is constant and each cell is represented by its cell center point. We assume that the cells (corresponding to the chosen spatial resolution) are chosen in a way such that each cell of the grid is being passed at least twice from different directions with the corresonding trajectory points r(t k ) hitting the cell. In the following, we say that a time sample t k belongs to cell i, if r(t k ) is in this cell. For each cell i we collect the signal data s(t k i ) from those times samples t k i , such that t k i belongs to the cell I (i.e., r(t k i ) hits the cell i), and gather them in a matrix S i . Accordingly, we collect all the velocity vectorsṙ(t k i ) = v(t k i ) and gather them in a matrix V i . Associating all points that fall in the same cell with the cell center point x i , and then using the signal encoding model, we have for each cell i, the following matrix fitting problem w.r.t. A i , (4.2) By our assumptions above it is ensured that V i has full rank. (For sure, it is possible to relax the above assumption considerably by involving spatial regularity assumptions; see our discussion on future work.) We solve the system (4.2) by least squares fitting : we let V T i = Q i R i be the reduced QR-decomposition [8] with a full rank 2 × 2 matrix R i . We obtain A i and the corresponding trace data u i on cell i (corresponding to α h [ρ](x i )) by

Reconstruction Algorithm.
Our reconstruction scheme may be subdivided into two major steps: in the first step we produce spatial data from the given time data on the trajectory by (4.3). This is an intermediate step which yields the trace data and represents the particle distribution ρ convolved with the kernel κ h according to (3.27). In the second step we perform deconvolution, i.e. inversion of κ h , which we do in a regularized way since we learned in Theorem 3.8 that this problem is ill-posed. In this second step, we employ the analytic results of Section 3.
Step 1: Deriving Trace Data on a Spacial Grid from the Raw Input. We obtain a grid function representing trace data u i in each pixel (grid cell) as explained in Section 4.1.
Step 2: Reconstruction of the Signal from the Derived Trace Data by Deconvolution. Having obtained data u i related to the trace α h [ρ](x i ) from Step 1, we now take the convolution description of α h [ρ](r) derived in Section 3 into account: if the data were noiseless, we would have to solve the problem of finding ρ in κ h * ρ = u given u. However, normally the data s, and therefore also u, are noise-contaminated, and we have learned in Section 3 that this reconstruction problem is severely ill-posed. Therefore, regularization is needed; cf. [4,26,2,14]. As regularization technique, we here input : Time dependent samples s k = s(t k ) along trajectory r k = r(t k ) with tangent v k =ṙ(t k ) at times t k ; regularization parameter µ. output: Reconstructed particle density ρ (solution of (4.5) with preprocessing given by (4.3),(4.2)).
if r k in cell i then // Regularized deconvolution of the trace data using (4.5) by // solving (4.6) with conjugate gradients (CG). ρ = CG(−µ L + K 2 h , K h u); use a variational approach based on classical Tychonov regularization. More precisely, we consider the following (spatially continuous) problem: The convolution with κ h is discretized on the considered N 1 × N 2 grid using the summed-up midpoint rule. We denote the resulting matrix by K h . Moreover, we discretize the gradient with forward finite differences to obtain a corresponding matrix D. The discrete Tychonov-regularized problem is now of the form with vectorized ρ, u ∈ R N 1 ·N 2 . We solve the corresponding discrete Euler-Lagrange equation Here, −L = D T D is the five point stencil discretization of the Laplacian with zero Dirichlet boundary and K h = K T h is symmetric. The whole reconstruction algorithm is subsumed as Algorithm 1. 4.3. Numerical Results. We demonstrate the potential of our method by applying the developed method to synthetic data. Detailed numerical investigations are the topic of a forthcoming paper.
Generation of Synthetic Data. We generate synthetic data using Lissajous curves [25]. More precisely, we simulate the current signal I(t) as the Lissajous curve l(t) which is defined by l(t) = (sin(2π m 1 t), sin(2π m 2 t)). (4.7) Here m 1 , m 2 are integer frequencies which implies that the trajectory is closed. Then, we assume for simplicity, that the sensitivity pattern of the drive coils P is the identity matrix, i.e. the trajectory of our non-dimensional field-free point r(t) is thus the Lissajous curve. We note that the described framework does not depend on this particular choice of trajectories. Experimental Setup. For our experiment we used a square N × N = 100 × 100 grid, as the underlying spatial signal -the ground truth-we use Figure 4.1 (d).
To simulate a scan we use the above Lissajous curves with frequencies m 1 = 101 and m 2 = 102. We pick K = 20 · N 2 = 200000 equidistant time samples t k from the time interval [0, 1] to obtain spatial points r k = r(t k ) with tangents v k = v(t k ) on the Lissajous curve.
Finally, we simulate the MPI time signal: we start with generating a signal according to the model where we evaluate A h [ρ](r k ) discretely by approximating the matrix-valued convolution via the summed-up midpoint rule on our 100 × 100 grid. The resolution parameter h of the convolution kernel is set to h = 10 −2 and is in the range which was identified in Section 2. Next, we add noise to imitate the noise-contamination of signal data from a MPI scanner The N k are i.i.d normally distributed with zero mean and a standard deviation of one. The noise amplitude ε is set to 10 percent of the signal strength: (4.10) The perturbed signalŝ k is now the input time data for our reconstruction algorithm. With this experimental setup we obtain the trace data u i from the signalŝ k . The last step of the algorithm is regularized deconvolution. The convolution kernel κ h here has the same value for h = 10 −2 as used in the setup. The Tychonov regularization-parameter µ was set to µ = 3 · 10 −4 . Moreover, we used the CG-method from the Octave library with a relative tolerance τ of τ = 2 10 −3 to solve (4.6) for ρ. The last two parameters µ and τ have to be adjusted in general. With our choice we obtained good results for our test case and the CG-method converges within 29 iterations.
Finally, the experiments were conducted on a laptop with a Intel Core I5-3337U CPU, 1.8 GHz, and 8GB of RAM. The programs were implemented in Octave and run on Octave 3.8.1 under Ubuntu 14.04.
Discussion of the Results. The result of our method is shown Figure 4.1. The generated noisy time signalŝ is displayed Figure 4.1(a). By (4.3) we obtain the spatial data u i = α h [ρ](x i ) out of the time signal; u is shown in Figure 4.1(b). The spatial data u is the right-hand side for the deconvolution problem. According to the theory developed in Section 3 the data u is a smoothed version of the particle density ρ in the noiseless case. Here the data u inherits its noise from the noise in the time signalŝ. By solving the Euler-Lagrange equation (4.6) with a CG-iteration we obtain the reconstruction of ρ displayed in Figure 4.1(c). Finally, Figure 4.1(d) shows the ground truth, the true particle density ρ, of our experiment.
Our result is a very good approximation of the ground truth and demonstrates the potential of our approach. Besides smoothing effects due to classical Tychonov regularization (which is a wellknown effect and not particular to MPI) the original density function ρ is well-reconstructed despite the relatively high noise level. Model-based reconstruction for multivariate MPI. Algorithm 1, proposed in this paper, is used to reconstruct a spatial particle density function ρ from a noisy MPI time signal. (a) shows the nosiy time signal s k ∈ R 2 (first component in red, second component in green). In (b) we see an intermediate step obtained from (4.3): this is the data u for our deconvolution problem. In a noiseless scenario the data u would be a smoothed version of the true particle density function ρ. Image (c) shows the reconstructed ρ obtained from deconvolution with the data u. The underlying ground truth, the true particle density function ρ, is displayed in (d). We observe that, besides smoothing effects due to classical Tychonov regularization (which is a well-known effect), the ground truth (d) is well-reconstructed by (c) despite the relatively high noise level.