# Current Density Imaging Using Directly Measured HarmonicBz Data in MREIT

Authors: Chunjae Park; Oh In Kwon
Journal: Computational and Mathematical Methods in Medicine
 (2013)

---

## Abstract

Magnetic resonance electrical impedance tomography (MREIT) measures magnetic flux density signals through the use of a magnetic resonance imaging (MRI) in order to visualize the internal conductivity and/or current density. Understanding the reconstruction procedure for the internal current density, we directly measure the second derivative ofBz data from the measured k-space data, from which we can avoid a tedious phase unwrapping to obtain the phase signal of Bz. We determine optimal weighting factors to combine the derivatives of magnetic flux density data, ∇2Bz, measured using the multi-echo train. The proposed method reconstructs the internal current density using the relationships between the induced internal current and the measured ∇2Bz data. Results from a phantom experiment demonstrate that the proposed method reduces the scanning time and provides the internal current density, while suppressing the background field inhomogeneity. To implement the real experiment, we use a phantom with a saline solution including a balloon, which excludes other artifacts by any concentration gradient in the phantom.


---

## Body

## 1. Introduction

Magnetic resonance electrical impedance tomography (MREIT) visualizes a cross-sectional conductivity and/or current density inside the human body. The MREIT technique injects currents through attached electrodes in order to probe the imaging subject and then measures the induced magnetic flux density, inside the subject using an MRI scanner. The MRI scanner only measures thez-component of the induced magnetic flux density B=(Bx,By,Bz); therefore the MREIT techniques have focused on the reconstruction of the internal conductivity and/or current density by using the measurable Bz data instead of subject rotation [1–18].The MREIT techniques used to reconstruct the conductivity and/or the current density have been widely developed and have reached the stage of imaging experiments for live animals and the human body [18, 19]. Due to the poor SNR of measured Bz data in current MREIT experiments, it is critical to reduce the scan time in MREIT, while maintaining the spatial-resolution and sufficient contrast, for practical in vivo implementations of MREIT.In order to increase the quality of measuredBz data, a measurement technique called the injected current nonlinear encoding (ICNE) method was developed, which extends the duration of the injection current until the end of the read-out gradient in order to maximize the signal intensity of the magnetic flux density [20]. Motivated by the ICNE pulse sequence method, an ICNE-multiecho technique was developed and optimized by finding an optimal weighting factor for the multiple measured Bz data [21].The MREIT technique typically uses an interleaved acquisition, which scans each phase encoding consecutively by injecting two currents possessing positive and negative polarities with the same scan duration and amplitude to double theBz signal and cancel out the background field inhomogeneity. In order to reduce the scan time, for the measurement of Bz, [22] reconstructed the phase signal Bz by filling a partial k-space region using the interleaved measurement property.Functional MRI (fMRI) has been applied to a wide range of neuroscience researches by visualizing neural activities inside the brain in a fast and direct way [23–25]. A fast MREIT imaging technique has been proposed as a promising imaging technique for the continuous monitoring of internal electrical property inside the subject [26]. In this paper, we propose a method to monitor spatial and temporal internal current density changes in the subject by using a fast gradient multi-echo pulse sequence to maximize the measured Bz signal in a short scanning time. Moreover, we derive a direct method to measure ∇2Bz instead of Bz data from the measured k-space data. The proposed method can also avoid a tedious unwrapping procedure, which may introduce phase artifact in the recovered phase signal.For the recovery of the internal current density, we investigate the reconstruction procedure for the internal current density from the measured∇2Bz data. In the paper [27], a projected current  JP  was provided by the decomposition J=JP+JN, where J is the internal current density influenced by the injected current and JP is a determined component of J from the measured Bz data. The projected current JP is identical to the true current J when the z-component Jz of J=(Jx,Jy,Jz) is the same as Jz0 where Jz0 is the z-component of the background current J0.The projected currentJP can be determined in a concrete form which consists of the background current J0 and the solution of a two-dimensional harmonic equation with the Dirichlet condition that matches the external injection current on the surface of the subject. To recover the internal current density JP with the generated Bz caused by the injected current, we only use the second derivatives of Bz and ∇2Bz, which are required to solve the two-dimensional harmonic equation for JP.To reduce the noise artifact, we apply the ICNE-multi-echo train based on the fast gradient echo and solve an optimal weighting factor of∇2Bzℓ,  ℓ=1,…,NE, where NE denotes the number of echoes at each RF pulse.In order to verify how the proposed method works, we designed a phantom with a saline solution and fixed a balloon inside the phantom, inflating the balloon by injecting the same saline solution. The phantom was designed to provide a homogeneous magnitude image, but the recovered current density distinguishes inside the balloon. For the experiment, the total scan time was 12.36 seconds to obtain the completek-space data using the interleaved acquisition in order to implement the proposed method with a 128×128 spacial matrix size. The phantom experiment demonstrates that the proposed method reduces the scanning time and recovers internal current density, while suppressing the measured noise artifact.


