CoReHA 2.0: A Software Package for In Vivo MREIT Experiments

Abstract

Magnetic resonance electrical impedance tomography (MREIT) is a new medical imaging modality visualizing static conductivity images of electrically conducting subjects. Recently, MREIT has rapidly progressed in its theory, algorithm, and experiment technique and now reached to the stage of in vivo animal experiments. In this paper, we present a software, named CoReHA 2.0 standing for the second version of conductivity reconstructor using harmonic algorithms, to facilitate in vivo MREIT reconstruction of conductivity image. This software offers various computational tools including preprocessing of MREIT data, identification of 2D geometry of the imaging domain and electrode positions, and reconstruction of cross-sectional scaled conductivity images from MREIT data. In particular, in the new version, we added several tools including ramp-preserving denoising, harmonic inpainting, and local harmonic 𝐵_𝑧 algorithm to deal with data from in vivo experiments. The presented software will be useful to researchers in the field of MREIT for simulation, validation, and further technical development.

1. Introduction

Recently, a new imaging modality called magnetic resonance electrical impedance tomography (MREIT) has been introduced, which allows high resolution imaging of tomographic electrical conductivity distributions of biological objects [1, 2]. The technique involves (i) current injection into an electrically conducting object such as animal or human body through surface electrodes, (ii) measurement of induced internal magnetic flux density using an MRI system, typically only the 𝑧-component 𝐵_𝑧 of the induced magnetic flux density 𝐁=(𝐵_𝑥,𝐵_𝑦,𝐵_𝑧), where 𝑧 is the axis parallel to the main magnetic field of the MR scanner, and (iii) conductivity reconstruction by solving nonlinear boundary value problems with given injected currents and measured magnetic flux density, employing finite element methods. Although it looks straight forward to reconstruct conductivity distribution as described in the three steps, computation involves several innovative approaches including magnetic flux density estimation, data verification, segmentation, and solving forward/inverse problems. These specific computations cannot be handled through readily available finite element packages; thus, it is required to develop a user friendly software with graphics user interface (GUI) for those who wish to reconstruct conductivity distributions. Upon these requests, we developed a software package, called CoReHA [3, 4] which stands for conductivity reconstructor using harmonic algorithms, using VC++ MFC 6.0 (Microsoft Foundation Class Library 6.0) and OpenGL under the Microsoft Windows operating system. Based on the harmonic 𝐵_𝑧 algorithm [5], CoReHA supports all procedures from the preprocessing of raw data to conductivity imaging through the intuitively apprehensible graphic user interface, more specifically, data conversion of raw 𝑘-space data, data verification, segmentation tools for numerical computation, solvers of forward/inverse problems using finite element methods, and 2D/3D data view as well as histogram. This software has been a major tool to facilitate multilateral studies for MREIT and has brought out successful reconstruction of conductivity imaging in many researches [4, 6–8].

In the stage of in vivo MREIT experiments, we have to address the improvement of the signal-to-noise ratio (SNR) of measured 𝐵_𝑧 data, since 𝐵_𝑧 data has weak strength due to the low amount of injected current for the safety guide. Moreover, there may exist MR signal void regions in the animal or human body, where noise levels of 𝐵_𝑧 are excessively high, and as a result, uncertain effects are produced. For dealing with these technical problems, several methods were developed, called ramp preserving denoising [9], harmonic inpainting [10], and local harmonic 𝐵_𝑧 algorithm [11]. Incorporating these new features, we release the second version of CoReHA for in vivo MREIT experiments. The presented software will be useful to researchers in the field of MREIT experimental studies as well as in vivo animal/human experiments. CoReHA is available from the website http://iirc.khu.ac.kr/.

This paper is organized as follows. In Section 2, MREIT system is introduced and related works are explained. Brief explanation of the previous version of this software is given for comprehensive understanding of our works in Section 3. Details of new tools of the software for dealing with in vivo stage issues are provided in Section 4. Finally, we close in Section 5 with conclusions and a discussion of future works.

2. The Basic of MREIT and Related Works

