Determining anisotropic conductivity using diffusion tensor imaging data in magneto-acoustic tomography with magnetic induction

TL;DR

This paper introduces a novel method combining magneto-acoustic tomography and diffusion tensor imaging data to accurately reconstruct anisotropic electrical conductivity tensors in tissues.

Contribution

It develops a mathematical and numerical framework for reconstructing anisotropic conductivity tensors using an optimal control approach integrating DTI data.

Findings

The proposed algorithm converges and is Lipschitz stable.

Numerical examples demonstrate the method's accuracy.

The approach improves spatial resolution in conductivity imaging.

Abstract

In this paper we present a mathematical and numerical framework for a procedure of imaging anisotropic electrical conductivity tensor by integrating magneto-acoutic tomography with data acquired from diffusion tensor imaging. Magneto-acoustic Tomography with Magnetic Induction (MAT-MI) is a hybrid, non-invasive medical imaging technique to produce conductivity images with improved spatial resolution and accuracy. Diffusion Tensor Imaging (DTI) is also a non- invasive technique for characterizing the diffusion properties of water molecules in tissues. We propose a model for anisotropic conductivity in which the conductivity is proportional to the diffusion tensor. Under this assumption, we propose an optimal control approach for reconstructing the anisotropic electrical conductivity tensor. We prove convergence and Lipschitz type stability of the algorithm and present numerical examples to illustrate its accuracy and feasibility.

Equations