## 2. Methods

We inject the currentI through the attached electrodes on a three-dimensional cylindrical conducting object Ω with its conductivity distribution σ. The injection currents I produce the voltage distribution u satisfying the following elliptic partial differential equation:
(1)∇·(σ∇u)=0inΩ,-σ∇u·ν=gon∂Ω,∫∂Ω‍uds=0,
where ν is the outward unit normal vector and g is the normal component of the current density on ∂Ω. Clearly, ∫∂Ω‍gds=0 due to ∇·(σ∇u)=0 in Ω. The current density J=-σ∇u and the magnetic flux density B=(Bx,By,Bz) in Ω satisfy the Ampère law J=∇×B/μ0 and Biot-Savart law, where μ0=4π10-7 Tm/A is the magnetic permeability of the free space. We let Bz(x,y)=Bz(x,y,z0) where z0 is the center of a selected imaging slice.


### 2.1. Measurement ofBz Using Interleaved ICNE Acquisition

For the interleaved ICNE acquisition method, we inject the positive and negative currents,I+ and I-, through the attached electrodes by scanning each phase encoding consecutively. For a standard spin echo pulse sequence without current injection, the k-space MR signal can be expressed as
(2)S(n,m)=∫ℝ2‍ρ(x,y)eiδ(x,y)e-i2π(Δkxnx+Δkymy)dxdy,
where ρ(x,y) is the real transverse magnetization, δ denotes the phase artifact of background field inhomogeneity, and Δkx and Δky are the reciprocals of fields of view for the x direction and y direction, respectively. During the data acquisition, we set
(3)Δkx=γ2πGxΔt,Δky=γ2πΔGyTpe,
and sample the data in (2) finitely for
(4)n=-Nx2,…,Nx2-1,m=-Ny2,…,Ny2-1.
The names of the above parameters are  γ=26.75×107rad/T·s: the gyromagnetic ratio of the proton,Gx: the magnetic reading gradient strength,Δt: the sampling interval,ΔGy: the phase encoding step,Tpe: the phase encoding time interval,Nx: the number of sampling points,Ny: the number of phase encoding lines.For the conventional MREIT case, we inject the current for the duration ofTc0 from the end of the 90∘ RF pulse to the beginning of the reading gradient. In this case, the induced magnetic flux density Bz due to the injection current provides the additional dephasing of spins and consequently the extra phase is accumulated during Tc0. The corresponding k-space data for the injection currents I± can be represented as
(5)SI±(n,m)=∫ℝ2‍ρ(x,y)eiδ(x,y)e±iγTc0Bz(x,y)×e-i2π(Δkxnx+Δkymy)dxdy.
Using the notations
(6)ρ±(x,y)=ρ(x,y)eiδ(x,y)e±iγTc0Bz(x,y),
we can compute the magnetic flux density Bz as
(7)Bz(x,y)=12γTc0tan-1(α(x,y)β(x,y)),
where α and β are the imaginary and real parts of ρ+/ρ-, respectively.In the conventional MREIT case, the noise standard deviation of the measuredBz, sdBz, is given as a known quantity, which is inversely proportional to the current injection time Tc and the SNR of the MR magnitude image Υ as follows [28, 29]:
(8)sdBz(r)=12γTcΥ(r).
Since the ICNE MR pulse sequence injects the current until the end of a read gradient, the total current injection time of the ICNE case is TcT:=Tc0+Ts and the k-space data is represented as
(9)SC±(n,m)=∫ℝ2‍ρ(x,y)eiδ(x,y)e±iγ(Tc0+Ts/2)Bz(x,y)×e-i2π(Δkxnx+Δkymy)dxdy,
where Ts=ΔtNx is the data acquisition time. In the usual spin echo, the ICNE current injection method demonstrates better SNR in the measured magnetic flux density data than the conventional current injection method. The optimal data acquisition time Ts* has been calculated for the usual spin echo as
(10)Ts*=233TE-3τrf,
which optimally reduces the noise in the Bz data, where τrf is the time of RF pulse [30].In the ICNE MR pulse sequence case, the noise standard deviation of the measuredBz,  sdBzICNE, is given as follows [30]:
(11)sdBzICNE(r)=12γ(Tc0+Ts*/2)Υ(r).
The prolonged data acquisition time, however, may suffer from undesirable side artifacts such as blurring, chemical shift, and motion artifacts along the phase encoding direction. To reduce the undesirable side artifacts, we divide the prolonged data acquisition time into several short ones in the ICNE-multi-echo MR pulse sequence.


