Conductivity image enhancement in MREIT using adaptively weighted spatial averaging filter

Background In magnetic resonance electrical impedance tomography (MREIT), we reconstruct conductivity images using magnetic flux density data induced by externally injected currents. Since we extract magnetic flux density data from acquired MR phase images, the amount of measurement noise increases in regions of weak MR signals. Especially for local regions of MR signal void, there may occur excessive amounts of noise to deteriorate the quality of reconstructed conductivity images. In this paper, we propose a new conductivity image enhancement method as a postprocessing technique to improve the image quality. Methods Within a magnetic flux density image, the amount of noise varies depending on the position-dependent MR signal intensity. Using the MR magnitude image which is always available in MREIT, we estimate noise levels of measured magnetic flux density data in local regions. Based on the noise estimates, we adjust the window size and weights of a spatial averaging filter, which is applied to reconstructed conductivity images. Without relying on a partial differential equation, the new method is fast and can be easily implemented. Results Applying the novel conductivity image enhancement method to experimental data, we could improve the image quality to better distinguish local regions with different conductivity contrasts. From phantom experiments, the estimated conductivity values had 80% less variations inside regions of homogeneous objects. Reconstructed conductivity images from upper and lower abdominal regions of animals showed much less artifacts in local regions of weak MR signals. Conclusion We developed the fast and simple method to enhance the conductivity image quality by adaptively adjusting the weights and window size of the spatial averaging filter using MR magnitude images. Since the new method is implemented as a postprocessing step, we suggest adopting it without or with other preprocessing methods for application studies where conductivity contrast is of primary concern.

http://www.biomedical-engineering-online.com/content/13/1/87 Background Electrical conductivity is a passive material property of a biological tissue or organ providing diagnostic information of its physiological function and pathological state [1,2]. In magnetic resonance electrical impedance tomography (MREIT), we aim to visualize the internal conductivity distribution of the human body by injecting electrical current and measuring induced magnetic flux density data using an MRI scanner [3,4]. The MREIT technique uses the z-component B z of the magnetic flux density B = (B x , B y , B z ) to recover conductivity and/or current density images [5][6][7][8][9][10].
In experimental MREIT studies of animal and human subjects [11][12][13][14], the quality of reconstructed conductivity images highly depends on the noise level in measured B z data and the injection current amplitude. For a given current amplitude, we can improve the image quality by reducing the noise level. If we reduce the current amplitude to avoid adverse effects such as electrical stimulations of nerve and muscle, the range of B z decreases proportionally and the measured B z data become more vulnerable to noise. In order not to deteriorate the conductivity image quality, we should reduce the noise level as well.
There have been numerous studies to minimize the noise level in measured B z data by optimizing the data collection method including pulse sequences, RF coils, shimming, averaging, and so on. For example, the injected current nonlinear encoding (ICNE) method was introduced to reduce the noise level in B z data by extending the current injection time without increasing the scan time [15]. Various multi-echo and multicoil techniques have been developed together with optimization methods to combine multiple signals from multiple coils and echoes [16].
Once we acquire k-space data from a designed MREIT experiment, we should carefully process the data to suppress the noise and enhance the image quality. This requires proper understanding of the noise characteristics. The noise standard deviation of B z is inversely proportional to the signal-to-noise ratio (SNR) of the MR magnitude image and the current injection time when the magnitude image SNR is higher than about 3 [17,18]. Since the magnitude image SNR varies over the image, the noise level of B z also varies. If the magnitude image SNR is much lower, for example, less than 3, then we should expect an excessive amount of noise in B z . This may happen when we scan an animal or human subject since there exist local regions of weak MR signal such as the lungs, gas-filled organs, and outer layers of bones. In those regions with very low magnitude image SNRs, there occur excessive amounts of noise in measured B z data. It is, therefore, important to properly suppress the noise effects in preprocessing or postprocessing steps.
There are several denoising methods to process the acquired B z data [19][20][21][22][23]. These methods use either a PDE-based approach or a localized data processing approach and require a certain degree of manual adjustment of parameters. All of these methods are applied to measured B z data before conductivity image reconstructions as a preprocessing step.
In this paper, we present a new method to suppress the noise effects in reconstructed conductivity images as a postprocessing step. We suggest incorporating a prior information from the MR magnitude image which is always available in MREIT. For current density imaging using measured magnetic flux density data, Joy et al. suppressed noise effects by setting some appropriate threshold in the MR magnitude image SNR [24]. A similar approach based on the morphological enabled dipole inversion (MEDI) has been http://www.biomedical-engineering-online.com/content/13/1/87 proposed in quantitative susceptibility map (QSM) to overcome the ill-posedness in its inverse problem [25].
The proposed method enhances the quality of reconstructed conductivity images by using an adaptively weighted spatial averaging filter. To determine the weights, we will define a distance function with respect to the noise standard deviation in measured B z data. Since the noise is inversely proportional to the MR magnitude image, we will incorporate the magnitude image in the spatial filtering process. After describing the details of the proposed method, we will show how it performs with experimental data from a conductivity phantom and also animal subjects.