We briefly explain how MR scanner is used as a tool to capture internal magnetic flux density images, and the MR scanner has its magnetic field in 𝑧-direction. Let an electrically conducting subject occupy a three-dimensional domain Ω with its boundary 𝜕Ω and a conductivity distribution 𝜎. As shown in Figure 1, we attach pairs of electrodes ℰ⁺_𝑗 and ℰ⁻_𝑗 along 𝜕Ω in order to inject a current 𝐼_𝑗 in a form of pulses whose timing is synchronized with an MR pulse sequence for 𝑗=1,2. The injection current produces current density 𝐉_𝑗 satisfying the following elliptic equations: ∇⋅𝜎∇𝑢_𝑗[𝜎]=0inΩ,𝐼=_ℰ+𝑗𝜎𝜕𝑢_𝑗[𝜎]𝜕𝐧𝑑𝑠=−_ℰ−𝑗𝜎𝜕𝑢_𝑗[𝜎]𝜕𝐧𝑑𝑠,∇𝑢_𝑗[𝜎]×𝐧|_ℰ⁺_𝑗∪ℰ⁻_𝑗𝜎=0,𝜕𝑢_𝑗[𝜎]𝜕𝐧=0on𝜕Ω⧵ℰ⁺_𝑗∪ℰ⁻_𝑗, (1) where 𝐉_𝑗=−𝜎∇𝑢_𝑗. Also, magnetic flux density 𝐁 is given by Biot-Savart law as follows: 𝐁(𝐱)=𝜇₀4𝜋_Ω𝐉𝐱^′×𝐱−𝐱^′𝐱−𝐱^′3𝑑𝐱^′+ℋ(𝐱), (2) where ℋ is magnetic flux density from electrodes, wires, and others. Note that ∇²ℋ=0 in Ω. Then MR spectrometer provides the complex 𝑘-space data 𝒮 that is influenced by 𝐵_𝑧 in the following way. Given current injection 𝐼_𝑗, we have 𝒮_𝑗(𝑚,𝑛,𝑧)=𝑀𝑥,𝑦,𝑧𝑒^{𝑖𝛿(𝑥,𝑦,𝑧)}𝑒^{𝑖𝛾𝐵𝑧,𝑗(𝑥,𝑦,𝑧)𝑇𝑐}×𝑒^{−𝑖(𝑥𝑚Δ𝑘𝑥+𝑦𝑛Δ𝑘𝑦)}𝑑𝑥𝑑𝑦, (3) where 𝑀 is a conventional MR magnitude image, 𝛿 any systematic phase artifact, 𝛾=26.75×10⁷ rad/T·s the gyromagnetic ratio of hydrogen, and 𝑇_𝑐 the current pulse width in seconds. By the discrete inverse Fourier transformation of the 𝑘-space data 𝒮 in (3), we obtain the following complex images: ℳ⁺𝑥,𝑦,𝑧=𝑀𝑥,𝑦,𝑧𝑒^{𝑖𝛿(𝑥,𝑦,𝑧)}𝑒^{𝑖𝛾𝐵𝑧(𝑥,𝑦,𝑧)𝑇𝑐}. (4) To extract 𝐵_𝑧 data stably while eliminating systematic phase artifact 𝛿, we inject the counter directional current providing the counter part of ℳ⁺ in (4) so that we get ℳ⁻𝑥,𝑦,𝑧=𝑀𝑥,𝑦,𝑧𝑒^{𝑖𝛿(𝑥,𝑦,𝑧)}𝑒^{−𝑖𝛾𝐵𝑧(𝑥,𝑦,𝑧)𝑇𝑐}. (5) Division of ℳ⁺ by its counter part ℳ⁻ leads to 𝐵_𝑧=1𝑥,𝑦,𝑧2𝛾𝑇_𝑐ℳarg⁺𝑥,𝑦,𝑧ℳ⁻.𝑥,𝑦,𝑧 (6)

Figure 1 

Domain Ω for MREIT.

MRCDI [12] and early stage MREIT [13–15] use measurements of all three components of 𝐁=(𝐵_𝑥,𝐵_𝑦,𝐵_𝑧) which require subject rotations inside MR scanner. However, experiences show that these subject rotations are impractical and also cause other problems such as misalignments of pixels. Hence, in order to make the MREIT technique easily applicable to clinical situations, we should use only 𝐵_𝑧 data for the conductivity reconstruction.