### 2.2. Measurement of∇~2Bz Using ICNE-Multiecho Train

Using the ICNE-multi-echo MR pulse sequence, the measuredk-space data can be represented as
(12)Sℓ±(n,m)=∫ℝ2ρℓ(x,y)eiδℓ(x,y)e±iγTcℓBz(x,y)×e-i2π(Δkxnx+Δkymy)dxdy,ℓ=1,…,NE,
where NE is the echo number, Tcℓ is the ℓth time width of the injected current, and ρℓ and δℓ denote the ℓth transverse magnetization and phase artifact, respectively.Figure1 presents a schematic diagram for the ICNE-multi-echo MR pulse sequence based on a gradient echo pulse sequence. By taking the inverse fast Fourier transform, the ICNE-multi-echo sequence generates multiple complex images with different magnitude amplitudes depending on T2* decay and different widths of current injection time:
(13)ψℓ±(x,y):=ρℓ(x,y)eiδℓ(x,y)e±iγTcℓBz(x,y),ℓ=1,…,NE.
Using the relation (13), we derive a formula for ∇~2Bzℓ as
(14)∇~2Bzℓ(x,y)=1iγTcℓ∇~·(|ψℓ+(x,y)|ψℓ+(x,y)∇~ψℓ+(x,y)|ψℓ+(x,y)|-|ψℓ-(x,y)|ψℓ-(x,y)∇~ψℓ-(x,y)|ψℓ-(x,y)|),
where ∇~Bzℓ=(∂Bzℓ/∂x,∂Bzℓ/∂y) denotes the two-dimensional gradient of Bzℓ. The induced ∇~2Bzℓ in (14) removes the low-frequency phase artifact δℓ by subtracting ∇~2δℓ in (14).Figure 1
Diagram of the ICNE-multi-echo MR pulse sequence based on a gradient echo.The calculated vector(|ψℓ±|/ψℓ±)∇~(ψℓ±/|ψℓ±|) corresponding to iγTcℓ∇~(δℓ±Bzℓ) includes unavoidable measured noise. When we consider the decomposed form of (|ψℓ|/ψℓ)∇~(ψℓ/|ψℓ|)=∇~f+∇~×Ψ, where the curl term ∇~×Ψ is a part of unavoidable measured noise, the divergence procedure for ∇~2Bzℓ in (14) cancels ∇~×Ψ, and therefore the measured ∇~2Bzℓ includes a denoising procedure by suppressing a part of the measured noise.


### 2.3. Optimal Combination of Measured∇~2Bzℓ,  ℓ=1,…,NE

The measured∇~2Bzℓ,  ℓ=1,…,NE includes different amounts of unavoidable noise since the intensity of transverse magnetization and the width of injected current are different at each echo.The noise standard deviation of∇~2Bzℓ in (14) is given as
(15)sd∇~2Bzℓ(x,y)=CγTcℓΥℓ(x,y),
where the constant C only relates to the numerical differentiations for ∇~2Bzℓ, and Υℓ denotes the SNR of the ℓth MR magnitude image.Since the noise levels of the measured∇~2Bzℓ,  ℓ=1,…,NE in (14) are given as known quantities, we can utilize the known information sd∇~2Bzℓ to determine an optimized ∇~2Bz which combines the multiple ∇~2Bzℓ:
(16)∇~2Bz(x,y)=∑ℓ=1NE‍ωℓ(x,y)∇~2Bzℓ(x,y).
The problem of determining the weighting factors ωℓ for ∇~2Bzℓ can be formulated as
(17)minωℓ(x,y),ℓ=1,…,NE∑ℓ=1NE‍(ωℓ(x,y))2Var∇~2Bzℓ(x,y)subjectto∑ℓ=1NE‍ωℓ(x,y)=1,ωℓ(x,y)>0,
where Var∇~2Bzℓ denotes the noise variance of ∇~2Bzℓ,  ℓ=1,…,NE in (14).Following similar arguments in [21], the weighting factors ωℓ can be determined as
(18)ωℓ(x,y)=(Tcℓ)2|ψℓ±(x,y)|2∑m=1NE(Tcm)2|ψm±(x,y)|2,ℓ=1,…,NE,
where ψℓ± in (13) is the inverse fast Fourier transform of the measured k-space data Sℓ±.


