Magnetic resonance current density imaging using one component of magnetic flux density

Abstract

Current density distribution generated inside a volume conductor by externally applied currents can be calculated by using spatial distribution of its magnetic flux density, . The imaging modality used to reconstruct images of the current density distribution is known as magnetic resonance current density imaging (MRCDI). In MRCDI, spatial distribution of the current‐induced magnetic flux density is measured on a magnetic resonance imaging (MRI) platform. Calculation of current density distribution from magnetic flux density measurements requires all three components of . As only the component of which is in parallel with the main field of the MRI system can be measured, rotation of the subject is necessary to be able to acquire all three components of for a given current injection pattern. Rotating the subject inside the magnet is not trivial even for small‐sized objects and remains as a strong limitation to clinical applicability of the technique. In this study, a novel MRCDI reconstruction algorithm using only one component of is proposed to eliminate the need for subject rotation. Reconstruction performance of the proposed algorithm is evaluated on numerical and experimental models. Images obtained by the proposed and the conventional MRCDI algorithms, which utilizes three components of , are compared.

Keywords:

• magnetic resonance imaging
• current density imaging
• electrical impedance tomography

1. Introduction

Joy et al. 1 have developed magnetic resonance current density imaging (MRCDI) to image distribution of electrical current density, , generated by externally applied currents to volume conductors containing magnetic resonance (MR) active nuclei based on measurements of magnetic flux density distribution. Only the component of in the same direction with the static field of an MR system, B_z, can be measured by means of MR imaging (MRI) techniques. Each orthogonal component can be measured by rotating the subject to align its corresponding axis with the main imaging field 2. Magnetic resonance electrical impedance tomography (MREIT) is another imaging modality which is used to image conductivity distribution based on magnetic flux density or current density distribution 3,4. Similar measurement techniques are adopted both in MRCDI and in MREIT, which need subject rotation. Rotating the subject inside the MRI magnet is not trivial even for small‐size objects. Rotation may result in misalignment of successive measurements, hence, may cause distortion and artefact in the reconstructed image and remains as a strong limitation to clinical applicability of the technique.

After Joy et al. 1 conducted the first MRCDI experiments, several investigators have studied MRCDI problem to calculate current density by measuring all three components of magnetic flux density 1,2,5. Others have studied imaging current density distribution by measuring only one component of magnetic flux density to eliminate subject rotation 6–11. Oh et al. 6 have used a special geometry in their study, in which current has only one component. Hence, their study is applicable only to problems with specific geometries. Seo et al. 9 have calculated current density distribution by assuming an initial conductivity distribution and solving two Neumann boundary volume problems. They have achieved good results in low‐noise simulations. However, as the measured B_z is differentiated twice, their algorithm is more prone to noise. Therefore, their algorithm is not satisfactory with noisy measurements. Ider et al. 11 have reconstructed B_z outside the object iteratively and then calculated the current density distribution using a Fourier transform (FT)‐based algorithm. In their algorithm, quality of the reconstructed current density is highly dependent on the reconstructed B_z outside the object. Therefore, errors in the reconstructed B_z cause additional errors in the images of reconstructed current density. Pyo et al. 10 and Park et al. 12 have shown that one component of magnetic flux density is not sufficient to recover current density distribution uniquely in a volume. However, Park et al. 12 have defined projected current density, which is the recoverable component of current density from measured B_z.

In this study, a novel approach to reconstruct current density in a two‐dimensional (2D) field of view (FOV), using only the component of orthogonal to the imaging plane, is proposed. In the proposed algorithm, measured magnetic flux density inside the conductor region is sufficient. Moreover, it can be applied to problems with all types of geometry, and the reconstruction time is fast. Performance of the proposed algorithm is evaluated on both measured and simulated data and compared with the conventional MRCDI algorithm, which utilizes three components of . The need for subject rotation is eliminated by introducing the proposed algorithm for MRCDI reconstruction.

2. MRCDI reconstruction using one component of B

2.1. MRCDI reconstruction algorithm using only B_z

In this study, a novel current density reconstruction algorithm using only B_z (i.e. the z‐component of ) is proposed. The Biot–Savart integral is the foundation of the proposed algorithm. The magnetic flux density due to a current element can be written as (1)