In the harmonic 𝐵_𝑧 algorithm [5, 16, 17] which is the first constructive 𝐵_𝑧-based MREIT algorithm, we inject two independent electrical currents 𝐼₁ and 𝐼₂ through two pairs of surface electrodes ℰ^±₁ and ℰ^±₂, respectively, to get two raw 𝑘-space data. The amount of each injection current is determined by the product of current amplitude and time duration. A modified spin-echo pulse sequence synchronized with the current injection is typically used as described in experimental works [5, 18, 19], and 𝐵_𝑧,1 and 𝐵_𝑧,2 are obtained from the formula (6). Then the reconstruction algorithm is based on the following identity: 𝜕ln𝜎(𝐱)𝜕𝑥𝜕ln𝜎(𝐱)=1𝜕𝑦𝜇₀(𝔸[𝜎](𝐱))⁻¹∇²𝐵_𝑧,1(𝐱)∇²𝐵_𝑧,2(𝐱),𝐱∈Ω, (7) where 𝔸[𝜎](𝐱)=−𝐽_𝑦,1𝐽_𝑥,1−𝐽_𝑦,2𝐽_𝑥,2=𝜎𝜕_𝑦𝑢₁−𝜎𝜕_𝑥𝑢₁𝜎𝜕_𝑦𝑢₂−𝜎𝜕_𝑥𝑢₂,𝐱∈Ω (8) and 𝐽_𝑥,𝑗 and 𝐽_𝑦,𝑗 are the 𝑥 and 𝑦 components of induced current density 𝐉_𝑗, respectively. We now present the conductivity reconstruction procedure using the harmonic 𝐵_𝑧 algorithm in CoReHA, which consists of six steps.

Step 1. We impress electrical currents 𝐼₁ and 𝐼₂ through pairs of surface electrodes ℰ^±₁ and ℰ^±₂, respectively, and get the 𝑘-space data 𝒮_𝑗 for 𝑗=1,2 using an MR scanner.

Step 2. We produce an MR magnitude image 𝑀 and induced magnetic flux densities 𝐵_𝑧,𝑗 from the 𝑘-space data 𝒮_𝑗 for 𝑗=1,2 given by (6).

Step 3. Using the MR magnitude image 𝑀, we perform segmentation of 𝜕Ω and ℰ^±_𝑗. Here, we use level-set-based segmentations.

Step 4. We set the initial guess 𝜎⁰(𝐱,𝑧)=1 in Ω.

Step 5. We solve the problem (1) with 𝜎=𝜎⁰.

Step 6. We solve ∇²_𝑥𝑦ln𝜎¹=∇_𝑥𝑦⋅𝔸𝜎⁰⁻¹∇²𝐵_𝑧,1∇²𝐵_𝑧,2inΩ_𝑧 (9) with the boundary condition ∇_𝑥𝑦ln𝜎¹⋅𝜈=𝔸𝜎⁰⁻¹∇²𝐵_𝑧,1∇²𝐵_𝑧,2⋅𝜈on𝜕Ω_𝑧, (10) where ∇²_𝑥𝑦 and ∇_𝑥𝑦 are two-dimensional Laplacian and gradient operators on the 𝑥𝑦-plane, respectively, Ω_𝑧 a two-dimensional slice of Ω, which is perpendicular to the 𝑧-axis, and 𝜈 the two-dimensional outward normal vector to 𝜕Ω_𝑧.