### 2.4. Recovery of Internal Current Density Using the Optimized∇~2Bz

The internal current densityJ=-σ∇u and the magnetic flux density B=(Bx,By,Bz) in Ω satisfy the Ampère law J=∇×B/μ0 where μ0 is the magnetic permeability of the free space. The magnetic resonance current density imaging (MRCDI) technique, which allows the rotation of the object in the MRI scanner, directly visualizes the internal current density by measuring the full components of B [31].The MREIT techniques focus on visualizing the internal current density using onlyBz component of B without rotating the subject. A cylindrical imaging domain Ω can be represented as
(19)Ω=⋃t∈(-H,H)‍Ωt,whereΩt=Ω∩{(x,y,z)∈ℝ3∣z=t},
where Ω0 denotes the middle slice of the imaging subject Ω.In the paper [27], the only recoverable current from the measured Bz data can be represented as JP=J0+J*, where J0=∇α and J*=(∂β/∂y,-∂β/∂x,0). Here, α is a homogeneous voltage potential satisfying
(20)∇2α=0inΩ,∇α·ν=J·νon∂Ω,∫∂Ω‍αds=0,
and βt(x,y):=β(x,y,t) satisfies the following two-dimensional Laplace equation for each slice Ωt⊂Ω:
(21)∇~2βt=1μ0∇2BzinΩt,βt=0on∂Ωt,
where ∇=(∂/∂x,∂/∂y,∂/∂z) and ∇~=(∂/∂x,∂/∂y). From the optimized ∇~2Bz in (16) on each imaging slice Ωt, we can estimate ∇2Bz in (21).Equations (20) and (21) show that we can reconstruct the projected current JP from the optimized ∇~2Bz immediately, instead of Bz, by solving two-dimensional Laplace equations in the region of interest (ROI). The projected current JP provides an optimal approximation of the true current J and, moreover, the gap J-JP depends only on the longitudinal component Jz-Jz0 of  J-J0.


### 2.5. Experimental Setup

In order to demonstrate the proposed method, we performed a phantom with a saline solution including a balloon for the visualization of internal current density. The internal of the balloon was filled with the same saline solution and the volume of the balloon was controlled by injecting the saline solution, which excluded other artifacts by any concentration gradient in the phantom. Figure2(a) illustrates the used balloon for the phantom experiment, and Figures 2(b) and 2(c) show a phantom design to describe how to setup the balloon phantom.(a) Balloon used for the experiment, (b) and (c) balloon phantom design and the electrodes position, respectively.
(a)
(b)
(c)After positioning the phantom inside a 3.0T MRI scanner (Achieva, Philips), we collectedk-space data with 8-channel RF coil using the gradient multi-echo ICNE pulse sequence, which extends throughout the duration of the injection current until the end of a readout gradient [20]. The maximum amplitude of the injection current was 5 mA and the total imaging time was 12.36 seconds to measure the interleaved k-space Sℓ± data, ℓ=1,…,NE. The slice thickness was 5 mm, the number of axial slices was one, the repetition time TR=60 ms, the echo spacing ΔTE=6 ms, the flip angle was 40 degree, and the multi-echo time TEℓ=6+(ℓ-1)×6 ms for NE=9. The FOV was 160×160 mm2 with a matrix size of 128×128. The current injection time Tcℓ for each echo was almost the same as the multi-echo time TEℓ=6+(ℓ-1)×6,  ℓ=1,…,9 because the current was continuously injected until the end of the readout gradient.


## 2.1. Measurement ofBz Using Interleaved ICNE Acquisition