The primed variables are used to indicate the parameters corresponding to source points. Assume that for each element, current is localized at the centre of the corresponding element. Then, the current element can be written in terms of current density and the area of corresponding element as (2)

Also, the vector from source to field points, , can be written as (3)

Assuming that the current flows in the x–y plane in this 2D problem, magnetic flux density is generated only in the z‐direction, which can be written as (4)

Taking the integral of both sides in the above equation gives the desired magnetic flux density distribution. Effect of each current element on itself is neglected to overcome the singularity problem in evaluating the integral. Now, assume that 2D conductor region, Ω is divided into N finite elements. Hence, b_z can be written in matrix form as (5) where, b_z, j_x and j_y are N × 1 vectors and C_x, C_y are N × N matrices. C_x and C_y matrices only depend on . Hence, these matrices are constant for a given subject geometry.

If the current is injected to the 2D conductor region, Ω, through the electrodes placed on each side, the current density distribution can be expressed as (6) where g(x, y) is the injected current and it satisfies (7) If is assumed to be solenoidal, a differential equation relating to J_x and J_y can be obtained. However, this assumption is not valid on boundaries because of current injection. Ider et al. 11 propose use of difference currents as (8)

Here, and are the currents for uniform conductivity distribution. In order to obtain these currents, a simulation model, which has the same geometry as the conductor region but with uniform conductivity, is prepared. Then, finite element method (FEM) is used to obtain and . Note that the distribution of is independent from the chosen uniform conductivity value. The divergences of total current density distribution and uniform current density distribution are the same. Hence, the difference current density distribution is solenoidal and a relation between and can be obtained as (9)

The derivatives in Equation (9) are approximated by using central difference method as it yields a more accurate approximation than forward and backward difference methods. The discretization of Equation (9) shows that each element of can be expressed in terms of the elements of . Hence, an N × N matrix relating to the elements of and is formed as (10)

The elements of this A matrix is either zero or one depending on which elements of are related to the elements of , the subject geometry and number of finite elements, N. Equation (5) can be written for difference currents as (11)

Using Equations (10) and (11) together, the relation between and is obtained as (12) (13) is the difference magnetic flux density vector, which is calculated by subtracting from the measured magnetic flux density vector, b_z. is generated by using uniform currents, which are obtained from the simulation model. is calculated by using Equation (13) and then solution of is trivial. Note that, C_t is a non‐singular and well‐conditioned matrix. Lastly, j_x and j_y are obtained by adding difference currents and uniform currents. To summarize, the proposed algorithm has the following steps:

Step 1

determining the subject geometry and number of finite elements, N,

Step 2

calculating C_x, C_y and A matrices,

Step 3

calculating C_t,

Step 4

computing using Equation (13), and then compute .

2.2. Performance evaluation of the proposed MRCDI reconstruction algorithm

2.2.1. The experimental setup and the phantom

Magnetic flux density measurements are performed on 0.15T Middle East Technical University (METU) MRI system. A standard spin‐echo pulse sequence with an addition of a bipolar current pulse is used. Current is applied to the phantom with a current source synchronized with the RF pulse.

The experimental phantom is prepared in order to test the performance of the algorithm with measured data. Experimental phantom is made of Plexiglas® and its dimensions are 9 × 9 × 9 cm³. The 2D geometry is obtained with additional Plexiglas® walls which force the applied current to flow in a volume of 9 × 9 × 2 cm³. Electrodes are placed in the middle of each side and their dimensions are 2 × 2 cm². The phantom elements are prepared by using solidifying materials, Agar‐Agar, TX‐150 and TX‐151. Also, NaCl is used to adjust conductivity values of elements, and CuSO₄ is used to fix T1 relaxation time. In this experimental study, conductivity values of background and circle element are 0.2 and 2 S m⁻¹, respectively. Square element is prepared as a ‘pure insulator’. The geometry of the phantom elements and their dimensions are as shown in Figure 1. Composition of materials used to prepare the phantom elements is given in Table 1. A photograph of the phantom and its MR magnitude image are given in Figure 2.

Figure 1. The geometry of phantom elements.

