Douglas-Rachford algorithm for magnetorelaxometry imaging using random and deterministic activations

Abstract

Magnetorelaxometry imaging is a novel tool for quantitative determination of the spatial distribution of magnetic nanoparticle inside an organism. The use of multiple excitation patterns has been demonstrated to significantly improve spatial resolution. However, increasing the number of excitation patterns is considerably more time consuming, because several sequential measurements have to be performed.

In this paper, we use compressed sensing in combination with sparse recovery to reduce the total measurement time and to improve spatial resolution. For image reconstruction, we propose using the Douglas-Rachford splitting algorithm applied to the sparse Tikhonov functional including a positivity constraint. Our numerical experiments demonstrate that the resulting algorithm is capable to accurately recover the magnetic nanoparticle distribution from a small number of activation patterns. For example, our algorithm applied with 10 activations yields half the reconstruction error of quadratic Tikhonov regularization applied with 50 activations, for a tumor-like phantom.

Keywords

Compressed sensing magnetorelaxometry image reconstruction magnetic nanoparticles multiple excitation Douglas-Rachford splitting sparse recovery

1. Introduction

Magnetic nanoparticles (MNP) offer a variety of promising biomedical applications. For example, they can be used as agents for drug delivery or hyperthermia, where the aim is to heat up specific regions inside a biological specimen [32]. These applications require quantitative knowledge of the magnetic nanoparticles distribution for safety and efficacy monitoring. In this paper, we consider magnetorelaxometry (MRX) imaging, which is a novel and promising non-invasive technique to spatially resolve the location and concentration of MNPs in vivo [4,23]. It beneficially combines a highly sensitive magnetic measurement technology for MNP imaging with a broad range of parameters and has the potential to image particle distributions in a comparably large volume [9].

In MRX imaging, the magnetic moments of the magnetic nanoparticles are aligned by a magnetic field generated by excitation coils [24]. Therewith, a magnetization can be measured from the MNP. After switching off the excitation field, the decay of magnetization (relaxation) is recorded, typically by SQUIDs (superconducting quantum interference device), yielding information about the particle concentration and related properties. If the relaxation is measured by a sensor array, the particle distribution can be reconstructed by inverse imaging methods [3,23].

In multiple excitation MRX (ME-MRX) imaging, several inhomogeneous activation fields generated by an array of excitation coils are applied. With such a coil array, a variety of inhomogeneous excitation fields can be generated, for example, by switching on the coils in a sequential manner [29]. First experimental realizations of ME-MRX with sequential activation and least squares estimation for imaging have been obtained in [23]. ME-MRX imaging using sequential coil activation and standard image reconstruction techniques requires a large number of measurement cycles and thus a considerable time for data acquisition. Recently, different approaches have been proposed using advanced excitation schemes [1,8,11].

A strategy to maintain high resolution while reducing measurement time is via compressed sensing (CS), a new sensing paradigm [5,12] that allows to capture a high-resolution image (or signal) by using fewer measurements than predicted by Shannon’s sampling theorem. CS replaces point-measurements by general linear measurements, where each measurement consists of a linear combination of the entries of the image of interest. Without prior information, recovering the original image is highly under-determined. The standard theory of CS predicts that under suitable assumptions on the image (sparsity) and the measurement matrix (incoherence), stable image reconstruction is nevertheless possible. CS has led to several new sampling strategies in medical imaging, for example, for speeding up MRI data acquisition [25], accelerating photoacoustic tomography [18], or completing under-sampled CT images [7]. While these applications involve relatively mildly ill-conditioned problems, ME-MRX constitutes a severely ill-conditioned or ill-posed problem [2,13]. No standard approaches for image reconstruction combining compressive measurements and such severely ill-posed problems exist.

In this work, we investigate CS techniques for accelerating ME-MRX. For that purpose, we consider random activations of the coils as well as sparse deterministic activation schemes. In order to image the MNP distribution from the CS data we propose a sparse reconstruction framework. We propose ℓ ¹-Tikhonov regularization together with a positivity constraint on the set of reconstructed MNP distributions, and apply the Douglas-Rachford splitting algorithm [10] for its minimization. We evaluate the performance of the resulting sparse recovery algorithm using different phantoms in dependence of the number of activation patterns. In any of the cases, the sparse recovery algorithm reduces the root mean square error (RSME) by more than 50%. For two of the three phantoms (CS and smiley), the deterministic scheme outperforms the random pattern by around 10% in terms of reduced RSME. For the third phantom (tumor), all tested sampling patterns roughly give the same results. We therefore can propose the sparse sampling pattern among the tested schemes, in combination with the proposed sparse recovery framework.

2. The forward model in ME-MRX imaging

This section gives a precise formulation of the forward problem of ME-MRX imaging and CS. Throughout this text we work with a discrete model [4,27,31]. A continuous model has recently been presented in [13].

Fig. 1.

Illustration of a simplified two-dimensional ME-MRX setup consisting of an array of activation coils and two layers of sensors arranged around the volume of interest. The coil locations r _s are indicated by red circles, and the corresponding normal vectors $\boldsymbol{\nu}_s$ by black arrows.

2.1. Signal generation