Note that in the harmonic 𝐵_𝑧 algorithm in [5, 16], Steps 5 and 6 are repeated to get 𝜎^𝑛 from 𝜎^𝑛−1 until it converges. In order to ensure the convergence of the harmonic 𝐵_𝑧 algorithm, the conductivity values of the subject’s boundary should be homogeneous and the contrast of the conductivity distribution is assumed small enough [16, 17]. But, in general, the conductivity distribution of animal and human is quite inhomogeneous [11]. Also, for the completion of the iterative procedure, we have to resolve heavy cost computational issues for 3D forward problem (1), related to, for example, 3D segmentation, 3D mesh generation, and 3D forward problem solver. Fortunately, for the clinical purpose, it is enough to find the scaled conductivity images instead of the true ones because the reconstructed scaled ones reflect the fine details of the true conductivity contrast when the conductivity contrast is low [11]. Hence, in CoReHA, we iterate Steps 5 and 6 only once. Due to the same reason, in Step 5, we solve a simplified problem ∇_𝑥𝑦⋅𝜎∇_𝑥𝑦𝑢_𝑗[𝜎]=0inΩ_𝑧𝑢_𝑗=1onℰ⁺_𝑗∩Ω_𝑧𝑢_𝑗=−1onℰ⁻_𝑗∩Ω_𝑧𝜎𝜕𝑢_𝑗[𝜎]𝜕𝜈=0on𝜕Ω_𝑧⧵ℰ⁺_𝑗∪ℰ⁻_𝑗, (11) instead of (1). In addition, in Step 6, we adopt the boundary condition 𝜎¹=1on𝜕Ω_𝑧, (12) to replace (10). Note that three-dimensional Laplacian of 𝐵_𝑧 may not be invalid due to the circumstance of the experiments. Hence, CoReHA supports for both 2D and 3D Laplacian of 𝐵_𝑧.

3. CoReHA 1.0: For Multilateral Studies

CoReHA implements Step 2 through 6 of the harmonic 𝐵_𝑧 algorithm since Step 1 is about the experiment using an MR scanner and EIT equipment. In this section, we briefly describe each step of CoReHA 1.0 for understanding the software package on the whole.

After an MREIT experiment, one can perform the phase extraction process by applying fast Fourier transform in Section 2 as well as obtain the MR imaging by taking magnitude of ℳ^±. For this computation, we use FFTW [20], one of the well-known fast Fourier transform libraries, for obtaining (4) from (3). Since 𝐵_𝑧 in (6) is wrapped due to the branch cut of argument operator, we apply the Goldstein’s algorithm [21] for two-dimensional phase unwrapping. Note that even if 𝐵_𝑧 data is continuous in the 𝑥𝑦-plane given by phase unwrapping, we do not guarantee the continuity along the 𝑧-direction. In order to verify the continuity of 𝐵_𝑧 along the 𝑧-direction, CoReHA supports a verification tool. Hence, if the verification fails, one should drop (𝜕²/𝜕𝑧²)𝐵_𝑧 in ∇²𝐵_𝑧 for preventing the artifact coming from wrong Laplacian calculation.

In order to get the imaging domain, we have employed statistically reinstating method (SRM) [22] on MR image, which is a level-set-based segmentation method. After the boundary of the object including electrodes is segmented by SRM, one can operate additional manual work to configure the object without electrodes or to modify the local geometry. After segmentation, one can generate the triangle mesh for the numerical computation. For the triangulation, we adopt the Triangle [23] which is a well-known open source software for the two-dimensional triangulation.

For the computation of (11) and (9), we implement the standard 𝒫₁ finite element method and apply the conjugate gradient (CG) method for matrix inversion. The size of generated triangulation is automatically determined by given geometrical information and its size is small enough to cover the pixel size. Of course, one can adjust the size of triangulation for one’s own purpose.

For the details of CoReHA 1.0, see [3].

4. CoReHA 2.0: For Better Imaging Quality

In this section, we describe in detail new features in the second version of CoReHA for dealing with in vivo MREIT experiments.

4.1. Preprocessing: Ramp Preserving Denoising

For the safety of the subject in in vivo experiments, the amount of injected current 𝐼 should be reduced less than a few milliamperes. The SNR of 𝐵_𝑧 data is also directly affected by the noise of ℳ. Therefore, a proper denoising method is required for obtaining the reconstructed image with reduced noise artifact. Let us consider the reconstruction identity (7) for the derivation of denoising algorithm. Apparently, the change of the conductivity distribution ∇_𝑥𝑦ln𝜎 is directly proportional to the Laplacian of 𝐵_𝑧 data if the matrix 𝔸[𝜎](𝐱) is invertible. When we apply a denoising algorithm, the structure of ramp should be preserved in order to prevent the change of shape or wrong alignment of anomaly location. For our better understanding, we illustrate a simple example using a modified Shepp-Logan phantom. We compute the current density using (11) with the vertically identical conductivity distribution given in Figure 2(a). Applying the Biot-Savart law (2), we obtain the 𝐵_𝑧 data as shown in Figure 2(b). We can see that the ramp structure of 𝐵_𝑧 reflects the change of the conductivity distribution as shown in Figure 2(c). In order to keep the ramp structure during the denoising, we adopt a nonlinear diffusion equation method based on structure tensor [9, 24].

