Noise distribution and denoising of current density images

Abstract. Current density imaging (CDI) is a magnetic resonance (MR) imaging technique that could be used to study current pathways inside the tissue. The current distribution is measured indirectly as phase changes. The inherent noise in the MR imaging technique degrades the accuracy of phase measurements leading to imprecise current variations. The outcome can be affected significantly, especially at a low signal-to-noise ratio (SNR). We have shown the residual noise distribution of the phase to be Gaussian-like and the noise in CDI images approximated as a Gaussian. This finding matches experimental results. We further investigated this finding by performing comparative analysis with denoising techniques, using two CDI datasets with two different currents (20 and 45 mA). We found that the block-matching and three-dimensional (BM3D) technique outperforms other techniques when applied on current density (J). The minimum gain in noise power by BM3D applied to J compared with the next best technique in the analysis was found to be around 2 dB per pixel. We characterize the noise profile in CDI images and provide insights on the performance of different denoising techniques when applied at two different stages of current density reconstruction.


Introduction
Low frequency current density imaging (LFCDI) is a magnetic resonance imaging (MRI)-based technique which measures current density inside a subject while a current pulse is injected into the subject. The method was developed in the late 1980s and early 1990s by Joy et al. 1,2 This method has been tested on phantoms 1,3,4 and also applied to the in vivo and ex vivo tissues in order to obtain the current density map inside the subject. [5][6][7][8][9] The method works based on the fact that the magnetic flux induced by injected current can be measured using phase changes recorded in MRI phase images. 1 Although it is a very useful method to study the current density maps, like any other method it has some limitations in terms of image quality and implementation. The images can become very noisy and unusable if acquisition parameters are not chosen appropriately, the energy of the injected current pulse is not enough, or the size of the studied subject is very small. The noise and its effect on the current density images are studied in the literature in two main works by Scott et al. 10 and Sadlier et al. 11 The former studied the effect of MRI acquisition parameters like echo time (TE), injected current pulse width (Tc), and current amplitude on the noise level in the images. The latter used the fact that the noise distribution in real and imaginary images is zero mean Gaussian. These real and imaginary images are used to calculate the phase changes generated by induced magnetic flux. Sadlier et al.
subsequently studied the effect of this Gaussian noise on the magnetic flux induced in one direction and also investigated the effect of the strength of the main magnetic field of the MRI machine on the level of measurable currents. Another work by Lee et al. 12 suggested ramp-preserving denoising for CDI images by developing a new ramp preserving method based on the Perona-Malik denoising method. 13 In this work, the probability distribution of the measured phase for each pixel inside the subject is derived and the current density distribution is approximated. Furthermore, because the distribution of noise in the real and imaginary images is Gaussian and in the current density case the derived distribution will be shown to be close to a Gaussian, three different denoising methods are applied on a cylindrical phantom. These three methods are ramp-preserving Perona-Malik (PM) denoising, 13 two-dimensional adaptive Wiener filter, 14 and the state-of-theart block-matching and three-dimensional (BM3D) filtering. 15 The PM method is chosen as a ramp-preserving method already used for CDI images, Wiener as a smoothing method, and BM3D because it is shown to have a better performance compared with the traditional denoising methods. 15 The denoising is performed at two different stages: first on the initial noisy real and imaginary images recorded by MRI machine and second on the final current density maps, and the results are compared with the desired noiseless current density maps.