In MRX imaging, an activation field is applied to a region of interest containing MNPs. We assume this region to be divided into a number of N _v quadratic voxels, each represented by its midpoint located at r _v and containing a concentration n (r _v) of MNPs (compare Fig. 1). Due to the presence of the applied field, the MNPs align up and, after the activation field is switched off, produce a relaxation signal. According to the Biot-Savart law, the overall contribution of the MNP concentrations in the individual voxels to the measured signal at sensor location r _s is given by [23] $\begin{eqnarray}\displaystyle b^{\mathbf{H}}(\mathbf{r}_{s})=\displaystyle \frac{{\mu}_{0}}{4{\pi}}∼\boldsymbol{{\nu}}_{s}\boldsymbol{\cdot }\mathop{\sum }_{v=1}^{N_{v}}\left(\displaystyle \frac{3\mathbf{r}_{s,v}\otimes \mathbf{r}_{s,v}}{|\mathbf{r}_{s,v}|^{5}}-\displaystyle \frac{1}{|\mathbf{r}_{s,v}|^{3}}\right)∼\mathbf{H}(\mathbf{r}_{v})∼n(\mathbf{r}_{v})\quad \text{for }s\in \{1,\ldots ,N_{s}\}. & & \displaystyle\end{eqnarray}$ (1) Here H is the activation field, r _{s, v} = r _s − r _v is the vector joining r _s and r _v, $\boldsymbol{\nu}_s$ is the normal vector of the sensor at location r _s, and ⊗ and ⋅ are used to denote the tensor and inner product, respectively. Recall that the tensor product of two column vectors $\mathbf{a},\mathbf{b}\in \mathbb{R}^{3}$ is defined by $\mathbf{a}\otimes \mathbf{b}\triangleq \mathbf{a}\mathbf{b}^{\mathsf{T}}\in \mathbb{R}^{3\times 3}$ which results in a matrix, whereas the inner product $\mathbf{a}\boldsymbol{\cdot }\mathbf{b}\triangleq \mathbf{a}^{\mathsf{T}}\mathbf{b}\in \mathbb{R}$ results in a scalar.

Assuming the concentrations n (r _v) and the activation field H to be known, the measured data can be computed by (1). Collecting the individual measurements in a vector b ^H = (b ^H(r ₁), …, b ^H(r _{N
_s}))^T, we obtain the linear equation $\begin{eqnarray}\displaystyle \mathbf{b}^{\mathbf{H}}={\mathcal{L}}^{\mathbf{H}}∼\mathbf{n}\in \mathbb{R}^{N_{s}}. & & \displaystyle\end{eqnarray}$ (2) Here $\mathbf{n}=(n(\mathbf{r}_{1}),\ldots ,n(\mathbf{r}_{N_{v}}))^{\mathsf{T}}\in \mathbb{R}^{N_{v}}$ is the vector of MNP concentrations, and ${\mathcal{L}}\in \mathbb{R}^{N_{s}\times N_{v}}$ is the system matrix having entries $\begin{eqnarray}\displaystyle {\mathcal{L}}_{s,v}^{\mathbf{H}}\triangleq \displaystyle \frac{{\mu}_{0}}{4{\pi}}∼\boldsymbol{{\nu}}_{s}\boldsymbol{\cdot }\left(\displaystyle \frac{3\mathbf{r}_{s,v}\otimes \mathbf{r}_{s,v}}{|\mathbf{r}_{s,v}|^{5}}-\displaystyle \frac{1}{|\mathbf{r}_{s,v}|^{3}}\right)\mathbf{H}(\mathbf{r}_{v}) & & \displaystyle\end{eqnarray}$ (3) as derived from relation (1). The matrix ${\mathcal{L}}^{\mathbf{H}}$ is called Lead field matrix corresponding to the excitation field H. Equations (2), (3) constitute the standard discrete forward model of MRX using a single activation field.

Fig. 2.

Singular values of the Lead field matrix on a logarithmic scale. The Lead field matrix uses 110 × 120 = 13200 measurements for 75² = 5625 unknowns. The rapid decay of the singular values is due to the severe ill-conditioning of the MRX imaging problem. Increasing the number of measurements results in a similar rapid decay which indicates the inherent ill-posedness of the inverse problem.

2.2. Sequential coil activation

In ME-MRX imaging, the volume is (partially) surrounded by an array of excitation coils and sensor arrays (see Fig. 1). The coils can be steered to generate different types of activation fields. In the following we denote by H ₁, …, H _{N
_c} the magnetic fields induced by individual activation of the coils. Further, we write $\mathbf{b}_{c}\triangleq \mathbf{b}^{\mathbf{H}_{c}}\in \mathbb{R}^{N_{s}}$ and ${\mathcal{L}}_{c}\triangleq {\mathcal{L}}^{\mathbf{H}_{c}}\in \mathbb{R}^{N_{s}\times N_{v}}$ for the data and the Lead field matrix according to ((2)) and ((3)), corresponding to the activation field H = H _c. The measurements from all activations can be combined to a single equation of the form $\begin{eqnarray}\displaystyle \mathbf{b}={\mathcal{L}}∼\mathbf{n}\in \mathbb{R}^{N_{c}N_{s}}, & & \displaystyle\end{eqnarray}$ (4) where $\begin{eqnarray}\displaystyle \mathbf{b}\triangleq \left[\begin{array}{@{}c@{}}\mathbf{b}_{1}\\ \text{—}\\ \vdots \\ \text{—}\\ \mathbf{b}_{N_{c}}\end{array}\right]\in \mathbb{R}^{N_{c}N_{s}}\quad \text{ and }\quad ∼{\mathcal{L}}\triangleq \left[\begin{array}{@{}c@{}}{\mathcal{L}}_{1}\\ \text{—}\\ \vdots \\ \text{—}\\ {\mathcal{L}}_{N_{c}}\end{array}\right]\in \mathbb{R}^{N_{c}N_{s}\times N_{v}}. & & \displaystyle\end{eqnarray}$ (5) We will refer to b in ((4)) as full activation data. Evaluating ((4)) constitutes the forward model of ME-MRX. The corresponding inverse problem consists in determining the MNP distribution $\mathbf{n}\in \mathbb{R}^{N_{v}}$ from the data vector b that is additionally corrupted by noise. Note that Eq. ((4)) is known to be severely ill-conditioned as its singular values are rapidly decreasing (see [2]; an analysis in the infinite-dimensional setting has been performed in [13]). In Fig. 2 we plot the singular values of the forward matrix for full data. The singular values decay rapidly. This reflects the severe ill-conditioning and implies that only a small number of singular components can be robustly reconstructed.