Figure 2. (a) Top view of the phantom and (b) cross‐sectional MR magnitude image of the phantom at the plane of the electrodes.

Table 1. Composition of phantom elements.

	Square object (∼0 S m^−1)	Circular object (2 S m^−1)	Background (0.2 S m^−1)
Distilled water	100 mL	100 mL	100 mL
Agar–Agar	1 g	1 g	–
TX‐150	–	–	0.2 g
TX‐151	1 g	1 g	2 g
NaCl	–	1.75 g	–
CuSO4	0.1 g	0.1 g	0.1 g

Two current injection patterns are used. The first current injection pattern is vertical current injection pattern, where current is injected from the upper electrode and sunk from the lower electrode. The other one is horizontal current injection pattern, where current is injected from the electrode on the left and sunk from the electrode on the right. 20 mA current is injected to the phantom in both simulations and experiments.

2.2.2. Generation of simulated data

The forward problem of MRCDI is solved using COMSOL® Multiphysics FEM solver. The FEM model is implemented for the geometry given in Figure 2 with a thickness of d = 2 cm. Once the current density distribution is obtained from the FEM solution, B_z is calculated by using Equation (5).

2.2.3. Noise model for computer simulations

A noise model is prepared in order to evaluate the performance of the algorithm on noisy data. In this study, random Gaussian noise model is used for noise simulations 13. This noise model only depends on signal‐to‐noise ratio (SNR) of imaging system where magnetic flux density is measured. In this study, SNR 30 and SNR 13 levels are used in simulations because SNR 30 corresponds to an MRI system with 2 T magnet and an SNR value of 0.15 T METU‐EE MRI system is around 13 14.

To simulate noisy data, current density distribution is obtained using FEM. Then, magnetic flux density is calculated by the Biot–Savart law from these simulated current density distributions. Random Gaussian noise is added to magnetic flux density based on the SNR level of MRI System. Finally, noisy magnetic flux density is used in the MRCDI reconstruction algorithms.

2.2.4. Error calculation for simulated data

To evaluate the performance of the algorithm quantitatively, the error in reconstructed current density distribution is calculated as (14) Here, J_r and J_c represent the real and calculated values of current density distribution, respectively.

2.2.5. Error calculation for experimental data

In order to measure all three components of , the object to be imaged should be rotated inside the MRI scanner. The phantom rotation causes artefacts in the reconstructed due to variations in the internal geometry during rotation. Hence, obtained by using three components of deviates from the true current density. In other words, there is no ground truth current density distribution for experimental data. Therefore, quantitative evaluation of the reconstruction performance of the proposed algorithm on measured data, relative to the performance of the algorithms using all component of , is not possible. However, the reconstructed can be verified by using divergence theorem, and performance of the proposed algorithm can be evaluated in terms of error in the reconstructed MREIT conductivities, utilizing J‐based MREIT algorithms on the reconstructed current density distributions.

From the divergence theorem, it is known that the integral of over each slice, which is orthogonal to the current injection pattern, should be the same and equal to the externally applied current. Hence, standard deviation of the integrals in orthogonal slices should be zero. The standard deviations of the integrals give valuable information about the reconstruction performances of the algorithms.

The performances of the MRCDI reconstruction algorithms can be evaluated in terms of error in the reconstructed MREIT conductivities. For this purpose, J‐substitution algorithm, proposed by Kwon et al. 15, is adopted as it provides higher quality conductivity images compared with the other current density‐based reconstruction algorithms 16. This algorithm is an iterative algorithm, which uses only the magnitude of current density inside the FOV. The algorithm tries to minimize a cost function, which is defined as the difference between calculated and measured current density distributions. Error for the reconstructed conductivity images are calculated as (15) where N, σ_ir and σ_ic are number of pixels, real and calculated conductivity values, respectively.

3. Results

3.1. Simulation results

Figure 3 shows the z‐component of magnetic flux density, B_z, obtained from vertical and horizontal current injection patterns.

Figure 3. The z‐component of magnetic flux density, B_z (in Tesla): (a) for horizontal current injection pattern and (b) for vertical current injection pattern.

Reconstructed J_x and J_y for different noise levels for the horizontal current injection pattern are shown in Figure 4.

