Qun CHENJijun LIU

1 Introduction

Magnetic resonance electrical impedance tomography(MREIT,for short)is a new electrical conductivity imaging technique to visualize the cross-sectional images of a conductivity distributionσof biologic tissues.In contrast to the traditional electrical impedance tomography(EIT,for short)technique(see[1,7,18]),this new technique applies essentially the internal electrical current distribution to recover the conductivity,which weakens the ill-posedness of EIT problem and provides a higher resolution of conductivity image.

In MREIT,we place a subject inside a magnetic resonance imaging(MRI,for short)scanner and inject a currentIbetween two electrodes attached on its boundary.Then there exists the internal currentJ=(Jx,Jy,Jz)inside the subject,generating a magnetic flux densityB=(Bx,By,Bz).Herez-axis is the direction of the main magnetic field of the scanner.TheBzdata can be measured by using the MRI scanner,from which we try to reconstruct the bio-tissue conductivity,see Figure 1 for the configuration of this system.

Figure 1 MREIT system at impedance imaging research center and harmonic Bz algorithm mathlab toolkit

Recently,some reconstruction schemes usingBzdata as inversion input have been proposed,such as harmonicBzmethod,gradientBzmethod and variational gradientBzmethod(see[2,12,17]).It has been proven that the measurementsBz,jcorresponding to two incoherent injection currentsIjwithj=1,2 can determine the conductivity distributionσuniquely in 2-dimensional case(see[3,14])under some a priori assumptions.

The harmonicBzalgorithm was the first constructive MREIT imaging method based onBzdata(see[17]).From the Ampere law

and

we have

whereµ0is the magnetic permeability of the free space,u[σ]as a nonlinear function ofσis the induced electrical potential satisfying a nonstandard PDE problem specified in Section 2.

Taking thez-component of(1.3),it follows that

Corresponding to two incoherent injected currentsIj,j=1,2 through two pairs of surface electrodesand,it follows from(1.4)that

where

anduj[σ],Bz,jare the potential and the magnetic flux density,respectively,corresponding toIjwithj=1,2.

The harmonicBzalgorithm is an explicit iteration scheme to approximateσat each 2-dimensional slice Ωz0=Ω∩{z=z0}⊂R2based on the relation(1.5).Since the harmonicBzalgorithm was proposed,it has been improved rapidly in various numerical simulations and phantom experiments(see[8–11,15–16]).In[4],the authors proved that,for a relatively small contrast of the target conductivity,the iterative harmonicBzalgorithm based on(1.5)with a good initial guess is stable and convergent in the continuous norm.In[5],the author improved the convergence result based on the following equivalent equality of(1.5):

and derived a posteriori error estimate ofwhereσ∗is the true conductivity.However,these two convergence results are considered only for exact magnetic flux fieldBz.

In practical situations,we can only acquire the noisy dataofBzusing MRI equipment.However,the harmonicBzalgorithm applies in fact the Laplacian ofBzas inversion input,the noise contained inBzwill be amplified by such an operation and therefore has essential influence on the approximation accuracy of the iteration solution.Such an influence depends not only on the error level of noisy input data,but also on the regularizing strategy computing the Laplacian from the noisy data.So it is necessary to give an error estimate on the iterative harmonicBzalgorithm for the noisy input datacorresponding to some regularization scheme for the practical applications of harmonicBzalgorithm,which is the purpose of this paper.

This paper is organized as follows.In Section 2 we state the mathematical formulation of the harmonicBzalgorithm.Then for a relatively small contrast of the target conductivity,the iteration error of this algorithm is established in Section 3 for noisy input data,under the assumption that a stable numerical differentiation process has been applied.In Section 4,we propose a numerical regularizing scheme for the computation of Laplacian from the noisy measurement datawith error estimate,which provides the basis on computing the error ofiterative solution of conductivity.

2 Mathematical Model

Let Ω⊂R3be an electrically conducting subject with its smooth connected boundary∂Ω.In MREIT,we inject a currentIthrough a pair of surface electrodesε±,then it produces an internal current densityJ=(Jx,Jy,Jz)inside the subject Ω satisfying the following problem:

wherenis the outward unit normal vector on∂Ω and dsis the surface area element.

SinceJ=−σ∇u[σ],(2.1)can be converted to

This exact model(2.2)can be solved in terms of the following standard problem(see[4]):

More precisely,ifu[σ]and[σ]are the solution of problems(2.2)–(2.3),respectively,then

whereCis a constant decided by the electric potential specified at one point.

We now consider the magnetic field produced by the injection currentI.From the Biot-Savart law,it follows that