The multiple coil setup has been realized in [1,29] and shown to significantly improve the spatial resolution compared to a single coil activation [29]. However, consecutive activations lead to a more time consuming measurement process. In order to accelerate data acquisition and to improve imaging resolution, we use CS techniques and develop a Douglas-Rachford based sparse reconstruction algorithm.

2.3. Compressive coil activations

The basic idea to employ CS for ME-MRX imaging is to use m specific coil activations, with m ≪ N _c, instead of activating all N _c coils in a sequential manner. Because of linearity, the data corresponding to the jth random activation pattern are given by $\begin{eqnarray}\displaystyle \mathbf{y}_{j}=\mathop{\sum }_{c=1}^{N_{c}}a_{j,c}\mathbf{b}_{c}=\mathop{\sum }_{c=1}^{N_{c}}a_{j,c}({\mathcal{L}}\mathbf{n})_{c}+\boldsymbol{{\xi}}_{j}\quad \text{ for }j=1,\ldots ,m. & & \displaystyle\end{eqnarray}$ (6) Here 𝜉 _j models the error in the data and a _j,c is the contribution of coil N _c to the jth activation pattern. The measurement matrix $\begin{eqnarray}\displaystyle \mathbf{A}\triangleq \left(\begin{array}{@{}ccc@{}}a_{1,1} & \cdots ∼ & a_{1,N_{c}}\\ \vdots & & \vdots \\ a_{m,1} & \ldots ∼ & a_{m,N_{c}}\end{array}\right)\in \mathbb{R}^{m\times N_{c}} & & \displaystyle\end{eqnarray}$ (7) represents all activation patters. Typical choices for A are Bernoulli or Gaussian matrices, since they are known to guarantee stable recovery of sparse signals [14]. Such kind of measurements can be naturally realized for ME-MRX, by simultaneous activation of several coils. The corresponding image reconstruction problem consists in recovering the MNP distribution n from the data in ((6)).

2.4. Matrix formulation of the reconstruction problem

In order to write the reconstruction problem in a compact form we introduce some additional notation. First, we define the CS measurement vector as $\begin{eqnarray}\displaystyle \mathbf{y}\triangleq \left[\begin{array}{@{}c@{}}\mathbf{y}_{1}\\ \text{—}\\ \vdots \\ \text{—}\\ \mathbf{y}_{m}\end{array}\right]\in \mathbb{R}^{mN_{s}}. & & \displaystyle\end{eqnarray}$ (8) Second, we introduce the vectorization or reshaping operator vec_{N
_c}, that takes a matrix to a vector whose block entries are equal to the transposes of the rows of the matrix: $\begin{eqnarray}\displaystyle \text{vec}_{N_{c}}:\mathbb{R}^{N_{c}\times N_{s}}\rightarrow \mathbb{R}^{N_{c}N_{s}}:\mathbf{X}\mapsto \left[\begin{array}{@{}c@{}}\mathbf{X}_{1,-}^{\mathsf{T}}\\ \text{—}\\ \vdots \\ \text{—}\\ \mathbf{X}_{N_{c},-}^{\mathsf{T}}\end{array}\right]. & & \displaystyle\end{eqnarray}$ (9) We denote by mat_{N
_c} the inverse reshaping operation that maps a vector in $\mathbb{R}^{N_{c}N_{s}}$ to a matrix in $\mathbb{R}^{N_{c}\times N_{s}}$ . Further, we write Y ≜ mat_m y and B ≜ mat_{N
_c} b.

With the above notations, the CS data in (6) can be written in the compact form $\begin{eqnarray}\displaystyle \mathbf{Y }=\mathbf{A}∼\mathbf{L}∼\mathbf{n}+\boldsymbol{{\Xi}}\in \mathbb{R}^{m\times N_{s}}, & & \displaystyle\end{eqnarray}$ (10) where L is the reshaped Lead field matrix defined by the property $\mathbf{L}\mathbf{n}=\text{mat}_{N_{c}}({\mathcal{L}}\mathbf{n})$ and 𝛯 denotes the noise matrix. With the Kronecker (or tensor) product ${\mathcal{A}}\triangleq \mathbf{A}\otimes {\mathcal{I}}\in \mathbb{R}^{(mN_{s})\times (N_{c}N_{s})}$ between the CS activation matrix A and the identity matrix ${\mathcal{I}}\in \mathbb{R}^{N_{s}\times N_{s}}$ , equation (10) can further be rewritten in the form $\begin{eqnarray}\displaystyle \mathbf{y}={\mathcal{A}}∼{\mathcal{L}}∼\mathbf{n}+\boldsymbol{{\xi}}\in \mathbb{R}^{mN_{s}}. & & \displaystyle\end{eqnarray}$ (11) The image reconstruction task in CS ME-MRX consists in recovering the MNP distribution n from Eq. (10) or, equivalently, from Eq. (11).

3. Douglas-Rachford algorithm for CS MRX imaging

In order to recover the MNP distribution, the inverse problem ((10)) (or ((11))) has to be solved which is known to be severely ill-posed, see Fig. 2. We note that no standard approach for image reconstruction combining compressive measurements and severely ill-posed problems exist. This section, therefore, is also of interest from a general perspective on CS for inverse problems.

3.1. Sparse Tikhonov regularization