Figure 4. Reconstructed J_x and J_y for different noise levels (horizontal 20 mA current injection): (a) original J_x (A m⁻²) distribution, (b) original J_y (A m⁻²) distribution, (c) reconstructed J_x (A m⁻²) using the proposed algorithm for noise‐free case, (d) reconstructed J_y (A m⁻²) using the proposed algorithm for noise‐free case, (e) reconstructed J_x (A m⁻²) using the proposed algorithm for the SNR 30 case, (f) reconstructed J_y (Am⁻²) using the proposed algorithm for the SNR 30 case, (g) reconstructed J_x (A m⁻²) using the proposed algorithm for the SNR 13 case and (h) reconstructed J_y (A m⁻²) using the proposed algorithm for the SNR 13 case.

When reconstructed current density images are evaluated perceptually, it is seen that reconstruction performance of the proposed algorithm is good. The effect of noise can be distinguished in the SNR 30 and SNR 13 cases. However, the current density is reconstructed without major artefacts in the presence of noise. The quantitative evaluation of the reconstruction performance of the proposed algorithm is listed in Table 2. Quantitative evaluation of the reconstruction performance of the proposed algorithm shows that the proposed algorithm reconstructs with error less than 5% using the simulated data having the worst SNR. The reconstruction performance of the proposed algorithm for the component of current density orthogonal to the current injection pattern is poorer. This is expected since the component of current density in the direction of injected current is stronger.

Table 2. Errors in the reconstructed current density for simulated data.

	Horizontal current injection			Vertical current injection
	Jx	Jy	J	Jx	Jy	J
Noise‐free (%)	0.34	1.34	0.69	1.25	0.89	0.86
SNR 30 (%)	1.11	2.14	1.18	2.49	1.49	1.40
SNR 13 (%)	2.48	3.80	2.19	4.98	2.85	2.60

3.2. Experimental results

The proposed algorithm is also tested utilizing measurements performed using the experimental phantom described in Figures 6 and 7. 20 mA current is injected to the phantom as the horizontal and the vertical current injection patterns. The z‐component of total measured magnetic flux density is shown in Figure 5. In experimental studies, uniform magnetic flux density is calculated from the simulation model and the difference magnetic flux density is obtained by subtracting uniform magnetic flux density from the measured magnetic flux density.

Figure 5. Magnetic flux density in the z‐direction, B_z (in Tesla): (a) for the horizontal current injection pattern and (b) for the vertical current injection pattern.

Figure 6. Results for the horizontal current injection pattern: (a) reconstructed J_x (A m⁻²) using three components of magnetic flux density, (b) reconstructed J_y (A m⁻²) using three components of magnetic flux density, (c) reconstructed J_x (A m⁻²) using the proposed algorithm and (d) reconstructed J_y (A m⁻²) using the proposed algorithm.

Figure 7. Results for the vertical current injection pattern: (a) reconstructed J_x (A m⁻²) using three components of magnetic flux density, (b) reconstructed J_y (A m⁻²) using three components of magnetic flux density, (c) reconstructed J_x (A m⁻²) using the proposed algorithm and (d) reconstructed J_y (A m⁻²) using the proposed algorithm.

The conventional MRCDI algorithm needs all three components of to calculate . Therefore, other components of the magnetic flux density, B_x and B_y, were also measured in the experiments, which makes it possible to compare the proposed algorithm with the conventional MRCDI algorithm.

The results of the experimental model for the horizontal and the vertical current injection patterns are shown in Figure 6 and Figure 7, respectively. Note that the direction of positive y‐axis is taken as negative in J_y for vertical current injection to make the images consistent with J_x images for the horizontal current injection pattern.

As mentioned in Section 2.2.4, the reconstructed current density images are verified by the divergence theorem. In the conventional MRCDI algorithm, the standard deviations are found as 11 × 10⁻⁴ and 8 × 10⁻⁴ for horizontal and vertical current injection patterns, respectively. In the proposed algorithm, the standard deviations are found as 8.77 × 10⁻⁴ and 5 × 10⁻⁴ for horizontal and vertical current injection patterns, respectively.