Method
The experimental data used in this paper were obtained by performing a CDI experiment on a cylindrical phantom filled with 9 g∕L · NaCl and 0.64 g∕L · CuSO 4 solution with T2 ¼ 170 ms and T1 ¼ 200 ms. 16 The phantom had a diameter of 4.5 cm and length of 9 cm. Two datasets were obtained: for set 1, the current amplitude was 20 mA and for set 2 the current was 45 mA. The first dataset had six slices and the second one had five slices. The current pulses were generated by a pulse generator circuit and amplified using an HP Harrison 6824A power supply-amplifier which supports maximum peak-topeak voltage of 50 V for currents up to 1 A. The current injection system was constant current source. The current was maintained constant and was measured as a voltage across a 10 Ω resistor in series with the phantom. Imaging was conducted in a 1.5 T GE machine at Toronto General Hospital. The acquisition parameters were TE equal to 40 ms, Tc equal to 18 ms, number of excitations equal to 2, repetition time equal to 700 ms, a field of view of 15 cm, and the image size was 256 × 256 pixels. The current was injected through copper disc electrodes at the top and bottom of phantom. The current pulse sequence for LFCDI is shown in Fig. 1. 17 The axis of the cylinder was aligned with the z-direction and it was parallel to the direction of the current. The coordinate system was adopted from Scott et al. 17 in which the z-direction is perpendicular to the main magnetic field. The slices were equally spaced and were perpendicular to the z-axis. Two hollow bars of Plexiglas were positioned inside the phantom to track different orientations. Due to the symmetry of the cylinder, there was only one component of current density, which was in the z-direction. Therefore, only Bx and By were measured by having a 90-deg rotation of the cylinder around its axis. A spin-echo sequence was used for imaging and six slices were defined with equivalent distance and in the x − y plane. Figure 2 shows the empty phantom used for this experiment and the magnitude image for one slice.