For the interleaved ICNE acquisition method, we inject the positive and negative currents,I+ and I-, through the attached electrodes by scanning each phase encoding consecutively. For a standard spin echo pulse sequence without current injection, the k-space MR signal can be expressed as
(2)S(n,m)=∫ℝ2‍ρ(x,y)eiδ(x,y)e-i2π(Δkxnx+Δkymy)dxdy,
where ρ(x,y) is the real transverse magnetization, δ denotes the phase artifact of background field inhomogeneity, and Δkx and Δky are the reciprocals of fields of view for the x direction and y direction, respectively. During the data acquisition, we set
(3)Δkx=γ2πGxΔt,Δky=γ2πΔGyTpe,
and sample the data in (2) finitely for
(4)n=-Nx2,…,Nx2-1,m=-Ny2,…,Ny2-1.
The names of the above parameters are  γ=26.75×107rad/T·s: the gyromagnetic ratio of the proton,Gx: the magnetic reading gradient strength,Δt: the sampling interval,ΔGy: the phase encoding step,Tpe: the phase encoding time interval,Nx: the number of sampling points,Ny: the number of phase encoding lines.For the conventional MREIT case, we inject the current for the duration ofTc0 from the end of the 90∘ RF pulse to the beginning of the reading gradient. In this case, the induced magnetic flux density Bz due to the injection current provides the additional dephasing of spins and consequently the extra phase is accumulated during Tc0. The corresponding k-space data for the injection currents I± can be represented as
(5)SI±(n,m)=∫ℝ2‍ρ(x,y)eiδ(x,y)e±iγTc0Bz(x,y)×e-i2π(Δkxnx+Δkymy)dxdy.
Using the notations
(6)ρ±(x,y)=ρ(x,y)eiδ(x,y)e±iγTc0Bz(x,y),
we can compute the magnetic flux density Bz as
(7)Bz(x,y)=12γTc0tan-1(α(x,y)β(x,y)),
where α and β are the imaginary and real parts of ρ+/ρ-, respectively.In the conventional MREIT case, the noise standard deviation of the measuredBz, sdBz, is given as a known quantity, which is inversely proportional to the current injection time Tc and the SNR of the MR magnitude image Υ as follows [28, 29]:
(8)sdBz(r)=12γTcΥ(r).
Since the ICNE MR pulse sequence injects the current until the end of a read gradient, the total current injection time of the ICNE case is TcT:=Tc0+Ts and the k-space data is represented as
(9)SC±(n,m)=∫ℝ2‍ρ(x,y)eiδ(x,y)e±iγ(Tc0+Ts/2)Bz(x,y)×e-i2π(Δkxnx+Δkymy)dxdy,
where Ts=ΔtNx is the data acquisition time. In the usual spin echo, the ICNE current injection method demonstrates better SNR in the measured magnetic flux density data than the conventional current injection method. The optimal data acquisition time Ts* has been calculated for the usual spin echo as
(10)Ts*=233TE-3τrf,
which optimally reduces the noise in the Bz data, where τrf is the time of RF pulse [30].In the ICNE MR pulse sequence case, the noise standard deviation of the measuredBz,  sdBzICNE, is given as follows [30]:
(11)sdBzICNE(r)=12γ(Tc0+Ts*/2)Υ(r).
The prolonged data acquisition time, however, may suffer from undesirable side artifacts such as blurring, chemical shift, and motion artifacts along the phase encoding direction. To reduce the undesirable side artifacts, we divide the prolonged data acquisition time into several short ones in the ICNE-multi-echo MR pulse sequence.


## 2.2. Measurement of∇~2Bz Using ICNE-Multiecho Train

Using the ICNE-multi-echo MR pulse sequence, the measuredk-space data can be represented as
(12)Sℓ±(n,m)=∫ℝ2ρℓ(x,y)eiδℓ(x,y)e±iγTcℓBz(x,y)×e-i2π(Δkxnx+Δkymy)dxdy,ℓ=1,…,NE,
where NE is the echo number, Tcℓ is the ℓth time width of the injected current, and ρℓ and δℓ denote the ℓth transverse magnetization and phase artifact, respectively.Figure1 presents a schematic diagram for the ICNE-multi-echo MR pulse sequence based on a gradient echo pulse sequence. By taking the inverse fast Fourier transform, the ICNE-multi-echo sequence generates multiple complex images with different magnitude amplitudes depending on T2* decay and different widths of current injection time:
(13)ψℓ±(x,y):=ρℓ(x,y)eiδℓ(x,y)e±iγTcℓBz(x,y),ℓ=1,…,NE.
Using the relation (13), we derive a formula for ∇~2Bzℓ as
(14)∇~2Bzℓ(x,y)=1iγTcℓ∇~·(|ψℓ+(x,y)|ψℓ+(x,y)∇~ψℓ+(x,y)|ψℓ+(x,y)|-|ψℓ-(x,y)|ψℓ-(x,y)∇~ψℓ-(x,y)|ψℓ-(x,y)|),
where ∇~Bzℓ=(∂Bzℓ/∂x,∂Bzℓ/∂y) denotes the two-dimensional gradient of Bzℓ. The induced ∇~2Bzℓ in (14) removes the low-frequency phase artifact δℓ by subtracting ∇~2δℓ in (14).Figure 1
Diagram of the ICNE-multi-echo MR pulse sequence based on a gradient echo.The calculated vector(|ψℓ±|/ψℓ±)∇~(ψℓ±/|ψℓ±|) corresponding to iγTcℓ∇~(δℓ±Bzℓ) includes unavoidable measured noise. When we consider the decomposed form of (|ψℓ|/ψℓ)∇~(ψℓ/|ψℓ|)=∇~f+∇~×Ψ, where the curl term ∇~×Ψ is a part of unavoidable measured noise, the divergence procedure for ∇~2Bzℓ in (14) cancels ∇~×Ψ, and therefore the measured ∇~2Bzℓ includes a denoising procedure by suppressing a part of the measured noise.