which generates the following relation betweenBzandσfrom(1.2):

Recently,a new iteration scheme based on the nonlinear integral equation(2.5)was proposed in[6],which applies theBzdata as the inversion data directly in the algorithm,without the computation of Laplacian on the magnetic flux.

The harmonicBzalgorithm is an iterative scheme at each 2-dimensional slice Ωz0=Ω∩{z=z0}based on the identity(1.5).To give the complete iteration scheme,we introduce the fundamental solution of 2-dimensional Laplace operatorsatisfying∆r′Φ(r,r′)=δ(r−r′)forr∈R2,then at each 2-dimensional slice,it holds that

where∇=(∂x′,∂y′),

It has been noticed thatA[σ]−1(x,y,z0)may be large near∂Ωz0due to the fact that two induced currentsσ∇u1[σ],σ∇u2[σ]are probably almost parallel for some configuration.This phenomena may lead to the unconvergence of the iteration scheme.To avoid this difficulty,we assume as usual that the unknown true conductivity is constant infor some interior domainThen it is easy to see from(1.4)that

We denote byσ∗the true unknown conductivity and assume that its value on∂Ωz0,still denoted asσ∗,is known.LetBz,j,j=1,2 be the exact magnetic flux density corresponding toσ∗for two inject currents.For given initial guessσ0(x,y,z0)in Ωz0with exact values inthe harmonicBziteration algorithm constructs an approximation sequence{σn(x,y,z0):n=0,1,2,···}from

for(x,y)based on the relations(1.5)and(2.6),whereisH(σ∗)with∂Ωz0replaced byFor(x,y)it is obvious thatσn(x,y,z0)≡σ∗(x,y,z0)from the first equation in(2.7)since

To give the error estimate for the iteration solution with noisy input data in the next section,we need two regularity results for direct problems.

Lemma 2.1Denote by E the regular open subsurface of the boundary∂ΩofΩ⊂R2.Then for the boundary value problem

with σ∈L∞(Ω)satisfyingthe following estimateshold:

Iff∈(L2(Ω))2and σ∈C(Ω),then u∈H1(Ω)and

iff∈(H1(Ω))2and σ∈C1(Ω),then uand

iff∈(C0,α(Ω))2with α∈(0,1)and σ∈C1(Ω),then uand

iff∈(Lp(Ω))2with p>1and σ∈C(Ω),then uand

whereare regular domains,and Ci(Ω)have the following forms:

The functions Fi(i=1,2,3,4)are known bounded continuous functions with respect to the arguments.

This result can be found in[4].

Lemma 2.2Letbe the solution of the following problem:

Then there exists a constant C(σ)such that

where C(σ)=(CsC3(σ)+1)C1(σ)C2(σ),⊂⊂Ω.

ProofIt follows from(2.8)–(2.9)that

where⊂⊂⊂⊂Ω.By the Sobolev imbedding theorem,we can obtain

for everyα∈(0,1).

Finally,combining these two estimates with(2.10),we have

which leads to(2.14).

3 Error Estimate for Noisy Input Data

We consider the error estimate of harmonicBziteration algorithm in axially symmetric cylindrical sections.Let Ω be a cylinder along thezdirection with infinite length and the electrode pair be parallel to thezdirection.Moreover we assume that the conductivityσ∗inΩ does not change alongzdirection.Then the conductivity is actually reconstructed in the 2-dimensional domain.To unify the notations,we still use Ω in the sequel to represent the 2-dimensional domain Ωz0.

In this section,we consider the error estimate of harmonicBzalgorithm for noisy input dataIn this noise input data situation,for given initial guessσ0(x,y)in Ω with exact values inΩ,where⊂⊂Ω,the iterative sequence{σn,δ:n=1,2,···}is generated from

and then forn=1,2,···,

where the notationis the approximation to∇2Bz,for which we will propose a regularizing scheme in Section 4.

We firstly give two known results related to the convergence for the exact magnetic field,which will be applied in our error estimate.

Lemma 3.1For⊂⊂Ω,if σ lies in the set

whereσ0,λ,ϵ0are positive constants,then there exists a constant d∗−depending only on λ,ϵ0,Ω,dist(∂Ω,)and,such that

Lemma 3.2Assume that the target conductivitymeets the followingconditions:

(H1)for known constants

(H2)there exists⊂⊂Ωsuch that σ∗is a known constant inΩ;

(H3)|detwhereis a known constant.

Under these hypotheses,there exist constantssmall enough and θ=such that if we take the initial guess σ0as the constantthenthe sequence σngiven by the harmonic Bziteration algorithm using exact input data holds forthat