Distribution of Noise
As reported in the literature, the noise in the background of real and imaginary images is zero mean Gaussian. 11 This fact was experimentally verified for the phantom experiment. The histograms of background for real and imaginary images in a slice are shown in Figs , respectively. It can be seen that the distribution of background noise is almost zero mean Gaussian and it has approximately the same distribution for real and imaginary images. The standard deviation of this noise can be extracted from the background. 11,18 Therefore, if we assume that for each pixel inside the subject X R is the random variable representing the corrupted intensity in real images and X I is the random variable representing the corrupted intensity in imaginary images, we have for one pixel (1) N R is the noise random variable for the chosen pixel in the real image and N I is the random variable representing noise in the imaginary image for the same pixel. Therefore, X R and X I are Gaussian at each pixel with nonzero means R and I at that pixel. In this case, the joint distribution of X R and X I is a bivariate Gaussian distribution given as The distribution of the angle and magnitude of X R þ iX I can be found using the Jacobian of the transformation as The symbol arg represents the argument of (X R þ iX I ). Because the current density is measured from the phase changes, we find fðθÞ by integrating over r fðθÞ Defining Using the following identities the distribution can be simplified to While erf is Gaussian error function, Eq. (12) shows the wrapped phase distribution. Unwrapping will not change the distribution, it will only shift the mean by 2kπ. The unwrapped values are proportional to the induced magnetic flux (B). Therefore, we can reason that B at each point has a distribution similar to Eq. (12) just scaled or shifted.
To calculate each component of current density, derivatives of two perpendicular components of the induced B should be used. The derivatives are calculated using derivative templates. 10 These derivative templates rely on an averaging operation in a neighborhood to approximate the derivatives, therefore, their application compromises the spatial resolution and itself acts like denoising. If a 3 × 3 template is used, the nine values of Bs in the neighborhood of each point are used to calculate the derivative at each point and two derivatives. Therefore, the final distribution is the summation of 18 scaled distributions of the form provided in Eq. (12). Using the fact that the distribution of summation is a convolution of the distributions, we can find the final distribution of the current density (J). It is also known that the convolution of Gaussians is Gaussian, 19 therefore, if the distribution in Eq. (12) is close to Gaussian, we can conclude the noise distribution in the measured current density is also Gaussian. Figures 4 and 5 show the distribution for various values of R, I, and σ. The values chosen for R and I in these figures are arbitrarily chosen in the range of observed values for real and imaginary images. These figures show the mean is dependent on R and I , while the variance depends on R, I, and σ. We can also conclude that the distribution in Eq. (12) can be approximated by a Gaussian for each point. The Gaussian claim for fðθÞ can be quantified using kurtosis and skewness measures. 20 Furthermore, the normality of distribution based on its samples can be classified using the D'Agostino test 20 which is based on kurtosis and skewness.  4 and 5 except for one case, the distributions can be safely approximated by a Gaussian. However in order to find the ranges of R, I, and σ for which the Gaussian approximation holds, the kurtosis, skewness, and D'Agostino test results for various ranges of R and I (between 1 and 100 with step of 5) and σ are shown in Fig. 6. The significance level used for the D'Agostino test is 0.01. It can be seen that the ratio of magnitude ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2 þ I 2 p to σ determines whether the Gaussian assumption holds. As this ratio increases, the distribution gets closer to a Gaussian. In both of our datasets, this ratio is close to 50 for a σ of approximately 10, which means the distribution at each point can be approximated by a Gaussian. We can verify the distribution of noise in the phantom by performing a baseline analysis. 21 In this case, we record real and imaginary images twice when there is no current in each orientation. The current is then calculated by considering repeated images in each orientation as two phase cycles of the CDI method. 21 The resulting current density map in A∕m 2 and its histogram are shown in Figs. 7(a) and 7(b). This  This matches the fact that the convolution of 18 Gaussians (for the case of a 3 × 3 derivative template) is Gaussian itself and, therefore, the final J has a noise which is Gaussian while its mean and variance depend on the R, I, and noise in the real and imaginary images. In all three cases in Fig. 7, the background noise has almost the same zero mean Gaussian distribution with a standard deviation of approximately 10 and the noise in the current density images is also zero mean Gaussian with a standard deviation around 3. For the cases where a current was injected, R and I are not constant because of phase changes induced by the current. However, the magnitude ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2 þ I 2 p is constant for all three cases in Fig. 7 because the phantom has the same homogeneous structure. Looking back at Eq. (12), we can see that if the background noise σ is constant and the magnitude ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2 þ I 2 p is constant, then the distribution fðθÞ only depends on C. Functions C in Eq. (6) for various Rs and Is satisfying a constant R 2 þ I 2 are the same; they are only shifted. This can be shown mathematically as R cosðθ þ ϕÞ þ I sinðθ þ ϕÞ ¼ R 0 cosðθÞ þ I 0 sinðθÞ; (13) where R 0 ¼ R cosðϕÞ þ I sinðϕÞ, I 0 ¼ I cosðϕÞ − R sinðϕÞ, and R 02 þ I 02 ¼ R 2 þ I 2 . This can also be verified by plotting fðθÞ for a constant R 2 þ I 2 which is shown in Fig. 8. Therefore, we can conclude that the standard deviation of the noise in CDIs only depends on magnitude and for a material with an approximately constant magnitude, the noise distribution in the current density has the same Gaussian distribution for all the points. The only difference between such a subject and our phantom is that the current density values are not constant due to the lack of symmetry which translates into Gaussian distributions for current densities with the same standard deviation and the mean determined by the nonnoisy current density value at that point. The standard deviation of this noise can be experimentally found for each magnitude level.

Denoising
As shown in Sec. 2.1, the noise in the final J could be approximated to be Gaussian, to remove the noise from CDI images we evaluated three denoising methods, adaptive Wiener, 14 ramppreserving PM 13 which has been previously used for CDI image denoising, 12 and the BM3D 15 method which was originally designed to remove Gaussian noise. The methods are applied on the initial imaginary and real images (I-R) or on the final current density images (J) and the results are compared. For the Wiener filter and BM3D when applied to real and imaginary images, we use the knowledge of the noise standard deviation estimated from the background of the real and imaginary images. Different block sizes were checked for the Wiener filter when applied on I-R and the best choice which provides the lowest sum of squared error (SSE) was 5 for Dataset 1 and 3 for Dataset 2. For PM smoothing, the best performance was obtained by 10 iterations for Dataset 1 and 6 iterations for Dataset 2, using a gradient modulus threshold of 30. The BM3D algorithm is implemented based on Dabov et al. 15