## 2.3. Optimal Combination of Measured∇~2Bzℓ,  ℓ=1,…,NE

The measured∇~2Bzℓ,  ℓ=1,…,NE includes different amounts of unavoidable noise since the intensity of transverse magnetization and the width of injected current are different at each echo.The noise standard deviation of∇~2Bzℓ in (14) is given as
(15)sd∇~2Bzℓ(x,y)=CγTcℓΥℓ(x,y),
where the constant C only relates to the numerical differentiations for ∇~2Bzℓ, and Υℓ denotes the SNR of the ℓth MR magnitude image.Since the noise levels of the measured∇~2Bzℓ,  ℓ=1,…,NE in (14) are given as known quantities, we can utilize the known information sd∇~2Bzℓ to determine an optimized ∇~2Bz which combines the multiple ∇~2Bzℓ:
(16)∇~2Bz(x,y)=∑ℓ=1NE‍ωℓ(x,y)∇~2Bzℓ(x,y).
The problem of determining the weighting factors ωℓ for ∇~2Bzℓ can be formulated as
(17)minωℓ(x,y),ℓ=1,…,NE∑ℓ=1NE‍(ωℓ(x,y))2Var∇~2Bzℓ(x,y)subjectto∑ℓ=1NE‍ωℓ(x,y)=1,ωℓ(x,y)>0,
where Var∇~2Bzℓ denotes the noise variance of ∇~2Bzℓ,  ℓ=1,…,NE in (14).Following similar arguments in [21], the weighting factors ωℓ can be determined as
(18)ωℓ(x,y)=(Tcℓ)2|ψℓ±(x,y)|2∑m=1NE(Tcm)2|ψm±(x,y)|2,ℓ=1,…,NE,
where ψℓ± in (13) is the inverse fast Fourier transform of the measured k-space data Sℓ±.


## 2.4. Recovery of Internal Current Density Using the Optimized∇~2Bz

The internal current densityJ=-σ∇u and the magnetic flux density B=(Bx,By,Bz) in Ω satisfy the Ampère law J=∇×B/μ0 where μ0 is the magnetic permeability of the free space. The magnetic resonance current density imaging (MRCDI) technique, which allows the rotation of the object in the MRI scanner, directly visualizes the internal current density by measuring the full components of B [31].The MREIT techniques focus on visualizing the internal current density using onlyBz component of B without rotating the subject. A cylindrical imaging domain Ω can be represented as
(19)Ω=⋃t∈(-H,H)‍Ωt,whereΩt=Ω∩{(x,y,z)∈ℝ3∣z=t},
where Ω0 denotes the middle slice of the imaging subject Ω.In the paper [27], the only recoverable current from the measured Bz data can be represented as JP=J0+J*, where J0=∇α and J*=(∂β/∂y,-∂β/∂x,0). Here, α is a homogeneous voltage potential satisfying
(20)∇2α=0inΩ,∇α·ν=J·νon∂Ω,∫∂Ω‍αds=0,
and βt(x,y):=β(x,y,t) satisfies the following two-dimensional Laplace equation for each slice Ωt⊂Ω:
(21)∇~2βt=1μ0∇2BzinΩt,βt=0on∂Ωt,
where ∇=(∂/∂x,∂/∂y,∂/∂z) and ∇~=(∂/∂x,∂/∂y). From the optimized ∇~2Bz in (16) on each imaging slice Ωt, we can estimate ∇2Bz in (21).Equations (20) and (21) show that we can reconstruct the projected current JP from the optimized ∇~2Bz immediately, instead of Bz, by solving two-dimensional Laplace equations in the region of interest (ROI). The projected current JP provides an optimal approximation of the true current J and, moreover, the gap J-JP depends only on the longitudinal component Jz-Jz0 of  J-J0.


## 2.5. Experimental Setup