where K:=diam(Ω)+1.

These two results can be found in[4–5],respectively.

From Lemma 3.2,for true conductivityσ∗lying in the set

the iterative sequence{σn:n=1,2,···}using exact input data lies in

for any

inand

Noticing thatσn≡σ∗in Ω,we get(3.5).

For practical measurement data with noise,the input data for the iteration scheme is in fact the Laplacian operationfrom(3.1)–(3.2).When presenting our error estimate,we must firstly analyze the errorρ(δ)of computing∇2Bzfrom the noisy measurement datawhich depends on the regularizing scheme.Since we generally measure the error of magnetic field itselfinL2-norm,while we need the error estimate of Laplacian inC-norm in our iteration,we give the following approximation for our computation on Laplacian:

(H4)For the noisy datasatisfying

a stable differentiation scheme is used to computesuch that

wherej=1,2.

(H5)is understood such thatin Ω,which is a natural condition if the conductivity is assumed to be a known constant in the domain Ωimplying

An implementa√ble regularizing scheme to compute∆Bzapproximately fromto reach(3.6)withas well as the regularity requirement on the target conductivity will be given in Section 4.

Now we can state the main result of our work as follows.

Theorem 3.1Assume that the target conductivitymeets the three hy-potheses in Lemma3.2and(H4)–(H5).Then there exist constantssmallenough andsuch that if we take the initial guess σ0as the constantthe sequence{σn,δ}given by(3.1)–(3.2)with noisy input data holds for≤ϵ and δ small enough that

ProofLet us takeDenote byandthe solutions of the direct problem

withσ=σ∗andσ=σn,respectively.It follows from Lemma 2.2 that

and

However,Lemma 3.2 says that{σn:n=1,2,···}⊂S2.So it follows from the expressions ofC(σ∗)andC(σn)in Lemma 2.2 that the constants are of a uniform upper bound:

with

Step 1 Estimate

Firstly,expand the initial guessσ0atσ∗asσ0=σ∗+e0.Sinceandσ0=σ∗in Ω,it follows that

Hence,(diam(Ω)+1)ϵ=:Kϵ.

We expandatas

Noticing thatσ0=σ∗in Ω,meets

Sinceande0=0 in Ω,it follows from(3.8)that the right-hand side of the first equation in(3.11)satis fies

Therefore it follows from Lemma 2.1 and the Sobolev imbedding theorem that

According to the above estimate and(2.10),we have for⊂⊂⊂⊂Ω that

Using the same arguments as those in deriving(3.8),we can get

Therefore we have

Denote by

a known function due to Lemma 2.1.Forwe introduce the constant

which is well defined.Noticing thatwe havefor 0<ϵ

Now it follows from(3.12)–(3.15)that

Sinceit follows that

On the other hand,it follows from(2.7)and(3.1)that

which can be written as

due to the definition of the matrixA[σ0]and(3.10).

However,it is obvious from(3.17)that

A direct computation leads tofrom which we deduce

due to(3.8)and(H3).

Now we takesmall enough such that

which implies from(3.20)that

Now it follows from(3.19),(3.21)–(3.22)that

where the Sobolev embedding theorem[0,1)on∇2ejbased on(3.6)is applied.This last estimate generates

Introduce a new constant

then the estimate(3.23)becomes

On the other hand,it follows from(3.1)and(H5)that∇σ1,δ=0 in Ω,and thereforeσ1,δ=σ∗in Ω.

Takeδsmall enough such thatMρ(δ)Then it follows from the above estimate and Lemma 3.2 that

for anySoσ1,δlies in the set

Now we can apply the induction argument to prove the theorem.That is,assume that the following properties

hold fork=n,which specially yields that

noticingσn∈S2.Then we need to prove that these properties are also true fork=n+1.Step 2 Expandσn,δatσn.

Let(j=1,2)be the solutions of the problem(3.7)withσ=σn,δ.It follows from Lemma 2.2 that

Sinceσn,δ∈S3,similarly to the estimate ofC(σn)in(3.9),there exists a constantdepending only on the upper and lower boundsofσ∗,Kand domain,still denoted byfor the simplicity of notation,such that

We expandand

Sincesatisfies the following problem:

Similarly to the derivation of(3.16),we have

where the definition ofis similar todue to S2⊂S3.

On the other hand,it follows fromand Lemma 3.1 that

whereis a constant depending only ondistand

Step 3 Estimate

It follows from(3.2)and(2.7)that

Hence we have

Firstly,we estimate II.From the definition ofA[σ]in(1.6),we have