We consider ((11)) as a single inverse problem with system matrix ${\mathcal{M}}={\mathcal{A}}{\mathcal{L}}$ . To address instability and non-uniqueness, and to incorporate prior knowledge, we use a sparse regularization approach. For that purpose, we minimize the generalized Tikhonov functional [28] $\begin{eqnarray}\displaystyle {\mathcal{T}}_{{\alpha},{\beta}}^{{\mathcal{M}}}(\mathbf{n})\triangleq \displaystyle \frac{1}{2}\Vert \mathbf{y}-{\mathcal{M}}\mathbf{n}\Vert _{2}^{2}+{\alpha}\Vert \boldsymbol{{\rm\Psi}}\mathbf{n}\Vert _{1}+{\beta}{\mathcal{R}}(\mathbf{n}). & & \displaystyle\end{eqnarray}$ (12) Here $\boldsymbol{{\rm\Psi}}:\mathbb{R}^{N_{v}}\rightarrow \mathbb{R}^{d}$ is a transform that sparsifies the MNP concentration (with $d\in \mathbb{N}$ denoting the dimension of the transformed domain), ${\mathcal{R}}:\mathbb{R}^{N_{v}}\rightarrow [0,\infty )$ is an additional regularizer that incorporates additional prior knowledge about the MNP distribution (such as positivity and other convex constraints) and ∥ ⋅ ∥_p is the standard ℓ ^p-norm defined by $\begin{eqnarray}\displaystyle \Vert \mathbf{x}\Vert _{p}\triangleq \sqrt[p]{\mathop{\sum }_{i=1}^{d}|x_{i}|^{p}}\quad \text{ for }\mathbf{x}=(x_{1},\ldots ,x_{d})\in \mathbb{R}^{d}. & & \displaystyle \nonumber\end{eqnarray}$ The non-negative parameters 𝛼 and 𝛽 allow to balance between the data consistency term $\frac{1}{2}\Vert \mathbf{y}-{\mathcal{M}}\mathbf{n}\Vert _{2}^{2}$ , the sparsity prior ∥Ψn∥₁ and the additional regularizer ${\mathcal{R}}(\mathbf{n})$ .

We call Ψ a sparsifying transform if Ψn can well be approximated by k-sparse vectors x, defined by the property that $\{i\mid x_{i}\neq 0\}$ has at most k elements. This approximation be can be quantified using the k-term approximation error, $\begin{eqnarray}\displaystyle {\sigma}_{k}(\mathbf{n})\triangleq \inf \{\Vert \mathbf{n}-\mathbf{n}_{k}\Vert _{1}\mid \mathbf{n}_{k}\text{ is }k\text{-sparse}\}. & & \displaystyle\end{eqnarray}$ (13) The k-term approximation error σ_k(n) appears as additive term in standard error estimates in compressed sensing theory. For a mathematically precise discussion of the k-term approximation error we refer to [14].

Remark 1 (Recovery theory for (12)).

The k-restricted isometry property (k-RIP) of ${\mathcal{M}}={\mathcal{A}}{\mathcal{L}}$ (after appropriate scaling) is defined as the smallest number δ_k such that for all k-sparse vectors Ψn we have $\begin{eqnarray}\displaystyle (1-{\delta}_{k})\Vert \mathbf{n}\Vert _{2}^{2}\leq \Vert {\mathcal{M}}\mathbf{n}\Vert _{2}^{2}\leq (1+{\delta}_{k})\Vert \mathbf{n}\Vert _{2}^{2}. & & \displaystyle\end{eqnarray}$ (14) Roughly spoken, CS theory [5,14] predicts uniform stable recovery with ((12)) for 𝛽 = 0, in the sense that it stably recovers any s-sparse vector provided that the s-RIP constant of ${\mathcal{M}}$ is sufficiently small.

Due to the severe ill-conditioning of the Lead field matrix ${\mathcal{L}}$ , the s-RIP is expected to be at most satisfied when the number N _v of voxels is small. Figure 2 shows that the singular values of ${\mathcal{L}}$ with N _v = 5625 are rapidly decaying and therefore the estimate $(1-{\delta}_{k})\Vert \mathbf{n}\Vert _{2}^{2}\leq \Vert {\mathcal{M}}\mathbf{n}\Vert _{2}^{2}$ for reasonable δ_k can only be satisfied for signals mainly formed by the first singular vectors. This implies that even in the case of full measurement data, where ${\mathcal{A}}$ is the identity matrix, uniform stable reconstruction of all sparse vectors is impossible. To obtain quantitative error estimates for (12), stable recovery results for individual elements [15–17] are a promising alternative. Such theoretical investigations are beyond the scope of this paper and an interesting line of future research.

Standard CS theory requires orthogonality of the sparsifying transform, which means that $\boldsymbol{Ψ}^{\mathsf{T}}\boldsymbol{Ψ}$ equals the identity matrix. Total variation (TV)-regularization [28] uses the gradient as sparsifying transform which is non-orthogonal. As in other compressed sensing imaging applications [21,26], we observed TV to outperform orthogonal sparsifying transforms. We will therefore focus on TV and, due to space limitations, no results for orthogonal bases like wavelets are included.

3.2. Douglas-Rachford algorithm

In order to minimize ((12)), we propose using the Douglas-Rachford minimization algorithm, which is a backward-backward type splitting method for minimizing the sum F + G of two functionals F and G. For our purpose we take $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle F(\mathbf{n})\triangleq \frac{1}{2}\Vert \mathbf{y}-{\mathcal{M}}\mathbf{n}\Vert _{2}^{2}\end{eqnarray}$ (15) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle G(\mathbf{n})\triangleq {\alpha}\Vert \boldsymbol{{\rm\Psi}}\mathbf{n}\Vert _{1}+{\beta}{\mathcal{R}}(\mathbf{n}).\end{eqnarray}$ (16) Minimizing ((12)) by the Douglas-Rachford algorithm, generates a sequence $(\mathbf{n}_{k})_{k\in \mathbb{N}}$ of estimated MNP distributions and auxiliary sequences $(\mathbf{z}_{k})_{k\in \mathbb{N}}$ , $(\tilde{\mathbf{z}}_{k})_{k\in \mathbb{N}}$ as described in Algorithm A. See Section 4.1, for a discussion on the choice of the parameter 𝜇.