Results and Discussion
In evaluating the denoising performance of the techniques, SSE on the normalized current values was used. The measured currents were normalized for ease of comparison between the two datasets. To remove the unreliable points from the edges, the final measured Js were eroded on the boundaries. This erosion causes the loss of some area of the phantom in the imaging slices. In this analysis, we lost 19% of the area due to erosion, therefore, the total measured (expected) current was around 16.2 mA for Dataset 1 and 36.5 mA for Dataset 2. Figures 9  and 10 show the box plots of the SSE on the normalized current results for all scenarios in Dataset 1 and Dataset 2, respectively. The data in each box plot are from multiple CDI acquisitions (multiple slices). The measured current densities are normalized by the injected currents' amplitudes in milliamperes and SSE is calculated between the normalized measured current densities and the normalized expected current density. The normalized expected current density is calculated by dividing 1 mA (normalized current) over the area of the phantom slice because of the uniformity of the current density in the phantom. From the SSE results, we can see that overall the BM3D method performs better in denoising CDI images than the compared existing techniques. It performs well both in applying denoising on the real and imaginary images before calculating the phase as well as on the final calculated current density map (J). The performance is better when the BM3D method was applied to the final measured J than when it was only applied to R and I for both datasets. It can be noted that the normalized SSEs have higher values for lower current because of the normalization by the current amplitude; for higher currents the values will be lower. We also computed the P-values to demonstrate the pairwise statistical differences between SSE distributions between the techniques. The statistical significance along with a lower SSE demonstrates the pairwise comparative performance of the techniques. The P-values calculated between different scenarios are shown in Table 2. The P-values were calculated using an ANOVA 1 test 22 while the null hypothesis is that the data of the compared scenarios were obtained from same dataset with the same mean. It can be seen that for both datasets the P-values show a significant statistical difference between PM and BM3D and also between the Wiener and BM3D scenarios. The Pvalues show more distinction for lower current values where we have a lower signal-to-noise ratio (SNR) (Dataset 1). This is expected since at lower SNRs, the real gain of denoising is emphasized and this separates the SSE distributions of the techniques.
To further quantify the gains obtained by applying the BM3D method on J images, we considered the BM3D result as the reference and calculated the average SSE per pixel between BM3D and the other two methods (Wiener and PM). The results are shown in Figs. 11 and 12 for both datasets. The data in each   box plot are from multiple CDI acquisitions (multiple slices). The comparative gain (dB per pixel) in denoising is evident from the results. We also note that the noise suppression is higher for PM in lower SNRs (Dataset 1), while it is higher for Wiener in the higher SNRs (Dataset 2). The performance of BM3D could be attributed to the fact that it uses a sparse representation of images by exploiting the fact that there is structural redundancy and similarity in the images and that it also uses a 3-D denoising method. It reconstructs the denoised image by emphasizing the groups of patches with higher similarity. Considering that our data were acquired using a homogeneous cylindrical phantom, the structural redundancy and similarity might have aided the BM3D to some extent; further analysis in future studies using real world scenarios should be done. However, because BM3D is shown to outperform other denoising methods on standard images, 15  Comparing BM3D with Wiener filtering, we can see that although BM3D also relies on Wiener filtering, it also exploits the image pattern similarities which Wiener filtering lacks. All these denoising methods compromise the spatial resolution for the sake of combating noise. Nevertheless, the approximation of noise distribution in J as Gaussian has provided with the motivation in exploring appropriate denoising methods that would optimally remove the noise influence from the current measurements.

Conclusion
The accuracy of current measurements using the CDI technique is affected by the inherent MRI noise. This noise, which is present in the initial real and imaginary images recorded by the MR scanner, is reflected in the final measured current density map. In this work, we have shown that the current density noise distribution at each pixel can be approximated by a Gaussian that depends on the noiseless values of real and imaginary images at the pixel and the variance of the background noise. Three denoising methods were evaluated at two different stages of the CDI current density measurement algorithm. Based on the obtained results using two CDI datasets of phantom images, the BM3D method which was originally designed for Gaussian noise, performs better than the other two techniques compared in this work. Furthermore, the BM3D method also performs better when applied directly on the final calculated J which simplifies the denoising operation in comparison to applying it on R and I images.