On the other hand,the Sobolev imbedding theorem and(H4)yield

for the simplicity of notation.So it follows from(3.28),(3.30)and(3.34)that

Then we estimate I in(3.32).Again using the definition ofA[σ],we get that

A direct computation leads to

which yields

from(3.9),(3.28)and(3.30).Therefore it follows from(3.28)and(3.37)that

Again from(3.28),we have

By(3.36),(3.38)–(3.39),we obtain that

So it follows from(3.29)and(3.40)that

On the other hand,(1.5)yields

from(3.8)and the condition

Finally combining(3.41)and(3.42)together yields

Inserting(3.35)and(3.43)into(3.32),we get

This estimate together withleads to

Now we takesmall enough such that

and then it follows from(3.44)that

Inserting(3.26)fork=ninto(3.45),we can get

Moreover,we conclude that∇σn+1,δ=0 in Ωfrom(3.2)and(H5),and thereforein Ω

Then it follows from Lemma 3.2,(3.46)and the triangle inequality that

The proofis complete.

The important conclusion derived from Theorem 3.3 is that,different from using the exact input magnetic field,the iteration solutionσn,δof harmonicBzalgorithm using noisy magnetic fieldcan only approximate the exact conductivityσ∗up to a finite accuracy.More precisely,it follows from Theorem 3.3 that

for any fixed error levelδ>0.This is reasonable from the general iteration scheme based on nonlinear integral equation of the second kind that the accuracy of the kernel determines the accuracy of solution,which can not be improved by increasing the iteration times.

To havewe must choosen→∞andρ(δ)→0 simultaneously.The total errorconstitutes of two parts:Mρ(δ)andKθnε.The former depends on the noise levelδand the strategy computing the Laplacian such thatρ(δ)→0;while the later describes the iteration error which can be improved by increasingn.Notice that we need to distinguish two different input errors:Input data errorδfor our MREIT problem and input errorρ(δ)for the harmonicBzalgorithm of MREIT problem.The efficient realization of harmonicBzalgorithm depends on decreasing both the iteration errorKθnand the input errorρ(δ)for the algorithm.In the next section,we will analyze the input data errorρ(δ).

4 Stable Computation for Laplacian of

We propose a regularizing scheme for computing 2-dimensional∇2Bzapproximately from the noisy input datain.Noticing that we can take⊂⊂Ω such thatlocates in the domain,whereσ∗is a known constant,we can assume thatBz(x)is exactly specified nearwhich means

Since we need theestimate for the Laplacian computation in Theorem 3.3,we assume that the exact magnetic field is approximated by its noisy measurement data in the sense

due to(4.1).Notice that the spacecan be characterized in terms of boundary conditions for(see[13,Theorem 7.41]):

Denote bythe fundamental solution of−∆operator,i.e.,

Then for exactBz(x),its Laplacian∆Bz(x):=f(x)inmeets

Moreover,by the Newtonian potential method for Poisson’s equation,we know thatBz(x)∈noticingf(x)≡0 nearfrom(1.4)and the choice of

Now let us define a linear bounded mapby

Then(4.4)can be written as

from(4.3)andf(x)≡0 near

For noisy input datasatisfying(4.1)–(4.2),we defineas the solution to the following integral equation of the second kind:

for some suitable regularizing parameterα=α(δ)>0,whereis the adjoint operator ofObviously,(4.7)is the Tikhonov regularizing equation to the following integral equation of the first kind:

The unique solvability of(4.7)follows from the standard Tikhonov regularizing theory.Notice thatis neither self-adjoint nor nonnegative.

The first result in this section is the error estimate onfδ−f.

Theorem 4.1Assume thatIf we choose the regularizing α=δ,then wehave for the noisy input datasatisfying(4.2)that

ProofIt follows from(4.6)and(4.7)that

noticing(4.1).Sincefor somedue tothe above equation becomes

Then we have

from the standard estimate on the Tikhonov regularizing operator.The proofis complete by takingα=δ.

Remark 4.1This result is a classical a priori choice strategy for Tikhonov regularization.The source conditionimplies the requirement thatBzshould be very smooth,not necessarily inL2(Ω).The other a posterior strategies such as discrepancy principle can also be considered under the framework of Tikhonov regularizing.However,sinceis not nonnegative,the scheme applying the Lavrentive regularizing for computing the Laplacian(see[19])does not work.

To avoid the explicit expression of the adjoint operatorwe consider(4.7)directly.Using the property of definition ofinner product,it follows for allχ∈andthat