The Douglas-Rachford algorithm is a splitting type algorithm for minimizing F + G which alternatingly performs updates according for F (in our case, the residual functional or data fitting term) and G (in our case, the regularizer). Because the regularizing term is non-smooth, standard methods such as Newton type methods cannot be applied. Splitting type methods are a natural choice in this case, and the implicit step in line 5 in Algorithm A accounts for the non-smoothness of G. Other splitting algorithms that are well suited to treat the non-smooth regularizer include the forward-backward splitting algorithm [10] and the Chambolle-Pock algorithm [6]. These algorithms use explicit gradient updates for the data fitting term requiring a matrix inversion. Therefore, they are not well suited for ME-MRX imaging. The Douglas-Rachford algorithm performs an implicit update for the data fitting term (line 4 Algorithm A), which is recognized as quadratic Tikhonov regularization step for the given equation and potentially accelerates the iteration.

For our numerical experiments, we take Ψ = 𝛻 as the discrete gradient operator as appropriate sparsifying transform for piecewise smooth MNP distributions. The additional regularizer is taken as the indicator function ${\mathcal{R}}=I_{C}$ of the convex set $\begin{eqnarray}\displaystyle C\triangleq \{\mathbf{n}\mid \forall \mathbf{r}_{v}:0\leq n(\mathbf{r}_{v})\leq n_{\text{ max}}\}\quad \text{for some bound }n_{\text{max}}, & & \displaystyle \nonumber\end{eqnarray}$ defined by I _C(n) = 0 if n ∈ C and I _C(n) = ∞. It guarantees non-negativity and boundedness. With the above choices, the Tikhonov functional (12) reduces to total variation (TV) regularization with positivity constraint. The minimization of $\frac{{\mu}}{2}\Vert \mathbf{n}-\mathbf{z}\Vert _{2}^{2}+\Vert {\nabla}\mathbf{z}\Vert _{1}+i_{C}(\mathbf{z})$ required for implementing Algorithm A is a TV denoising step and again performed by the Douglas-Rachford algorithm using the decomposition in $F(\mathbf{z})=\frac{{\mu}}{2}\Vert \mathbf{n}-\mathbf{z}\Vert _{2}^{2}+\Vert {\nabla}\mathbf{z}\Vert _{1}$ and G (z) = i _C(z).

Under the reasonable assumptions that the regularizer G is lower semicontinuous and convex, and that the sparse Tikhonov functional ${\mathcal{T}}_{{\alpha},{\beta}}^{{\mathcal{M}}}$ is coercive, the sequence $(\mathbf{n}_k)_{k\in \mathbb{N}}$ generated by Algorithm A is known to converge to a minimizer of (12); see [10,30]. Recall that the functional G is called lower semicontinuous if G(n) ≤ lim_k→∞ G(n _k) for any n and any sequence $(\mathbf{n}_k)_{k\in \mathbb{N}}$ converging to n. Note that Algorithm A performs implicit steps with respect to the residual functional $F=\frac{1}{2}\Vert {\mathcal{M}}\mathbf{n}-\mathbf{y}\Vert _{2}^{2}$ , which we found to have much faster convergence than forward-backward splitting algorithm [10] $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \mathbf{z}_{k}\triangleq \mathbf{n}_{k}-{\mathcal{M}}^{\mathsf{T}}({\mathcal{M}}\mathbf{n}_{k}-\mathbf{y})\end{eqnarray}$ (17) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle mnp_{k}\triangleq \text{argmin}_{\mathbf{n}}\frac{{\mu}}{2}\Vert \mathbf{z}_{k}-\mathbf{n}\Vert _{2}^{2}+G(\mathbf{n})\end{eqnarray}$ (18) where z _k are auxiliary quantities and n _k the reconstructed MNPs.

Fig. 3.

Full data setup for the numerical simulations. The phantom is contained in a square region of interest, the sensors are arranged in two parallel layers (blue dots), and the coils for full activation are arranged around the phantom in U-form (red dots).

4. Numerical results

For the following numerical simulations, we use a data setup similar to the realization in [1]. For simplification, we consider a two-dimensional setup representing one voxel plane, containing two parallel arrays of detector elements (one measuring the horizontal (1,0) and one measuring the vertical (0,1) component) located outside a quadratic region of interest. Circular shaped activation coils are arranged in U-form around the region of interest, all having the same normal vector (0,1). The used arrangement of sensor and coil locations for the full measurement data setup is shown in Fig. 3. For our simulations, we choose a discretization of imaging space into N _v = 75² voxels covering a region of interest of [−5,5] × [−5,5] cm². The data are generated for N _s = 110 positions and full measurement data correspond to N _c = 120 activation coils outside the region of interest.

Fig. 4.

Phantoms, data and full measurement data reconstruction. Top row: CS-phantom (left), smiley-phantom (middle) and tumor-phantom mimicking cancerous tissue with included blood vessel (right). Middle row: Corresponding full measurement data. Bottom row: Reconstructions from full measurement data using standard quadratic Tikhonov regularization.

4.1. Forward computations