Noise in magnetic flux density data
We inject current into an imaging object, of which timing is synchronized with a chosen MR pulse sequence. The usual current injection period is from the end of the exciting RF pulse to the beginning of the readout gradient. One may choose different current injection times and patterns depending on designed pulse sequences. For example, the ICNE-multi-echo pulse sequence extends the duration of current injection until the end of multiple read-out gradients.
The externally injected current produces an internal magnetic flux density distribution and its z-component B z in a voxel results in an extra phase in the MR phase image. To remove systematic phase artifacts, we sequentially inject positive I + and negative I − current pulses to obtain the following k-space data: M(x, y)e iδ(x,y) e ±iγ B z (x,y)T c e i2π(k x x+k y y) dxdy (1) where M is the MR magnitude image, δ is any systematic phase artifact, γ = 26.75 × 10 7 rad/T · s is the gyromagnetic ratio of the proton, and is a field-of-view (FOV). Here, the superscript of S ± denotes a brief notation for S + and S − . We extract the magnetic flux density B z by where ± = Me iδ e ±iγ B z T c . The noise standard deviation of the measured magnetic flux density B z is inversely proportional to the current injection time T c and the SNR of the MR magnitude image ϒ as for ϒ > 2.8 [17,18]. If we reduce the current amplitude I ± , the range of B z decreases proportionally. To reduce the noise standard deviation, we have to increase the current injection time T c and the SNR ϒ simultaneously. However, this is not possible mainly due to the T 2 or T * 2 decay of the MR signal.

Conductivity image reconstruction
The externally injected current I induces distributions of voltage u, current density J = (J x , J y , J z ), and magnetic flux density B = (B x , B y , B z ) inside the imaging object with its boundary ∂ . The voltage satisfies the following partial differential equation: http://www.biomedical-engineering-online.com/content/13/1/87 where r = (x, y, z), n is the outward unit normal vector on ∂ , and g denotes the Neumann boundary data subject to the injection current. From the Biot-Savart law, the magnetic flux density B z is related with the current density J = −σ ∇u as where μ 0 is the magnetic permeability of the free space. There exists numerous image reconstruction algorithms to visualize the conductivity σ or the current density J in the imaging object from measured B z data [5,9,16,26]. In this paper, we adopted the transversal J-substitution algorithm [27] since it differentiates the noisy B z data once and does not propagate the noise effects from one region to another. We first assume the imaging object with a homogeneous conductivity distribution σ H . Solving (4) with σ H in place of σ , we can compute the voltage, current density, and magnetic flux density, which are denoted as u H , J H , and B H z , respectively. In MREIT experiments, we measure two B z data subject to two orthogonal injection currents. Denoting them as B z,1 and B z,2 , the transversal J-substitution algorithm produces an image of the conductivity σ as where∇ ⊥ f := ∂f ∂y , − ∂f ∂x for a given scalar function f and ·, · denotes the scalar inner product. We used the computed values of two magnetic flux densities B H z,1 and B H z,2 corresponding to two computed voltages u H 1 and u H 2 , respectively, for the homogeneous case of σ H . Note that the noisy data of B z,1 and B z,2 may deteriorate the quality of the reconstructed image of σ . More details about the image reconstruction algorithm are described in [27].

