CoReHA 2.0: a software package for in vivo MREIT experiments.

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 B z 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 z-component B of the induced magnetic flux density B = (B , B , B ), where z 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 B 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 k-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 B data, since B 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 B 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 B 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 z-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 I in a form of pulses whose timing is synchronized with an MR pulse sequence for j = 1,2. The injection current produces current density J satisfying the following elliptic equations: where J = −σ∇u . Also, magnetic flux density B is given by Biot-Savart law as follows: where ℋ is magnetic flux density from electrodes, wires, and others. Note that ∇2 ℋ = 0 in Ω. Then MR spectrometer provides the complex k-space data 𝒮 that is influenced by B in the following way. Given current injection I , we have where M is a conventional MR magnitude image, δ any systematic phase artifact, γ = 26.75 × 107 rad/T·s the gyromagnetic ratio of hydrogen, and T the current pulse width in seconds. By the discrete inverse Fourier transformation of the k-space data 𝒮 in (3), we obtain the following complex images: To extract B data stably while eliminating systematic phase artifact δ, we inject the counter directional current providing the counter part of ℳ + in (4) so that we get Division of ℳ + by its counter part ℳ − leads to

Figure 1

Domain Ω for MREIT.

MRCDI [12] and early stage MREIT [13-15] use measurements of all three components of B = (B , B , B ) 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 B data for the conductivity reconstruction.

In the harmonic B algorithm [5, 16, 17] which is the first constructive B -based MREIT algorithm, we inject two independent electrical currents I 1 and I 2 through two pairs of surface electrodes ℰ 1 ± and ℰ 2 ±, respectively, to get two raw k-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 B and B are obtained from the formula (6). Then the reconstruction algorithm is based on the following identity: where and J and J are the x and y components of induced current density J , respectively. We now present the conductivity reconstruction procedure using the harmonic B algorithm in CoReHA, which consists of six steps.

Step 1

We impress electrical currents I 1 and I 2 through pairs of surface electrodes ℰ 1 ± and ℰ 2 ±, respectively, and get the k-space data 𝒮 for j = 1,2 using an MR scanner.

Step 2

We produce an MR magnitude image M and induced magnetic flux densities B from the k-space data 𝒮 for j = 1,2 given by (6).

Step 3

Using the MR magnitude image M, we perform segmentation of ∂Ω and ℰ ±. Here, we use level-set-based segmentations.

Step 4

We set the initial guess σ 0(x, z) = 1 in Ω.

Step 5

We solve the problem (1) with σ = σ 0.

Step 6

We solve with the boundary condition where ∇ 2 and ∇ are two-dimensional Laplacian and gradient operators on the xy-plane, respectively, Ω a two-dimensional slice of Ω, which is perpendicular to the z-axis, and ν the two-dimensional outward normal vector to ∂Ω .

Note that in the harmonic B algorithm in [5, 16], Steps 5 and 6 are repeated to get σ from σ until it converges. In order to ensure the convergence of the harmonic B 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 instead of (1). In addition, in Step 6, we adopt the boundary condition to replace (10). Note that three-dimensional Laplacian of B may not be invalid due to the circumstance of the experiments. Hence, CoReHA supports for both 2D and 3D Laplacian of B .

3. CoReHA 1.0: For Multilateral Studies

CoReHA implements Step 2 through 6 of the harmonic B 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 B 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 B data is continuous in the xy-plane given by phase unwrapping, we do not guarantee the continuity along the z-direction. In order to verify the continuity of B along the z-direction, CoReHA supports a verification tool. Hence, if the verification fails, one should drop (∂2/∂z 2)B in ∇2 B 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 𝒫 1 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 I should be reduced less than a few milliamperes. The SNR of B 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 B data if the matrix 𝔸[σ](x) 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 B data as shown in Figure 2(b). We can see that the ramp structure of B 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].

Figure 2

Explanation of ramp structure. (a) Conductivity distribution. (b) Corresponding B data generated by the current injection along vertical direction. (c) 1D profiles along the line y = 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 B data. Locations of conductivity change match the ramp structure which is characterised by changes of slopes of the B data.

For the sake of completeness, we summarise the results in [9, 24]. Let us consider the structure tensor U defined by where To reduce the noise, we solve the following nonlinear PDE cooperating with the structure tensor U: with where R is a two-dimensional region for the slice of B data, u is the (i, j) element of U, Λ and λ are maximum and minimum eigenvalues of U, respectively, v Λ and v are corresponding normalized eigenvectors, and U ≡ G ∗U is the (element by element) convolution of U with the two-dimensional Gaussian kernel G with a standard deviation s. In the new version of CoReHA, the user only considers the parameter T 1, 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 T 2. The noise level of B affects only T 1. Consequently, we fix T 2 = 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 B since the diffusion process is on the xy-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, B 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 ∇2 B = 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.

Figure 3

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

4.3. Reconstruction: Local Harmonic B Algorithm

We support a modified local harmonic B algorithm [11] for the successful reconstruction by excluding the defected region directly. Let D 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 D instead of the whole domain Ω. In order to solve (9) with (10), we use a 𝒫 1 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 B algorithm.

Figure 4

Local harmonic B algorithm: (a) MR image, (b) selected interested region, (c) conductivity reconstruction on the whole domain, and (d) conductivity reconstruction using the local harmonic B 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 k-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, B , 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 B data in the in vivo experiment. Also, we provided a user interface for the local harmonic B 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.