In order to demonstrate the proposed method, we performed a phantom with a saline solution including a balloon for the visualization of internal current density. The internal of the balloon was filled with the same saline solution and the volume of the balloon was controlled by injecting the saline solution, which excluded other artifacts by any concentration gradient in the phantom. Figure2(a) illustrates the used balloon for the phantom experiment, and Figures 2(b) and 2(c) show a phantom design to describe how to setup the balloon phantom.(a) Balloon used for the experiment, (b) and (c) balloon phantom design and the electrodes position, respectively.
(a)
(b)
(c)After positioning the phantom inside a 3.0T MRI scanner (Achieva, Philips), we collectedk-space data with 8-channel RF coil using the gradient multi-echo ICNE pulse sequence, which extends throughout the duration of the injection current until the end of a readout gradient [20]. The maximum amplitude of the injection current was 5 mA and the total imaging time was 12.36 seconds to measure the interleaved k-space Sℓ± data, ℓ=1,…,NE. The slice thickness was 5 mm, the number of axial slices was one, the repetition time TR=60 ms, the echo spacing ΔTE=6 ms, the flip angle was 40 degree, and the multi-echo time TEℓ=6+(ℓ-1)×6 ms for NE=9. The FOV was 160×160 mm2 with a matrix size of 128×128. The current injection time Tcℓ for each echo was almost the same as the multi-echo time TEℓ=6+(ℓ-1)×6,  ℓ=1,…,9 because the current was continuously injected until the end of the readout gradient.


## 3. Results

Figure3(a) shows the acquired magnitude images |ρℓ|,  ℓ=1,…,9, where ρℓ was the ℓth measured T2* weighted complex image. Figure 3(b) shows the measured ∇~2Bzℓ images using (14) corresponding to the ℓth k-space data Sℓ±,  ℓ=1,…,9. Inside and outside of the balloon, the MR magnitude images are almost the same because of the same saline solution, but the measured ∇~2Bzℓ images show distinguishable signals reflecting the conductivity changes inside and outside of the balloon.(a) Acquired magnitude images|ρℓ|,  ℓ=1,…,9, where ρℓ was the ℓth measured T2* weighted complex image, (b) measured ∇~2Bzℓ images using (14) corresponding to the ℓth k-space data Sℓ±,  ℓ=1,…,9.
(a)
(b)Since both sides, inside and outside of the balloon, are homogeneous, the∇~2Bzℓ≈-μ0∇u×∇σ should be near zero except the boundary of the balloon without noise effect because the conductivity value is constant in each region. To evaluate the noise level of Bzℓ, we calculated the discrete L2-norm:
(22)Err2∶=∥∇~2Bzℓ∥L2(Ω∖∂D)=∑(xi,yj)∈Ω∖∂D‍∇~2Bzℓ(xi,yj)2|Ωij|,
where Ω is the imaging ROI region, ∂D denotes the boundary of balloon, and |Ωij| is the pixel size.Table1 shows the L2-norm, Err2, in which the values depends on the T2* decay rate and the width of injected current.Table 1
Noise level estimation of the measured∇~2Bzℓ, ℓ=1,…,9, by calculating (22).
1st
2nd
3rd
4th
5th
6th
7th
8th
9th
Err
2
129.1
142.1
115.0
110.4
203.1
161.7
195.6
183.5
166.6The estimated noise levels were reduced up to the 4th echo, but increased in the following echoes because the intensity of magnitude images follows the exponentialT2* decay, and the width of the injected current linearly increases.Figures4(a) and 4(b) show the measured ∂Bzℓ/∂x and ∂Bzℓ/∂y images, respectively, where
(23)∇~Bzℓ(x,y)=|ψℓ+(x,y)|ψℓ+(x,y)∇~ψℓ+(x,y)|ψℓ+(x,y)|-|ψℓ-(x,y)|ψℓ-(x,y)∇~ψℓ-(x,y)|ψℓ-(x,y)|,ℓ=1,…,9.
Since the currents were transversally injected, the measured ∂Bzℓ/∂y reflected dominant internal current flows.(a) Measured∂Bzℓ/∂x images, (b) measured ∂Bzℓ/∂y images, ℓ=1,…,9.
(a)
(b)Figures5(a) and 5(b) display the reconstructed ∇~2Bzavrg=(1/NE)∑ℓ=1NE‍∇~2Bzℓ and ∇~2Bzopt=∑ℓ=1NE‍ωℓ∇~2Bzℓ images, respectively, where ωℓ is the weighting factor by solving (17).(a) Reconstructed∇~2Bzavrg=(1/NE)∑ℓ=1NE‍∇~2Bzℓ images, (b) reconstructed ∇~2Bzopt=∑ℓ=1NE‍ωℓ∇~2Bzℓ images, where ωℓ is the weighting factor by solving (17).
(a)
(b)Figure6 shows the recovered current density images, Jxℓ and Jyℓ, corresponding to the ℓth echo. To obtain the current density images, we solved (20) for the background homogeneous current and the two-dimensional harmonic equation (21) to reflect the measured ∇~2Bzℓ data.Recovered current density by solving (20) and (21). (a) Recovered Jxℓ images, (b) recovered Jyℓ images, ℓ=1,…,9.
(a)
(b)We recovered the current densityJopt by solving (20) and (21) using the optimized ∇~2Bzopt=∑ℓ=1NE‍ωℓ∇~2Bzℓ, where ωℓ is the weighting factor by solving (17). The recovered Jxopt and Jyopt are displayed in Figure 7.Recovered current density by solving (20) and (21). (a) Reconstructed Jxopt images, (b) Jyopt images by using Bzopt.
(a)
(b)Table2 shows the estimated noise level of the recovered ∇~2Bzavrg and ∇~2Bzopt, by calculating (22). The estimated noise levels validate the proposed method because the inside and outside of the balloon in the phantom should be homogeneous.Table 2
Noise level estimation of the recovered∇~2Bzavrg and ∇~2Bzopt, by calculating (22).
∇
~
2
B
z
avrg
=
(
1
/
N
E
)
∑
ℓ
=
1
N
E
∇
~
2
B
z
ℓ
∇
~
2
B
z
opt
=
∑
ℓ
=
1
N
E
ω
ℓ
∇
~
2
B
z
ℓ
Err
2
145.1
41.5