To set up the forward model ((2)), ((4)), ((5)), we first have to compute the magnetic fields generated by each coil. For that purpose, we follow the approach of [1,3], where any circular activation coil is approximated by short line segments, each carrying the same current I ₀. Every coil is modeled by 45 line segments for numerical calculations of the magnetic field. This is considered to be a sufficiently accurate approximation since the deviation of only 36 line segments from the analytical solution was shown to be already below 1% in [22]. Using this approximation, the induced magnetic field at the voxel center r _v can be computed by (see [20]) $\begin{eqnarray}\displaystyle \mathbf{H}(\mathbf{r}_{v})\simeq \displaystyle \frac{1}{4{\pi}}\mathop{\sum }_{i=1}^{45}\frac{|\mathbf{r}_{1,v,i}|+|\mathbf{r}_{2,v,i}|}{|\mathbf{r}_{1,v,i}|∼|\mathbf{r}_{2,v,i}|}\cdot \displaystyle \frac{\mathbf{r}_{1,v,i}\times \mathbf{r}_{2,v,i}}{|\mathbf{r}_{1,v,i}|∼|\mathbf{r}_{2,v,i}|+\mathbf{r}_{1,v,i}\boldsymbol{\cdot }\mathbf{r}_{2,v,i}}∼I_{0}, & & \displaystyle\end{eqnarray}$ (19) where r _1,v,i and r _2,v,i are the distance vectors between the voxel center r _v and the beginning and end points of the ith line segment, respectively. For the presented numerical computations, we use a coil diameter of 1 μm, illustrating an almost point like coil. Having computed the activation field H(r _v), we compute the entries of the Lead field matrix according to ((3)). By activating the coils sequentially, we obtain full measurement data and the full Lead field matrix.

Fig. 5.

Selection of the regularization parameter by the L-curve method. Each image shows a log–log plot of the residual functionals versus the solution norm in dependence of the regularization parameter. The L-curve method predicts an L-shaped curve and advices to select a regularization parameter close to the corner of the L-curve. In our results we empirically found that taking the regularization parameter slightly smaller than advised by the L-curve method yields accurate results.

For the presented results, we use three different magnetic particle distributions which, together with the corresponding measurements full data, are shown in Fig. 4. Any column in the measurement data (in Fig. 4 and below) corresponds to a single activation pattern and contains the data of all detectors. The phantoms are rescaled to have maximum value 1 and the forward operator ${\mathcal{L}}$ is rescaled to have matrix norm $\Vert {\mathcal{L}}\Vert _{2}=1$ . To all data we have added additive Gaussian noise amounting to a signal-to-noise (SNR) of 80 dB. The corresponding reconstructions $\begin{eqnarray}\displaystyle \mathbf{n}_{{\mu}}=({\mathcal{L}}^{\mathsf{T}}{\mathcal{L}}+{\mu}{\mathcal{I}})^{-1}{\mathcal{L}}^{\mathsf{T}}\mathbf{b} & & \displaystyle\end{eqnarray}$ (20) with standard quadratic Tikhonov regularization (penalized least squares) using regularization parameter 𝜇 = 10⁻¹² are shown in the bottom row of Fig. 4. Recall that ${\mathcal{L}}\in \mathbb{R}^{N_{c}N_{s}\times N_{v}}$ is the full Lead field matrix and $\mathbf{b}\in \mathbb{R}^{N_{c}N_{s}}$ the full activation data. Reconstructions with ((20)) can be seen as benchmark for the more sophisticated CS reconstructions applied to less data presented below. The regularization parameter has been chosen empirically as the trade-off between stability and smoothing. Taking a smaller regularization parameter would not suppress the noise (amplification) well enough, whereas a larger regularization parameter yields to overshooting. There are many strategies for selecting the regularization parameter automatically [19]. An example is the L-curve method, where the residual $\Vert {\mathcal{L}}\mathbf{n}_{{\mu}}-\mathbf{y}\Vert _{2}^{2}$ is plotted against the solution norm ∥n _𝜇∥ both on a logarithmic scale. The L-curve method predicts the graph to be an L-shaped curve and the regularization parameter should be taken close to the corner of the L-shaped curve. The L-curves for the considered phantoms are plotted in Fig. 5. Indeed, the selected 𝜇 =10⁻¹² are close to the corner in each case. Strictly taken, the L-curve method predicts a larger regularization which we however found to yield to a slight over-smoothing and a larger reconstruction error.

The CS forward matrix ${\mathcal{M}}={\mathcal{A}}{\mathcal{L}}$ is computed by multiplying the full Lead field matrix with ${\mathcal{A}}$ . CS measurements have been generated in two random ways and one deterministic way. In the random case, A is taken either as Bernoulli matrix having entries ±1 appearing with equal probability, or a Gaussian matrix consisting of i.i.d. ${\mathcal{N}}(0,1)$ -Gaussian random variables in each entry. The deterministic sparse sampling is performed by choosing m equispaced coil activations.

Fig. 6.

Reconstruction results with Algorithm A from 40 coil activations. Top row: Gaussian activation matrix. Middle row: Bernoulli activation matrix. Bottom row: Sparse deterministic activation pattern.

Table 1

Evaluation metrics for Gaussian measurements (row 1) Bernoulli measurements (row 2) and deterministic activation pattern (row 3) using 40 activation patterns. The bottom row shows the same evaluation metrics for Tikhonov regularization for full activations

	Relative RMSE (%)	SNR (dB)	Correlation
40 Gaussian activations (proposed algorithm)
CS-phantom	0.75	2.51	0.92
smiley-phantom	0.40	8.05	0.87
tumor-phantom	0.16	15.69	0.98
40 Bernoulli activations (proposed algorithm)
CS-phantom	0.78	2.15	0.88
smiley-phantom	0.39	8.11	0.87
tumor-phantom	0.17	15.56	0.98
40 deterministic activations (proposed algorithm)
CS-phantom	0.66	3.59	0.92
smiley-phantom	0.35	9.16	0.90
tumor-phantom	0.17	15.38	0.98
120 activations (Tikhonov regularization)
CS-phantom	1.20	--1.59	0.69
smiley-phantom	0.47	6.45	0.80
tumor-phantom	0.35	8.99	0.89

4.2. Reconstruction results