whereK∗is the adjoint operator ofKin the sense mappingto itself.Therefore the regularizing equation(4.7)inhas the following equivalent weak form:

However,it is easy to verify thatis self-adjoint.In fact,for anywe have

Therefore(4.12)is equivalent to

with

Now we consider how to solve(4.13).We define the inner product in Hilbert spaceby

which yields the equivalent norm toIntegrating by parts yields

Therefore(4.13)is equivalent to the following integral-differential system:

noticing thatis dense in

We can solve this well-posed system to get the approximation of∇2Bzfrom noisy dataOnce we generate∇2Bzfrom this system,which is a good approximation to∇2Bz,the harmonicBzalgorithm can be implemented efficiently.

[1]Cheney,M.,Isaacson,D.and Newell,J.C.,Electrical impedance tomography,SIAM Rev.,41,1999,85–101.

[2]Kwon,O.,Park,C.J.,Park,E.J.,et al.,Electrical conductivity imaging using a variational method inBz-based MREIT,Inverse Problems,21,2005,969–980.

[3]Kwon,O.,Pyo,H.C.,Seo,J.K.and Woo,E.J.,Mathematical framework for Bz-Based MREIT model in electrical impedance imaging,Comput.and Math.with Appl.,51(5),2006,817–828.

[4]Liu,J.J.,Seo,J.K.,Sini,M.and Woo,E.J.,On the convergence of the harmonicBzalgorithm in magnetic resonance electrical impedance tomography,SIAM J.Appl.Math.,67,2007,1259–1282.

[5]Liu,J.J.,Seo,J.K.and Woo,E.J.,A posteriori error estimate and convergence analysis for conductivity image reconstruction in MREIT,SIAM J.Appl.Math.,70,2011,2883–2903.

[6]Liu,J.J.and Xu,H.L.,Reconstruction of biologic tissue conductivity from noisy magnetic field by integral equation method,Appl.Math.and Comput.,218,2011,2647–2660.

[7]Metherall,P.,Barber,D.C.,Smallwood,R.H.and Brown,B.H.,Three-dimensional electrical impedance tomography,Nature,380,1996,509–512.

[8]OH,S.H.,Lee,B.I.,Woo,E.J.,et al.,Conductivity and current density image reconstruction using harmonicBzalgorithm in magnetic resonance electrical impedance tomography,Phys.Med.Biol.,48,2003,3101–3116.

[9]OH,S.H.,Lee,B.I.,Woo,E.J.,et al.,Electrical conductivity images of biological tissue phantoms in MREIT,Physiol.Meas.,26,2005,279–288.

[10]Park,C.,Kwon,O.,Woo,E.J.and Seo,J.K.,Electrical conductivity imaging using gradientBzdecomposition algorithm in magnetic resonance electrical impedance tomography(MREIT),IEEE Trans.Med.Imag.,23,2004,388–394.

[11]Park,C.,Park,E.J.,Woo,E.J.,et al.,Static conductivity imaging using variational gradientBzalgorithm in magnetic resonance electrical impedance tomography,Physiol.Meas.,25,2004,257–269.

[12]Park,C.,Woo,E.J.,Kwon,O.and Seo,J.K.,Static conductivity imaging using varia-tional gradientBzalgorithm in magnetic resonance electrical impedance tomography(MREIT),Physiol.Meas.,25,2004,257–269.

[13]Renardy,M.and Rogers,R.C.,An Introduction to Partial Differential Equations,TAM,Vol.13,Springer-Verlag,New York,2004.

[14]Seo,J.K.,A uniqueness result on inverse conductivity problems with two measurements,J.Fourier Anal.Appl.,2,1996,515–524.

[15]Seo,J.K.,Kwon,O.and Woo,E.J.,Magnetic resonance electrical impedance tomography(MREIT):conductivity and current density imaging,J.Phys.:Conf.Ser.,12,2005,140–155.

[16]Seo,J.K.,Pyo,H.C.,Park,C.,et al.,Image reconstruction of anisotropic conductivity tensor distribution in MREIT:Computer simulation study,Phys.Med.Biol.,49,2004,4371–4382.

[17]Seo,J.K.,Yoon,J.R.and Woo,E.J.,Reconstruction of conductivity and current density images using only one component of magnetic field measurements,IEEE Trans.Biomed.Eng.,50,2003,1121–1124.

[18]Webster,J.G.,Electrical Impedance Tomography,Adam Hilger,Bristol,1990.

[19]Xu,H.L.and Liu,J.J.,On the Laplacian operation with applications in magnetic resonance electrical impedance imaging,Inverse Problems in Science and Engineering,21,2013,251–268.