## 4. Discussion

Since the MREIT technique conventionally used the interleaved phase encoding acquisition scheme to measure the magnetic flux density by alternating two currents with positive and negative polarities, we could obtain the coil sensitivity information without additional scans by product ofψcj+ and ψcj-:
(24)Ψcj(x,y)∶=ψcj+(x,y)ψcj-(x,y)=(ρcj(x,y))2e2iδcj(x,y),j=1,…,NC,
where δci is the jth coil sensitivity and NC is the number of coils. For a fast MRI, using the a priori spatial information from the multiple receiver coils, the sensitivity encoding (SENSE) technique enables one to reduce the number of Fourier encoding steps while preserving the spatial resolution [32]. For a temporal variation of the internal conductivity, if we estimate the reference coil sensitivity using (24), which is independent of the injected current, the SENSE technique can be applicable to the proposed method to visualize the internal current density combining the multi-echo train.In this paper, we directly measure∇~2Bz, which is sufficient to reconstruct the internal current density using the injected current information. The proposed method to measure ∇~2Bz in (14) can avoid a tedious unwrapping procedure. The proposed method may exhibit potential to be applied for conventional phase imaging techniques.The optimal combination of multiple echoes by determining optimal weighting factor in (17) effectively reduces the noise level of measured ∇~2Bz. Since the decay rate of magnitude and the width of injected current can be determined pixel by pixel, we can determine a pixel-wise noise level of the optimized ∇~2Bz data. Since most algorithms for the MREIT technique visualize the internal conductivity and/or current density in an entire imaging region due to the relationships between the external injection current and the internal measured magnetic flux density data, the estimated noise level of Bz can be used to determine the denoising level of the measured data in defective regions.To optimize the multiple echoes, we consider only the uniformly distributed random noise effect, but unavoidable spike or different nonuniform noise may deteriorate a combined measured data. Thus it is important to develop a method to discard the non-uniform noises in the optimizing process in order to enhance the quality ofBz.Our future studies will focus on reducing the imaging time with a feasible noise level to produce conductivity images for the application of functional MREIT imaging to animal brains in order to visualize the rapidly changing conductivity associated with neural activation.


## 5. Conclusion

We have visualized the internal current density using a fast ICNE-multi-echo MR pulse sequence based on a gradient echo by two measurements in the interleaved acquisition. The interleaved acquisition method in MREIT is a conventional method to suppress the background field inhomogeneity phase artifact and to increase the SNR ofBz by doubling the accumulated phase signal. We used the multi-echo pulse sequence, which acquires multiple sampling points within each repetition time. The proposed method directly measures the Laplacian of Bz from the measured k-space data, which can avoid a tedious unwrapping procedure and include a denoising effect by removing a part of the measured noise. We determined an optimal combination of the magnetic flux densities from the multi-echo in order to reduce the noise level. Using the optimization of ∇2Bz, the proposed method visualized the internal current density using the relationships between the induced internal current and the measured ∇2Bz data, while suppressing the background field inhomogeneity. A real phantom experiment with a saline solution including a balloon was carried out to verify that the proposed method can be feasibly applied in real experiments. The total scan time in the phantom experiment was less than 13 seconds to visualize the current density with a 128×128 spacial matrix size.



---

*Source: 381507-2013-03-20.xml*