Adaptively weighted spatial averaging
To determine a neighborhood of a pixel, we define the following distance function: where M(r) is the MR magnitude image, B r (η(r)) is a neighborhood of the point r with a radius η(r), and h(r) is a function of the noise level in measured B z data. For each pixel s in the neighborhood of the pixel r, that is, s ∈ B r (η(r)), we define the following weighting factor w r (s): where ζ r := s e −D(r,s) is a normalization constant ensuring that s w r (s) = 1. We perform the adaptive spatial averaging of the reconstructed conductivity values as the following weighted sum: ( 9 ) http://www.biomedical-engineering-online.com/content/13/1/87 For the pixel at r, the distance D(r, s) measures the similarity of M(r) and M(s) in the surrounding region of r. Since the noise standard deviation of B z is inversely proportional to the magnitude image SNR, the weights are in effect adjusted by the noise level in measured B z data. In (7), two parameters of η(r) and h(r) regulate the extent of the spatial averaging filter. We now design the denominator h(r) and the radius η(r) in (7) as We can choose the proportionality parameters of h(r) and η(r) in (10) depending on the quality of the reconstructed conductivity image to be filtered.
The designed distance function D(r, s) preserves the conductivity value where the noise level of B z is low. To determine the characteristics of the weighting factor w r (s), we decompose the weighted conductivity σ w (r) to the true conductivity σ t (r) and the noise term: where σ t (s) and N r (s) denote the noiseless true conductivity and the noise term at s, respectively. The weighted conductivity σ w (r) and the true conductivity σ t (r) satisfy the following relation: where and E 2 (r) = 1 ζ r s∈B r (η(r)) e −D(r,s) N r (s) .
For a small amount of noise, both h(r) and η(r) are small and, therefore, the relation (12) implies that both E 1 (r) and E 2 (r) are small. This means that the method does no harm to the conductivity image when the noise level is low.
For a large amount of noise, the first error term E 1 (r) becomes small if there was no significant variations of the conductivity values within the neighboring region B r (η(r)) of the pixel at r. This is true only when the pixels with similar magnitude image values also have similar conductivity values. We will discuss implications of this restriction later. The second error term E 2 (r) is considerably reduced by the spatial averaging of the random noise.

Imaging experiments Phantom imaging
To evaluate the performance of the proposed method, we scanned a cylindrical phantom with four carbon-hydrogel electrodes (HUREV Co. Ltd, Korea). The phantom was filled with 0.4 S/m saline and included three cylindrical objects. The first (D 1 ) was a TX-151 http://www.biomedical-engineering-online.com/content/13/1/87 object (1.5 S/m), the second (D 2 ) was an agar object (0.1 S/m), and the third (D 3 ) was a TX-151 object (1.5 S/m) wrapped with an agar object (0.3 S/m) as shown in Figure 1(a). We placed the phantom inside the bore of our 3 T MRI scanner (Magnum 3, Medinus Co. Ltd., Korea). Using a custom-designed MREIT current source [28], we injected 10 mA currents in the horizontally and vertically directions. The parameters of the spinecho-based imaging sequence were as follows: repetition time T R = 1200 ms, echo time T E = 15 ms, number of echoes N E = 5, field of view (FOV) was 240 × 240 mm 2 , number of excitations (NEX) was 1, imaging matrix was 128 × 128, and total imaging time was 10.24 min. Figure 1(b) shows the injected current pulses synchronized with RF pulses. We combined the acquired multiple echoes to optimize the SNR in measured B z data [16]. Figure 1(c) shows the measured B z image subject to the vertical current injection.

Animal imaging
The animals were laboratory beagles (2-3 years old, weighing 8-15 kg) without having any known disease. To prevent dribbling, we injected 0.1 mg/kg of atrophine sulfate. Ten minutes later, we anesthetized the dog with intramuscular injection of 0.2 ml/kg Tiletamine and Zolazepam (Zoletil 50, Virbac, France). Twenty minutes later, we sacrificed it with an intravenous injection of 80 mg/kg KCL (Entobar, Hanrim Pharmacy, Korea). After clipping the hair, we attached four carbon-hydrogel electrodes (HUREV Co. Ltd., Korea) around the imaging area. The size of each electrode was 80 × 80 × 6 mm 3 . The procedure was approved by the Institutional Animal Care and Use Committee (IACUC). We injected currents in two mutually orthogonal directions between two pairs of electrodes facing each other. The injection current amplitude ranged from 5 to 10 mA. We adopted the ICNE pulse sequence. The imaging parameters were as follows: repetition time T R = 1200 ms, echo time T E = 30 ms, FOV was 280 × 280 mm 2 , slice thickness was 4 mm, number of slice was 8, NEX was 6, imaging matrix was 128 × 128, and total imaging time was 60 min. Figure 2 shows two data sets we used in the paper. Figure 3(a) shows the MR magnitude image at the middle slice of the phantom. We segmented three regions of D 1 (TX-151), D 2 (agar), and D 3 (TX-151 wrapped with agar). Depending on the noise level in measured B z data, we determined the denominator h in (7) to compute the distance function D(r, s). Figure 3(b) shows the denominator h in (7). the short T 2 relaxation time of the objects, the reconstructed conductivity values of the objects show strong noise defects. Figure 3(d) shows the filtered conductivity image using the proposed method. In the noisy regions of D i , i = 1, 2, 3, the values of the denominator h were relatively large to produce smaller distance values in (7). Since the values of η were also large in D i , more pixels in the neighboring region were included in the spatial averaging with relatively larger weights.