The system matrix ${\mathcal{M}}$ is rescaled to have matrix norm 1 and we use the parameter setting 𝜇 = 4 ⋅10⁻¹³, 𝛼 = 10⁻¹⁴, s = 1, n _max = 1 and N _iter = 50. Note that the parameter 𝜇 has a similar role for each iterative step as in quadratic Tikhonov regularization (20). Choosing it too small causes noise amplification. Choosing it too large causes oversmoothing, which however will be reduced during the iteration. Therefore, taking it somewhat larger than the value found in Tikhonov case seems reasonable. The regularization parameter 𝛼 has been selected empirically to yield good numerical performance.

We observed that after 50 iterations the reconstructed MNP stagnates which indicates that the iterates are close to the minimizer of the Tikhonov functional. The reconstructed MNP distributions using Algorithm A from 40 coil activations are shown in Fig. 6. Each reconstruction takes about 50 seconds in Matlab R 2017a on a MacBook Pro (2016) with 2.9 GHz Intel Core i7 processor.

To quantitatively evaluate the reconstruction results, we compute the relative root mean squared error (RMSE) ∥n − n _rec∥∕∥n∥, the SNR 20lg(∥n∥∕∥n − n _rec∥), and the Pearson correlation coefficient cov(n, n _rec)/ $\sqrt{\mathrm{var}(\mathbf{n})\mathrm{var}(\mathbf{n}_{\mathrm{rec}})}$ with cov denoting the covariance and var the variance. The results are shown in Table 1. For all phantoms and for any evaluation metric, the proposed algorithm clearly outperforms quadratic Tikhonov regularization. More specifically, the RSME for the proposed algorithm applied with 40 activations yields half of the RSME compared to quadratic Tikhonov regularization with 120 activations. For the CS-phantom and the smiley-phantom, the deterministic scheme performs best and decreases the reconstruction error by more than 10% compared to the random schemes. For the tumor phantom, all sampling schemes yield comparable performance.

Fig. 7.

Quality evaluation in dependence of the number of coil activations. Left: correlation between reconstructions and the true phantoms. Right: RMSE of the reconstructions.

Figure 7 shows a convergence study (for the correlation and the RMSE) in dependence of the number of activations. We see that the Douglas-Rachford algorithm constantly outperforms Tikhonov regularization and significantly reduces the RSME and increases the correlation. For the CS and the smiley phantom, the deterministic sampling pattern has a significantly smaller RMSE than the random schemes. We address this behavior to the strong smoothing effect of the forward operator, which removes most high frequency components in the data. In the full forward matrix A L, the incoherence contained in the measurements matrix A is annihilated by the smoothing effect of A. Finding particular activation patterns that outperform the sparse sampling scheme is a nontrivial issue and will be investigated in future work.

5. Conclusion

The use of multiple coil activation patterns in magnetorelaxometry imaging is time consuming and requires performing several consecutive measurements. It is therefore desirable to make the number of coil activations as small as possible, while keeping high spatial resolution. For that purpose, we investigated CS strategies in this paper. We compared Gaussian random activations, Bernoulli activations and deterministic sparse sampling schemes. For actual image reconstruction, we proposed a Douglas-Rachford splitting algorithm applied to the sparse Tikhonov functional; see Algorithm A. For a small number of coil activations, the random schemes slightly outperform the deterministic scheme for the smiley-phantom and tumor phantom. For a large number of activations, the deterministic sampling scheme clearly performed better than the random sampling patterns. We explain this behavior by the severe ill-conditioning of the Lead field matrix, which is more severe for a large number of activations. Only for a small number of activations, where the resolution is poor anyway, the incoherence at the low frequencies is not destroyed by the smoothing effect of the Lead field matrix. The reconstruction results clearly demonstrate the proposed Algorithm A significantly outperforms standard algorithms such as Tikhonov regularization. Depending on the complexity of the phantom between 6 (tumor) and 30 activations (CS) are sufficient for the used framework and further increasing the number of activations only decreases the RMSE by a few percent; see Fig. 7.

Several interesting research directions following this work are possible. First, we can replace the inner TV minimization in the single-stage approach by a different algorithm which should accelerate Algorithm A. Second, the derivation of theoretical error estimates for the single-stage approach is of significant interest. Results in that direction can advise which type of MRX measurements are best to obtain accurate results for certain phantom classes. Moreover, the derivation of adaptive compressed sensing strategies for online monitoring is of significant interest. Advising optimal coil activations given previous activations is a practically important challenge that will benefit from theoretical error estimates, numerical simulations as well as real-world experiments. Such issues will be addressed in future work. In this paper, we investigated standard random compressed sensing schemes (using Gaussian and Bernoulli activation patterns and sparse sampling). Using more problem adapted and task oriented measurement design we expect to derive improved coil activation schemes.

Footnotes

Acknowledgements

This work has been supported by the Austrian Science Fund (FWF), project P 30747-N32, and the German Science Foundation (DFG) within the priority program COSIP, project Cos-MRXI (BA 4858/2-1).

References

Baumgarten

Braune

Supriyanto

and Haueisen

, Plane-wise sensitivity based inhomogeneous excitation fields for magnetorelaxometry imaging of magnetic nanoparticles, J. Magn. Magn. Mater. 380 (2015), 255–260.

Baumgarten

Eichardt

Crevecoeur

Supriyanto

and Haueisen

, Magnetic nanoparticle imaging by random and maximum length sequences of inhomogeneous activation fields, in: 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), 2013, pp. 3258–3260.

Baumgarten

and Haueisen

, A spatio-temporal approach for the solution of the inverse problem in the reconstruction of magnetic nanoparticle distributions, IEEE Trans. Magn. 46(8) (2010), 3496–3499.

Baumgarten

Liehr

Wiekhorst

Steinhoff

Münster

Miethe

Trahms

and Haueisen

, Magnetic nanoparticle imaging by means of minimum norm estimates from remanence measurements, Med. Biol. Eng. Comput. 46(12) (2008), 1177–1185.