$$\left\{\begin{array}[]{lcll}\nabla\times E_{\gamma}&=&B_{1}&\mbox{in }\Omega,\\ \nabla\cdot(\gamma E_{\gamma})&=&0&\mbox{on }\Omega,\\ \gamma E_{\gamma}\cdot\nu&=&0&\mbox{on }\partial\Omega.\end{array}\right.$$

$$\gamma (x)=\sigma (x)D(x).$$

$$c_1\le \sigma (x)\le c_2,\forall x\in \Omega ,$$

$$\lambda \Vert \xi {\Vert }_2^2\le {\xi }^{{}^{\prime }}D\xi \le \Vert \xi {\Vert }_2^2,\forall \xi \in R^3,$$

$$def:solution{eqn:N-bd}\left\{\begin{array}[]{lll}\nabla\cdot(\sigma D\nabla u)&=-\nabla\cdot{E},&\qquad\mbox{ in }\Omega,\\ (\sigma D\nabla u+{E})\cdot\nu&=0,&\qquad\mbox{ on }\partial\Omega,\end{array}\right.$$

$$def:solutiondef:solution{\int }_{\Omega }\sigma D\nabla u\cdot \nabla \phi \mathrm{d}x=-{\int }_{\Omega }E\cdot \nabla \phi \mathrm{d}x,\forall \phi \in H^1.$$

$$\Vert \nabla u{\Vert }_{L^2}\le c_1^{-1}{\lambda }^{-1}\Vert E{\Vert }_{L^2}.$$

$$c_1\lambda \Vert \nabla u{\Vert }_{L^2}^2\le {\int }_{\Omega }\sigma \nabla u\cdot D\nabla u\mathrm{d}x.$$

$${\int }_{\Omega }\sigma D\nabla u\cdot \nabla u\mathrm{d}x=-{\int }_{\Omega }E\cdot \nabla u\mathrm{d}x.$$

$$c_1\lambda \Vert \nabla u{\Vert }_{L^2}^2\le \left|-{\int }_{\Omega }E\cdot \nabla u\mathrm{d}x\right|\le \Vert \nabla u{\Vert }_{L^2}\Vert E{\Vert }_{L^2},$$

$$\Vert E_{\sigma D}{\Vert }_{L^2}\le C_1,$$

$$\Vert E_{\sigma D}{\Vert }_{L^{\infty }(\Omega )}\le C_2,$$

$$\Vert E_{\sigma D}{\Vert }_{C^{1,\gamma }(\Omega )}\le C_3.$$

$$\Vert E_{{\sigma }_{1}D}-E_{{\sigma }_{2}D}{\Vert }_{L^2(\Omega )}\le C\Vert {\sigma }_1-{\sigma }_2{\Vert }_{L^2(\Omega )}.$$

$$\Vert E_{{\sigma }_{1}D}-E_{{\sigma }_{2}D}{\Vert }_{H^1(\Omega )}\le C\Vert {\sigma }_1-{\sigma }_2{\Vert }_{H^1(\Omega )}.$$

$$\left\{\begin{array}[]{lll}\nabla\cdot(\sigma D\nabla u)&=-\nabla\cdot(\sigma D{\tilde{E}}),&\qquad\mbox{ in }\Omega,\\ (\sigma D\nabla u+\sigma D{\tilde{E}})\cdot\nu&=0,&\qquad\mbox{ on }\partial\Omega.\end{array}\right.$$

$$\Vert E_{\sigma D}{\Vert }_{L^2}\le C_1.$$

$$\nabla u=E_1-E_2.$$

$$\left\{\begin{array}[]{lll}\nabla\cdot(\sigma_{1}D\nabla u)&=-\nabla\cdot((\sigma_{1}-\sigma_{2})D{E}_{2}),&\qquad\mbox{ in }\Omega,\\ D\nabla u\cdot\nu&=0,&\qquad\mbox{ on }\partial\Omega.\end{array}\right.$$

$$\Vert \nabla u{\Vert }_{L^2}\le c_1^{-1}{\lambda }^{-1}\Vert ({\sigma }_1-{\sigma }_2)D{E}_2{\Vert }_{L^2}.$$

$$\Vert E_1-E_2{\Vert }_{L^2}\le C\Vert ({\sigma }_1-{\sigma }_2){\Vert }_{L^2}.$$

$$\Vert u{\Vert }_{H^2(\Omega )}\le C\Vert \nabla \cdot (({\sigma }_1-{\sigma }_2)D{E}_2){\Vert }_{L^2(\Omega )},$$

$$\Vert \nabla {\sigma }_i{\Vert }_{L^{\infty }}\le K$$

$$|1-\lambda |\le \eta$$

$$\Vert \nabla ({\sigma }_1-{\sigma }_2){\Vert }_{L^2(\Omega )}\le L\Vert {\sigma }_1-{\sigma }_2{\Vert }_{L^2(\Omega )},$$

$$c\Vert {\sigma }_1-{\sigma }_2{\Vert }_{L^2(\Omega )}\le \Vert F({\sigma }_1)-F({\sigma }_2){\Vert }_{L^2(\Omega )}$$

$$F({\sigma }_1)-F({\sigma }_2)$$

$$I_1$$

$$I_2$$

$$I_3$$

Full text

Keywords: hybrid imaging, magneto-acoustic tomography with magnetic induction, diffusion tensor imaging, anisotropic conductivity.

1 Introduction

In this paper, we describe a method of reconstructing images of an anisotropic conductivity tensor distribution inside an electrically conducting object using Magneto-acoustic Tomography with Magnetic Induction (MAT-MI). MAT-MI is a new noninvasive modality for imaging electrical conductivity distributions of biological tissue [25, 18]. In the experiments, the biological tissue is placed in a static magnetic field. A pulsed magnetic field is applied to induce an eddy current inside the conductive tissue. In the process, the tissue emits ultrasound waves which can be measured around the tissue. The first step in the MAT-MI is to recover the acoustic source in the scalar wave equation from measurements at a set of locations around the object. This problem has been studied in many works, including [4, 12, 15, 16, 21]. The second step in the MAT-MI is to reconstruct the electrical conductivity distribution from knowledge of the acoustic source.

Biological tissues are known to have anisotropic conductivity values [19, 24]. However, up to now, all MAT-MI techniques have been devised to image isotropic conductivity distribution. The fundamental questions arising in this case have been addressed in [1, 20]. However, for the case of anisotropic conductivity, many basic questions remain. For instance, it is not known how many MAT-MI measurements are needed to uniquely determine the conductivity tensor.

DTI is a non-invasive technique for characterizing the diffusion properties of water molecules in tissues (see e.g. [10] and the references therein). Imaging conductivity tensors in the tissue with DTI is based on the correlation property between diffusion and conductivity tensors [24]. This linear relationship can be used to characterize the conductivity tensor. Once the conductivity directions of anisotropy are determined, one needs only to reconstruct a cross-property factor which is a scalar function. In [7, 8], it is shown how to recover this factor in Current Density Impedance Imaging. In [3], a multifrequency electrical impedance approach is developed for estimating the ratio between the largest and the lowest eigenvalue of the electrical conductivity tensor. An iterative procedure for reconstructing anisotropic conductivities from internal current densities has been proposed in [22].

In the process of the MAT-MI experiment, the tissue is placed in a constant static magnetic background field $B_{0}=(0,0,1)$. A pulsed magnetic stimulation of the form $B_{1}u(t)$ is applied, where the vector field $B_{1}$ is constant and $u(t)$ is the time variation. Let $\gamma$ denote the conductivity, $E_{\gamma}$ denote electric field, and $\Omega$ be the domain to be imaged. Then the electric field satisfies the following Maxwell equations

[TABLE]

The second step of MAT-MI is to reconstruct $\gamma$ from the known internal data $\nabla\cdot(\gamma E_{\gamma}\times B_{0})$ on $\Omega$.

In this paper, we will consider the anisotropic conductivity $\gamma$. We assume that DTI has been performed on the tissue being imaged, that is, the diffusion tensor $D(x)$ has been found. Next we follow [24] and make the assumption that conductivity is proportional to $D(x)$ and assume that the conductivity tensor $\gamma$ is of the form

[TABLE]

The cross-property factor $\sigma$ is a scalar function to be reconstructed. We will focus on the second step of MAT-MI combined with DTI, i.e., on reconstructing the cross-property factor $\sigma$ from the internal data given by $\nabla\cdot(\gamma E_{\gamma}\times B_{0})$ with known conductivity tensor $D(x)$.

In the following, we assume that $D(x)$ is a positive definite symmetric matrix everywhere and write it as $D=T^{{}^{\prime}}\Sigma T$, where $D=\mathrm{diag}(e_{1},e_{2},e_{3})$, $e_{1}\geq e_{2}\geq e_{3}$ are the eigenvalues of $D(x)$. The columns of $T^{{}^{\prime}}$ are the corresponding eigenvectors. As we can always write $\sigma=\sigma_{0}e_{1}T^{{}^{\prime}}\mathrm{diag}(1,e_{2}/e_{1},e_{3}/e_{1})T$, we assume that $e_{1}=1$ hereinafter.

2 Notation and preliminaries

In this section, we introduce the notation for the mathematical analysis. Let $\Omega$ be a bounded Lipschitz domain in $\mathbb{R}^{3}$. A typical point $x=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}$ denotes the spatial variable. Throughout this paper, the standard notation for Hölder and Sobolev spaces and their norms is used. If there is no confusion, we omit the dependence on the domain.

Assumption 1**.**

Let $\sigma$ and $D$ be positive functions belonging to $C^{1,\beta}$, $\beta>0$ and assume that

[TABLE]

and

[TABLE]

for some constants $\lambda,c_{1},c_{2}>0$.

Here we state several useful results on elliptic partial differential equations with Neumann boundary conditions.

We say that $u\in H^{1}$ is a weak solution of the Neumann boundary value problem

[TABLE]

if

[TABLE]

We give a brief proof of the following regularity result and standard energy estimate.

Proposition 1**.**

Suppose that $\sigma$ and $D$ satisfy Assumption 1. For field ${E}\in L^{2}$, the Neumann problem (LABEL:eqn:N-bd) has a solution $u\in H^{1}$. The solution $u$ is unique up to an additive constant and satisfies the estimate

[TABLE]

Proof.

The proof of the existence and uniqueness up to an additive constant is a standard result by the Lax-Milgram Theorem. We refer the readers to [23]. In the following, we prove the gradient estimate (2.5).

It follows from the ellipticity condition (2.2) that

[TABLE]

Taking the test function $\varphi$ in Definition LABEL:def:solution to be the solution $u$, we have that

[TABLE]

Consequently, applying the Cauchy-Schwarz inequality, we obtain that

[TABLE]

and (2.5) follows. ∎

Proposition 2**.**

Let $\sigma$ and $D$ satisfy Assumption 1. Then the system (1.1) is uniquely solvable and there exists a constant $C$ and $C_{i}$ $(1\leq i\leq 3)$ depending on $\lambda$, $c_{1}$, $c_{2}$ and $\Omega$, such that

[TABLE]

[TABLE]

[TABLE]

Moreover, if $\sigma_{1}$ and $\sigma_{2}$ satisfy Assumption 1, we have the following bound for the electric field difference,

[TABLE]

[TABLE]

Proof.

Let us first reduce the system (1.1) to a Neumann boundary value problem. Let ${\tilde{E}}=\frac{1}{2}(-y,x,0)$. We can readily check that $\nabla\times{\tilde{E}}={B}_{1}$. Hence $\nabla\times({E}_{\sigma D}-{\tilde{E}})=0$ and we can write ${E}_{\sigma D}={\tilde{E}}+\nabla u$. Substituting this into (1.1), we have that $u$ solves the Neumann boundary value problem

[TABLE]

With the help of proposition 1, we get

[TABLE]

From the standard $L^{p}$ estimate for elliptic equations [5, Chapter 9] and the Sobolev Embedding Theorem, we know that ${E}_{\sigma D}$ is a bounded function in $W^{2,p}(\Omega)$ for any $p>2$. Hence, ${E}_{\sigma D}$ is uniformly bounded by a constant $C$, which depends only on $r_{0}$, $\lambda$, $c_{1}$, $c_{2}$, and $\Omega$. Then (2.7) is proved.

With the assumption of $C^{1,\gamma}$ property, we would obtain the $C^{2,\gamma}$ Hölder continuity[6] of $u$, i.e., the $C^{1,\gamma}$ continuity of $E_{\sigma D}$. Estimate (2.8) has been proven.

Next, we estimate the electric field difference. We denote $E_{i}=E_{\sigma_{i}D}$, for $i=1,2$. Note that ${E}_{1}-{E}_{2}$ is curl-free. We set

[TABLE]

Then, $u$ satisfies the equation

[TABLE]

With the same argument for proving (2.6), we obtain that

[TABLE]

Thus, we conclude from (2.7) that

[TABLE]

From the standard theory of elliptic equations

[TABLE]

which implies (2.10). ∎

3 Uniqueness and stability

With the notation of the previous section, we will show the well-posedness of the inverse problem in a certain functional space.

We prove the following theorem on the stability of the inverse problem.

Theorem 1**.**

Let $F(\sigma)=\nabla\cdot(\sigma DE_{\sigma D}\times B_{0})$. Suppose Assumption 1 is satisfied and $\sigma_{1}-\sigma_{2}\in W^{1,\infty}_{0}$. If there exist constants $K$, $L$ and $\eta$ such that

[TABLE]

[TABLE]

[TABLE]

then

[TABLE]

holds for some positive constant $c$.

Proof.

We denote $E_{i}=E_{\sigma_{i}D}$, $i=1,2$ and write the data difference as

[TABLE]

Then, we can rewrite $F(\sigma_{1})-F(\sigma_{2})$ as

[TABLE]

where

[TABLE]

where $I$ is the identity matrix.

Next we multiply $F(\sigma_{1})-F(\sigma_{2})$ by $\sigma_{1}-\sigma_{2}$ and integrate over $\Omega$. For $I_{i}$, $i=1,2,3,4$, we can estimate the integrals $\int_{\Omega}(\sigma_{1}-\sigma_{2})I_{i}$ separately. We have

[TABLE]

Here we use the equality $\nabla\cdot(E_{1}\times B_{0})=1$ which can be easily checked from the identity $\nabla\cdot(E_{1}\times B_{0})=B_{0}\cdot(\nabla\times E_{1})-E_{1}\cdot(\nabla\times B_{0})=1$. On the other hand,

[TABLE]

Here the assumption (3.1) and inequality (2.9) have been used.

Now we turn to the terms $I_{3}$ and $I_{4}$. We have

[TABLE]

In the last inequality we have used estimate (2.7) together with the assumptions (3.2) and (3.3). Finally, we have

[TABLE]

Here we have used the assumptions (3.2), (3.3) and inequality (2.9).

Let $K$ and $\eta$ be such that $KC+2\eta LC<\frac{1}{2}$. We obtain

[TABLE]

for some constant $c$, which proves the theorem. ∎

Now we are ready to introduce a functional framework for which the inverse problem is well defined. We assume that $\sigma$ is known on the boundary of $\Omega$. In what follows, we let $\sigma_{*}$, the true cross-property factor of $\Omega$, belong to a bounded convex subset of $C^{1,\beta}(\Omega)$ given by

[TABLE]

where $\sigma_{0}$ is some positive function satisfying Assumption 1 and

[TABLE]

with $c_{1},c_{2},c_{3}$ and $c_{3}$ being positive constants. In other words, we can write $\widetilde{\mathcal{S}}=\sigma_{0}+\mathcal{S}.$

It is clear that the distribution of the electric field ${E}_{\sigma D}$ depends nonlinearly on the factor $\sigma$ and $\nabla\cdot(\sigma D{E}_{\sigma D}\times{B}_{0})$ is nonlinear with respect to $\sigma$. We first examine the Fréchet differentiability of the forward operator $F$. Then, some useful properties of the Fréchet derivative at $\sigma$, $DF[\sigma]$, are presented.

To introduce the Fréchet derivative, we consider the following Neumann boundary value problem

[TABLE]

and

[TABLE]

where $h\in\mathcal{S}$ is the increment to the factor $\sigma$.

By the same arguments as those in [20], together with Theorem 1, it is natural to conclude the following result that insures the well-posedness of the inverse problem.

Theorem 2**.**

For $\sigma$ and $D$ satisfying Assumption 1 and (3.2), the operator $F$ is bounded and Fréchet differentiable at $\sigma\in\widetilde{\mathcal{S}}$. Its Fréchet derivative at $\sigma$, $DF[\sigma]$, is given by

[TABLE]

where $\varphi_{h}$ solves (3.10). Meanwhile, $DF[\sigma]^{*}$, i.e., the adjoint of $DF[\sigma]$ is defined as below,

[TABLE]

where $\Phi_{g}$ solves (3.11). Furthermore, we have the following stability result,

[TABLE]

for some constant $C$ which depends on $\lambda,c_{1},c_{2}$ and $\Omega$ and the constant $c$ is defined in Theorem 1.

Proof.

The definition of $DF[\sigma](h)$ and (3.14) follow from [20] and Theorem 1. Here we only give a brief proof of the formulation of $DF[\sigma]^{*}$.

First by multiplying (3.10) by $\Phi_{g}$ and (3.11) by $\varphi_{h}$, we get after an integration by parts,

[TABLE]

Then we are ready to compute $DF[\sigma]^{*}(g)$. We have

[TABLE]

This proves (3.13). ∎

4 The reconstruction method

4.1 Optimization scheme

It is natural to formulate the reconstruction problem for $\sigma_{*}$ as a least-square problem. To find $\sigma_{*}$ we minimize the functional

[TABLE]

over $\sigma\in\widetilde{\mathcal{S}}$.

We can now apply the gradient descent method to minimize the discrepancy functional $J$. Define the iterates

[TABLE]

where $\mu>0$ is the step size and $T[f]$ is any approximation of the Hilbert projection from $L^{2}(\Omega)$ onto $\overline{\widetilde{\mathcal{S}}}$ with $\overline{\widetilde{\mathcal{S}}}$ being the closure of $\widetilde{\mathcal{S}}$ (in the $L^{2}$-norm).

The presence of the projection $T$ is necessary because $\sigma_{n}$ might not be in $\widetilde{\mathcal{S}}$.

Using the definition of $J$ we can show that the optimal control algorithm (4.1) is nothing else than the following projected Landweber iteration [9, 11, 14] given by

[TABLE]

For completeness, we state the convergence result of Landweber scheme here without proof. We refer to [2, 14] for details.

Theorem 3**.**

The sequence defined in (4.2) converges to the true cross-property factor $\sigma_{*}$ of $\Omega$ in the following sense: there exists $\epsilon>0$ such that if $\|T[\sigma_{1}]-\sigma_{*}\|_{L^{2}(\Omega)}<\epsilon$, then

[TABLE]

4.2 A quasi-Newton method

It has been observed in [20] that the challenge of the Landweber iteration lies in the difficulty of evaluating the adjoint operator of the Fréchet derivative. To avoid taking too many derivatives, we introduce a more efficient way to reconstruct the conductivity. This is a generalization of the quasi-Newton method proposed in [20] for the anisotropic case with known conformal class. In the following, we describe this algorithm and prove its convergence in $\widetilde{\mathcal{S}}$.

Let $\sigma$ be the scalar conductivity distribution function and let $D$ be the known conformal class matrix-valued function. The forward operator is given by

[TABLE]

where $E_{\sigma}$ satisfies the system

[TABLE]

Algorithm 1**.**

Step 0. Select an initial conductivity $\sigma_{1}\in\widetilde{\mathcal{S}}$ and set $k=1$.

Step 1. Calculate the associated electric field $E_{k}$ by solving the boundary value problem

[TABLE]

Step 2. Calculate the updated conductivity by solving the stationary advection-diffusion equation with the inflow boundary condition:

[TABLE]

where $\partial\Omega^{-}=\{x\in\partial\Omega\mid DE_{k}(x)\times B_{0}\cdot\nu(x)<0\}$.

Step 3. Let $\sigma_{k+1}=T[\sigma_{k+1/2}]$, where $T$ is the Hilbert projection operator onto $\tilde{S}$. Set $k=k+1$ and go to (4.3).

4.3 Convergence analysis

In the algorithm above, one updates the electric field $E$ and then updates the cross-factor $\sigma$ later. Using the same argument as for proving the well-posedness, we could get the following convergence results.

Theorem 4**.**

Suppose that the cross-factor $\sigma_{*}\in\widetilde{\mathcal{S}}$ and $D$ satisfies Assumption 1 and (3.2). Let $\{\sigma_{k}\}$ be determined by the Algorithm 1. Then for proper constants $K,L$ and $\eta$ in (3.9) and (3.2), there exists a constant $c<1$ such that

[TABLE]

Proof.

Note that $T$ is a projection and $\widetilde{\mathcal{S}}$ is convex. Then, we have that $T$ is nonexpansive and $\|\sigma_{k+1}-\sigma_{*}\|\leq\|\sigma_{k+1/2}-\sigma_{*}\|$ follows. It is left to estimate $\|\sigma_{k+1/2}-\sigma_{*}\|$.

First we subtract $\nabla\cdot(\sigma_{*}DE_{k}\times B_{0})$ from both sides of (4.4)to get

[TABLE]

Multiplying by $\sigma_{k+1/2}-\sigma_{*}$ and integrating over $\Omega$ yields

[TABLE]

We split the terms into

[TABLE]

Integrating by parts gives $\int_{\Omega}(\sigma_{k+1/2}-\sigma_{*})\nabla\cdot((\sigma_{k+1/2}-\sigma_{*})E_{k}\times B_{0})=\frac{1}{2}\|\sigma_{k+1/2}-\sigma_{*}\|_{L^{2}(\Omega)}^{2}$.

The second term in the left hand side can be estimated as follows:

[TABLE]

Here the smallness of $D-I$ and the $C^{1}$ property of $E_{k}$ (2.8) have been used.

For the right hand side, we have

[TABLE]

and

[TABLE]

Here we have used property (2.10) and the fact that $\sigma_{k}\in\tilde{S}$.

With the above estimates, as we did in Theorem 1, let $KC+\eta(L+2)C<\frac{1}{2}$. We derive

[TABLE]

where $c$ is a constant smaller than $1$. Hence,

[TABLE]

which proves the theorem. ∎

4.4 Numerical experiments

In this section, we present some numerical experiments to validate the reconstruction method proposed in Algorithm 1 and evaluate its robustness to measurement noise. To simplify the computation, we convert this three-dimensional problem into an equivalent two-dimensional problem assuming that the domain of interest is the cube $[0,1]^{3}$ and the conductivity and the diffusion tensors are invariant along the third dimension. Moreover, we assume that the diffusion tensor $D$ is of form

[TABLE]

where $d_{ij}$’s are constant plus some perturbations as shown in Figure 1. The non-zero part of the perturbation functions are used to characterize the anisotropy.

We use a uniform finite element triangular mesh over the two-dimensional unit square. The number of cells is $256$ in each direction. The total number of triangles and vertices are $2^{17}$ and $257^{2}$, respectively. Both the elliptic equation with a Neumann boundary condition and the stationary advection-diffusion equation are solved using the finite element method of first order implemented with FEniCS [17]. The internal data $F(\sigma)$ used for the reconstruction are synthetic data that are generated using the same solver. These data are commonly used to refer to the “noise-free” data, although they may contain some numerical errors.

For all examples, we use the same initial guess, constant function $0.2$, and the same true cross-property factor (Figure 2a) given by

[TABLE]

where $r(x_{1},x_{2})=\sqrt{(x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}}$ and $s=\frac{0.46-r}{0.12}$. The internal data generated with the diffusion tensor as in (4.10) is shown in Figure 2b. We also produce the data in the isotropic case (Figure 2c). The effect of the anisotropy can be observed clearly. The error-decay of the reconstruction with the noise-free data is shown in Figure 4a. The final error is smaller than $2\times 10^{-3}$. We only display the last iterate here.

This inverse problem bears a Lipschitz type stability and we avoid lowering the regularity of the cross-property factor using Algorithm 1. Therefore, the robustness of the reconstruction scheme to noisy data is expected. We perform the numerical tests with noisy data by perturbing the internal functional $g$ in the following way:

[TABLE]

where $w$ is a function taking values uniformly distributed in $[-1,1]$ and $\delta$ is the noise level.

Figure 3 shows the noisy data with noise level $\delta=24\%$ and the reconstructed cross-property factor. We do not use further regularization techniques since the regularization method may depend on the type of the noise in practical cases. But the projection onto the feasible space acts as a regularization scheme.

5 Concluding remarks

In this paper, we have considered the reconstruction of an anisotropic conductivity from MAT-MI data which is conformal to a known diffusion tensor measured from DTI. The data is the internal functional $\nabla\cdot(\sigma E_{\sigma}\times B_{0})$ throughout the domain. We have analyzed the linearization of the problem and the stability of the inversion. A local Lipschitz type stability estimate has been established for a certain class of anisotropic conductivities. A quasi-Newton type reconstruction method with projection has been introduced and its convergence has been proved. Numerical experiments demonstrate the effectiveness of the proposed approach and its robustness to noise.

In light of the numerical experiments, we have the following observations.

The effect of the electrical anisotropy is remarkably significant and can not be neglected in the reconstruction of electrical conductivity in MAT-MI. 2. 2.

There is still a room for improvement of the admissible class of conductivities. The convergence of the proposed algorithm has been observed for more general cases. 3. 3.

For the inversion with noisy data, oscillation in the reconstructed conductivity is observed. Regularization methods prompting sparsity, such as total variation regularization may be employed for a more stable reconstruction.

The present work leads to many important unanswered questions. First, is it possible to determine an anisotropic conductivity $\gamma(x)$ from MAT-MI data alone? Recall that each static field $B_{0}$ leads to a new measurement. It would be important to know how many static fields and their associated measurements are needed to determine $\gamma(x)$ uniquely. Second, how can the present approach be modified to address this very question? It is clear that the current analysis relies on vector identities which do not easily generalize to handle anisotropic conductivity. Finally, if we can answer these questions, then how do we devise stable and accurate computational methods to image anisotropic conductivity? We therefore must accept the conclusion that the present work is the first step towards a fully MAT-MI method of imaging anisotropic conductivity.