J‐substitution algorithm is used to reconstruct MREIT images from the reconstructed current density with five iterations. The reconstructed conductivity images are shown in Figure 8. Evaluating the performances of the current density reconstruction algorithms in terms of error in the reconstructed conductivity images is not possible because of the reconstruction artefacts at the corners of the FOV. These artefacts occur as a result of the low‐current density at the corners of the phantom. The reconstructed conductivity values of these corners dominate the conductivity of the remaining regions. Therefore, the reconstructed conductivity images can only be evaluated perceptually. The J‐substitution algorithm seems to reconstruct the better conductivity image by using obtained from the proposed algorithm. The circular object can be clearly distinguished and it is reconstructed with 26.3% error. The circular object can be distinguished in Figure 8(b) as well, it is reconstructed with 43.5% error.

Figure 8. Reconstructed conductivity images with J‐substitution algorithm after five iterations, (a) original conductivity distribution, (b) reconstructed conductivity distribution using current density obtained from the conventional MRCDI algorithm and (c) reconstructed conductivity distribution using current density obtained from the proposed algorithm.

4. Discussion and conclusion

A novel 2D MRCDI reconstruction algorithm has been proposed in this study. The performance of the proposed algorithm is evaluated with simulation and experimental models. It is also compared with the conventional MRCDI algorithm using three components of . The simulation results show that the proposed algorithm reconstructs with error less than 5% in the worst SNR case. In the experiment, the proposed algorithm gives satisfactory results when it is compared with the conventional MRCDI algorithm. It is important to note that using one component of magnetic flux density significantly reduces the time spent during data acquisition. This time can be spent to increase the number of measurements hence improve the SNR.

The reconstruction time is an important criterion in the performance evaluation of an MRCDI algorithm. Since the matrices used in the proposed algorithm only depend on the geometry of the subject and number of finite elements, these matrices can be calculated and stored for a subject with known geometry. Therefore, performing only the matrix inversion gives the current density. The elapsed time for matrix inversion is in the order of seconds for usual MR image resolutions.

In the proposed algorithm, magnetic flux density for a uniform conductivity distribution, , is calculated from a simulation model. Because, subtraction of from is required, misaligned electrode locations may introduce undesirable deformations in the reconstructed image in experimental studies. This is a disadvantage of the algorithm. By using low‐pass filters, the deformations in the reconstructed current density distribution can be suppressed at the expense of reduced spatial resolution in the reconstructed current density images.

It is expected that the quality of the reconstructed current density images obtained by using the proposed algorithm should be higher than those algorithms using FT 11. This is due to the fact that the algorithms using FT require the magnetic flux density outside the conductor region and must calculate it by processing the measured B_z data before solving the current density. Errors in calculated magnetic flux density outside the conductor region also cause errors in the reconstructed current density. As measured magnetic flux density inside the conductor region is sufficient for the proposed algorithm, error in reconstructed current density images is expected to be smaller. In Ersoz 17, the performance of the proposed algorithm was compared with the Iterative‐FT MRCDI algorithm, which is proposed by Ider et al. 11. Both algorithms were evaluated with the model in Figures 1 and 2. The results show that the error in the reconstructed current density reaches up to 4.98% for the proposed algorithm and 20.00% for the iterative FT‐MRCDI algorithm in the SNR 13 case. Ider et al. 11 have also calculated the error in the reconstructed current density as 11.8% for no noise case in their simulation model. The higher error in their algorithm is possibly due to processing of measured B_z data.

2D current density distributions are reconstructed in this study. However, most of the practical applications involve 3D distributions. One component of magnetic flux density is not sufficient to obtain current density in 3D distributions, uniquely 10. Application of this algorithm to image 3D current flows may introduce undesirable deformations in the reconstructed image due to non‐zero currents flowing in the z‐direction and it should be studied further.

Acknowledgements

This study is part of Ali Ersöz′s MSc thesis, and B.M. Eyüboğlu is the thesis supervisor. The authors would like to thank Evren Değirmenci for his help in performing experiments in the METU‐MRI system and processing of the experimental data. This study is partially supported by the Turkish Scientific and Technological Research Council (TUBITAK) under research grant 107E141.