(a)

(b)

(c)

Figure 2 

Explanation of ramp structure. (a) Conductivity distribution. (b) Corresponding 𝐵𝑧 data generated by the current injection along vertical direction. (c) 1D profiles along the line 𝑦=64 which is indicated by the red line in (a) and (b). Green dotted line is the conductivity distribution and the red one is the 𝐵𝑧 data. Locations of conductivity change match the ramp structure which is characterised by changes of slopes of the 𝐵𝑧 data.

For the sake of completeness, we summarise the results in [9, 24]. Let us consider the structure tensor 𝑈 defined by 𝑈∶=²_𝑖=1∇_𝑥𝑦𝑤_𝑖∇_𝑥𝑦𝑤^𝑇_𝑖, (13) where 𝑤₁𝑤₂∶=∇_𝑥𝑦𝐵_𝑧=𝜕_𝑥𝐵_𝑧𝜕_𝑦𝐵_𝑧. (14) To reduce the noise, we solve the following nonlinear PDE cooperating with the structure tensor 𝑈: 𝜕_𝑡𝐵_𝑧(𝐱,𝑡)=∇_𝑥𝑦⋅𝑔𝑈𝐱,𝑇₂∇_𝑥𝑦𝐵_𝑧(𝐱,𝑡)in𝑅×0,𝑇₁,𝑔𝑈𝐱,𝑇₂∇_𝑥𝑦𝐵_𝑧(𝐱,𝑡)⋅𝐧=0on𝜕𝑅×0,𝑇₁,𝐵_𝑧(𝐱,0)=𝐵_𝑧(𝐱)on𝑅, (15) with 𝜕_𝜏𝑢_𝑖𝑗(𝐱,𝜏)=∇_𝑥𝑦⋅𝑔𝑈_𝑠∇_𝑥𝑦𝑢_𝑖𝑗(𝐱,𝜏)in𝑅×0,𝑇₂,𝑔𝑈_𝑠∇_𝑥𝑦𝑢_𝑖𝑗(𝐱,𝜏)⋅𝐧=0on𝜕𝑅×0,𝑇₂,𝑢_𝑖𝑗(𝐱,0)=²_𝑘=1∇_𝑥𝑦𝑤_𝑘(𝐱,𝑡)∇_𝑥𝑦𝑤_𝑘(𝐱,𝑡)^𝑇_𝑖𝑗𝑔(𝑈)=1on𝑅,√𝑣1+Λ^Λ𝑣^𝑇_Λ+1√𝑣1+𝜆^𝜆𝑣^𝑇_𝜆, (16) where 𝑅 is a two-dimensional region for the slice of 𝐵_𝑧 data, 𝑢_𝑖𝑗 is the (𝑖,𝑗) element of 𝑈, Λ and 𝜆 are maximum and minimum eigenvalues of 𝑈, respectively, 𝑣_Λ and 𝑣_𝜆 are corresponding normalized eigenvectors, and 𝑈_𝑠≡𝐺_𝑠∗𝑈 is the (element by element) convolution of 𝑈 with the two-dimensional Gaussian kernel 𝐺_𝑠 with a standard deviation 𝑠. In the new version of CoReHA, the user only considers the parameter 𝑇₁, the total diffusion time of the denoising algorithm. Through extensive numerical studies [9, 24], we figured out that the denoising results are robust to the tensor regularization time 𝑇₂. The noise level of 𝐵_𝑧 affects only 𝑇₁. Consequently, we fix 𝑇₂=2 in this version, which seems to be enough time to regularize the diffusion tensor. Note that the denoising algorithm is only valid when the user considers two-dimensional Laplacian of 𝐵_𝑧 since the diffusion process is on the 𝑥𝑦-plane. For more details, see [9].

4.2. Preprocessing: Harmonic Inpainting Using MR Data