Results
We may quantify the effects of the filtering method by computing the amount of conductivity changes in a homogeneous region. We defined the following variance function to estimate how the reconstructed conductivity values varied in each region of interest (ROI): where |D i | denotes the volume of the ROI D i . Table 1 shows that the proposed method significantly improved the image quality. Figure 4(a) is the image of the denominator h in (7) for the upper abdominal region in Figure 2(a). The areas marked by the arrows have large values of h since their MR magnitude image SNRs were low. The average magnitude intensity in the marked areas was 3.78, where the intensity of the entire image ranged from 0 to 30. Since the noise level of B z is inversely proportional to the MR magnitude intensity, the recovered conductivity values of the marked regions included a relatively large amount of noise. Figure 4(b) shows the reconstructed conductivity image without using any filtering method and we can observe severe noise effects in those regions.  To compare the proposed method with a conventional denoising technique, we used the following reaction-diffusion iteration step [29]: where f is the reconstructed conductivity image shown in Figure 4(b). Figure 4(c) shows the denoised conductivity image using the iteration method (16). We can see that the iteration method blurred the entire image. Figure 4(d) is the conductivity image after applying the adaptively weighted spatial averaging filter. We can see that the proposed method enhanced the conductivity image in the marked regions without affecting other regions.
For the blurred image in Figure 4(c), we used the following parameters: iteration number = 200, α = 0.1, and fidelity term β = 0.01. These parameters were chosen to make the variance of the conductivity values in the marked regions of the image in Figure 4(c) to be equal to that of the image in 4(d). Figure 5 shows similar results for the lower abdominal region shown in Figure 2(c). The average magnitude intensity of the marked regions was 1.97 for this case and the measured B z data were noisier than the case in Figure 4. Comparing two filtered conductivity images in 5(c) and (d), we can see that the method proposed in this paper is clearly advantageous.

Discussion and conclusion
For animal or human subjects, their internal structures are heterogenous and there often exist local regions of weak MR signals. This makes the noise level in measured magnetic flux density data vary significantly for different pixels. Since we use the data to reconstruct conductivity images, noise effects are conveyed to the conductivity images with spatially varying image quality. Conventional filtering methods without considering this property, therefore, unnecessarily blur the entire conductivity image. In this paper, we proposed a fast and simple method to enhance the conductivity image quality by utilizing the MR magnitude image. Noting that the noise level in measured magnetic flux density data is inversely proportional to the pixel value of the MR magnitude image, we could adaptively adjust the weights and window size of the spatial averaging filter. In designing the filter, we used two parameters of η(r) and h(r) in (7) to regulate the extent of the spatial averaging filter. Since we set their values to be proportional to the noise level in measured B z data, the proposed method adaptively changes the extent of the spatial averaging depending on the noise level at each pixel. It is important to properly choose the proportionality constants for η(r) and h(r) in (10). Though we heuristically chose the constants in this paper, we plan to rigorously investigate their effects on the filtered image and develop an automatic method to determine them.
Since the spatial averaging is performed on neighboring pixels with similar MR magnitude image values, it may remove any useful conductivity contrast among those pixels. Using the proposed method, we can apply the spatial averaging to local regions of severe noise where we can not trust reconstructed conductivity values. We can, therefore, suppress noise for the price of reduced image resolution within local regions with large amounts of noise. The proposed method negligibly influences the conductivity image in other regions with enough SNR.
For the phantom image in Figure 3, we found that the spatial averaging recovered the correct conductivity values in the noisy regions since there was no conductivity contrast within each region. For the animal images in Figures 4 and 5, it was difficult to assert that the recovered conductivity values of the noisy regions were correct since there could have been some conductivity contrast among those pixels with similar MR magnitude image values. In application studies, we may analyze both unfiltered and filtered conductivity images if quantitative conductivity values of noisy regions are of primary concern.
While the existing denoising methods in MREIT are applied to measured B z data before conductivity image reconstructions [19][20][21][22][23], the new method proposed in this paper is a postprocessing method, which can be applied to reconstructed conductivity images. Without requiring any image segmentation, one can easily implement the method without or with chosen preprocessing steps. We suggest adopting the proposed method in future experimental MREIT studies where conductivity contrast is of primary concern to extract diagnostic information.