Candès

E.J.

Romberg

J.K.

and Tao

, Stable signal recovery from incomplete and inaccurate measurements, Comm. Pure Appl. Math. 59(8) (2006), 1207–1223.

Chambolle

and Pock

, A first-order primal-dual algorithm for convex problems with applications to imaging, J. Math. Imaging Vision 40(1) (2011), 120–145.

Chen

G.-H.

Tang

and Leng

, Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Med. Phys. 35(2) (2008), 660–663.

Coene

Crevecoeur

and Dupre

, Adaptive control of excitation coil arrays for targeted magnetic nanoparticle reconstruction using magnetorelaxometry, IEEE Trans. Magn. 48(11) (2012), 2842–2845.

Coene

Leliaert

Liebl

Loewa

Steinhoff

Crevecoeur

Dupré

and Wiekhorst

, Multi-color magnetic nanoparticle imaging using magnetorelaxometry, Phys. Med. Biol. 62(8) (2017), 3139.

10.

Combettes

P.L.

and Pesquet

J.-C.

, Proximal splitting methods in signal processing, in: Fixed-Point Algorithms for Inverse Problems in Science and Engineering, Springer, 2011, pp. 185–212.

11.

Crevecoeur

Baumgarten

Steinhoff

Haueisen

Trahms

and Dupré

, Advancements in magnetic nanoparticle reconstruction using sequential activation of excitation coil arrays using magnetorelaxometry, IEEE Trans. Magn. 48(4) (2012), 1313–1316.

12.

Donoho

D.L.

, Compressed sensing, IEEE Trans. Inf. Theory 52(4) (2006), 1289–1306.

13.

Föcke

Baumgarten

and Burger

, The inverse problem of magnetorelaxometry imaging, Inverse Problems 34(11) (2018), 115008.

14.

Foucart

and Rauhut

, A Mathematical Introduction to Compressive Sensing, Birkhäuser, Basel, 2013.

15.

Fuchs

J.J.

, Recovery of exact sparse representations in the presence of bounded noise, IEEE Trans. Inf. Theory 51(10) (2005), 3601–3608.

16.

Grasmair

Haltmeier

and Scherzer

, Necessary and sufficient conditions for linear convergence of ℓ ¹-regularization, Comm. Pure Appl. Math. 64(2) (2011), 161–182.

17.

Haltmeier

, Stable signal reconstruction via ℓ ¹-minimization in redundant, non-tight frames, IEEE Trans. Signal Process. 61(2) (2013), 420–426.

18.

Haltmeier

Berer

Moon

and Burgholzer

, Compressed sensing and sparsity in photoacoustic tomography, J. Opt. 18(11) (2016), 114004 (12 pages).

19.

Hansen

P.C.

, Rank-deficient and discrete Ill-posed problems, in: SIAM Monographs on Mathematical Modeling and Computation, SIAM, Philadelphia, PA, 1998.

20.

Hanson

J.D.

and Hirshman

S.P.

, Compact expressions for the Biot–Savart fields of a filamentary segment, Phys. Plasmas 9(10) (2002), 4410–4412.

21.

Krahmer

and Ward

, Stable and robust sampling strategies for compressive imaging, IEEE Trans. Image Process. 23(2) (2014), 612–622.

22.

Liebl

, Bildgebung von magnetischen Nanopartikelverteilungen mit einem zylindrischen magnetisch geschirmten squid-system, 2011. Master's thesis, Ilmenau University of Technology, Germany.

23.

Liebl

Steinhoff

Wiekhorst

Haueisen

and Trahms

, Quantitative imaging of magnetic nanoparticles by magnetorelaxometry with multiple excitation coils, Phys. Med. Biol. 59(21) (2014), 6607.

24.

Liebl

Wiekhorst

Eberbeck

Radon

Gutkelch

Baumgarten

Steinhoff

and Trahms

, Magnetorelaxometry procedures for quantitative imaging and characterization of magnetic nanoparticles in biomedical applications, Biomed. Tech. 60(5) (2015), 427–443.

25.

Lustig

Donoho

D.L.

Santos

J.M.

and Pauly

J.M.

, Compressed sensing mri, IEEE Sig. Proc. Mag. 25(2) (2008), 72–82.

26.

Poon

, On the role of total variation in compressed sensing, SIAM J. Imaging Sci. 8(1) (2015), 682–720.

27.

Richter

Kettering

Wiekhorst

Steinhoff

Hilger

and Trahms

, Magnetorelaxometry for localization and quantification of magnetic nanoparticles for thermal ablation studies, Phys. Med. Biol. 55(3) (2010), 623–633.

28.

Scherzer

Grasmair

Grossauer

Haltmeier

and Lenzen

, Variational Methods in Imaging, Applied Mathematical Sciences Vol. 167, Springer, New York, 2009.

29.

Steinhoff

Wiekhorst

Baumgarten

Haueisen

and Trahms

, Imaging of magnetic nanoparticles based on magnetorelaxometry with sequential activation of inhomogeneous magnetization fields, Biomed. Tech. 55 (2010), 22–25.

30.

Svaiter

B.F.

, On weak convergence of the Douglas–Rachford method, SIAM J. Control Optim. 49(1) (2011), 280–287.

31.

Wiekhorst

Seliger

Jurgons

Steinhoff

Eberbeck

Trahms

and Alexiou

, Quantification of magnetic nanoparticles by magnetorelaxometry and comparison to histology after magnetic drug targeting, J. Nanosci. Nanotechnol. 6(9-10) (2006), 3222–3225.

32.

Wiekhorst

Steinhoff

Eberbeck

and Trahms

, Magnetorelaxometry assisting biomedical applications of magnetic nanoparticles, Pharm. Res. 29(5) (2012), 1189–1202.