Let us revisit (6). In a local region of Ω, where |ℳ^±|≈0, 𝐵_𝑧 is defected severely by the amplification of noise due to the division procedure. In the previous version of CoReHA to overcome this trouble, a manual segmentation tool was provided for the extraction of the problematic region with the harmonic inpainting algorithm [10] to solve the Poisson equation. However, the manual segmentation causes user-dependent reconstruction results. For dealing with this situation, we assume that the defected region has low MR magnitude and its conductivity distribution is homogenous so that ∇²𝐵_𝑧=0. Then, we automatically select the defected region by taking a threshold for the MR data. Figure 3 shows an example for the selection of the defected region by thresholding. In the new version of CoReHA, we consider the 10% of maximum magnitude of MR data as a threshold. Incorporating the ramp-preserving denoising, we additionally apply isotropic diffusion on the defected region within enough duration. The user does not need to consider how well the defected region is segmented; therefore, it guarantees the robustness of the conductivity reconstruction without user dependency.

(a)

(b)

Figure 3 

Defected region select using threshold (a) MR image and (b) defected region by the thresholding (purple color).

4.3. Reconstruction: Local Harmonic 𝐵_𝑧 Algorithm

We support a modified local harmonic 𝐵_𝑧 algorithm [11] for the successful reconstruction by excluding the defected region directly. Let 𝐷 be a local region without defection which attracts our interests. Then, we locally reconstruct the conductivity by solving (9) and (10) on the local region 𝐷 instead of the whole domain Ω. In order to solve (9) with (10), we use a 𝒫₁ finite element method on a selected triangle mesh and apply conjugate gradient (CG) method for matrix inversion. CoReHA 2.0 supports a user interface to draw and manipulate a polygonal local region by adding or removing points. It is useful when the user wants to see a small local change of conductivity on the local region such as brain or pelvis. In Figure 4(c), we see that the reconstructed conductivity is severely contaminated by locally magnified noise. But, Figure 4(d) shows the successful reconstruction of local brain region by the local harmonic 𝐵_𝑧 algorithm.

(a)

(b)

(c)

(d)

Figure 4 

Local harmonic 𝐵𝑧 algorithm: (a) MR image, (b) selected interested region, (c) conductivity reconstruction on the whole domain, and (d) conductivity reconstruction using the local harmonic 𝐵𝑧 algorithm.

4.4. Additional Tools

Due to the different physical properties between water and fat, the measured MR signal can be shifted in the fat region, since the MR scanner assumes that spin isochromats precess at the Larmor frequency of water. Therefore, the MR signal may become void, and overlapped images are produced in the animal/human experiments. To overcome this difficulty, new version of CoReHA contains the chemical shift correction using 3-points Dixon method [25]. In the 𝑘-space data conversion, one can manipulate the MR signal using the user interface which provides separated water and fat images of MR. For more details, see [25].

Furthermore, we implemented many useful features to help our understanding of MR, 𝐵_𝑧, and conductivity images, for example, magnification windows and advanced alignment and resizing of windows. Figure 5 is a snapshot of CoReHA 2.0, that shows newly introduced features.

Figure 5 

Snapshot of CoReHA 2.0.

5. Conclusions and Future Works

In this paper, we presented new features in a software package, called CoReHA 2.0, to deal with the data from the in vivo stage of MREIT. In this version, we enhanced the preprocessing procedure based on the ramp preserving denoising and harmonic inpainting with MR data to increase SNR of 𝐵_𝑧 data in the in vivo experiment. Also, we provided a user interface for the local harmonic 𝐵_𝑧 algorithm for the reconstruction on the locally interested region. The performance results of suggested algorithms and methods are given in [26]. Note that CoReHA reconstructs the contrast image of conductivity. For the medical diagnosis, the contrast image is enough. However, recently the necessity for the true conductivity came to the fore. For the true conductivity image reconstruction, we need a subject-dependent 3D segmentation and modeling tools, 3D denoising algorithm, and 3D finite element solver and higher order methods for more accurate computation.

Conflict of Interests

The authors declare that this paper has no conflict of interests with CoReHA 2.0 software.

Acknowledgment

This work was supported by NRF-2011-0015399.