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 phantoms1,3,4 and also applied to the in vivo and ex vivo tissues in order to obtain the current density map inside the subject.56.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-the-art 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.
The experimental data used in this paper were obtained by performing a CDI experiment on a cylindrical phantom filled with and solution with and .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-to-peak 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 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 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 -direction and it was parallel to the direction of the current. The coordinate system was adopted from Scott et al.17 in which the -direction is perpendicular to the main magnetic field. The slices were equally spaced and were perpendicular to the -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 -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 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. 3(a) and 3(b), respectively, with a Gaussian distribution fitted to the samples (solid curve). The Gaussian distribution fitted to Fig. 3 has standard deviations of 10.75 and 10.57 for Figs. 3(a) and 3(b), respectively, and means of and for Figs. 3(a) and 3(b), 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 is the random variable representing the corrupted intensity in real images and is the random variable representing the corrupted intensity in imaginary images, we have for one pixel
is the noise random variable for the chosen pixel in the real image and is the random variable representing noise in the imaginary image for the same pixel. Therefore, and are Gaussian at each pixel with nonzero means and at that pixel. In this case, the joint distribution of and is a bivariate Gaussian distribution given as
The distribution of the angle and magnitude of can be found using the Jacobian of the transformation as
The symbol arg represents the argument of (). Because the current density is measured from the phase changes, we find by integrating over
Changing variable ,
Using the following identities
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 . The unwrapped values are proportional to the induced magnetic flux (). Therefore, we can reason that 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 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 template is used, the nine values of 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 (). 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 , , and . The values chosen for and 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 and , while the variance depends on , , and . We can also conclude that the distribution in Eq. (12) can be approximated by a Gaussian for each point. The Gaussian claim for 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 test20 which is based on kurtosis and skewness. The results for Figs. 4 and 5 are shown in Table 1. The significance level used for the D’Agostino test in this table is 0.01. The table shows that for all the cases investigated in Figs. 4 and 5 except for one case, the distributions can be safely approximated by a Gaussian. However in order to find the ranges of , , and for which the Gaussian approximation holds, the kurtosis, skewness, and D'Agostino test results for various ranges of and (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 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 and its histogram are shown in Figs. 7(a) and 7(b). This histogram is reflective of the noise distribution because there is not much change in real and imaginary images in the absence of injected current, thus the real and imaginary images should be almost constant. Figures 7(c) and 7(d) and 7(e) and 7(f) show the calculated current density maps and their histograms when current is injected for Dataset 1 and Dataset 2, respectively. It can be seen that the noise distribution is only shifted when the current is injected. The red curves on the histograms show the Gaussian distributions fitted on the samples. For the no current case, the mean of the fitted Gaussian is 0.03 and its standard deviation is 3.62. The Gaussian fitted for Dataset 1 has mean of 13.24 and standard deviation of 3.67, while for Dataset 2 the Gaussian has mean of 29.45 and standard deviation of 3.19. This matches the fact that the convolution of 18 Gaussians (for the case of a derivative template) is Gaussian itself and, therefore, the final has a noise which is Gaussian while its mean and variance depend on the , , 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, and are not constant because of phase changes induced by the current. However, the magnitude 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 is constant, then the distribution only depends on . Functions in Eq. (6) for various s and s satisfying a constant are the same; they are only shifted. This can be shown mathematically asFig. 8.
Kurtosis, skewness, and D’Agostino test results for Figs. 4 and 5.
|70||50||24||2.8384||−6.5718 × 10−4||Normal|
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 non-noisy current density value at that point. The standard deviation of this noise can be experimentally found for each magnitude level.
As shown in Sec. 2.1, the noise in the final could be approximated to be Gaussian, to remove the noise from CDI images we evaluated three denoising methods, adaptive Wiener,14 ramp-preserving PM13 which has been previously used for CDI image denoising,12 and the BM3D15 method which was originally designed to remove Gaussian noise. The methods are applied on the initial imaginary and real images () or on the final current density images () 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 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 s 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 (). The performance is better when the BM3D method was applied to the final measured than when it was only applied to and 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 -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 -values calculated between different scenarios are shown in Table 2. The -values were calculated using an ANOVA 1 test22 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 -values show a significant statistical difference between PM and BM3D and also between the Wiener and BM3D scenarios. The -values 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.
P-values of SSEs for various scenarios.
|Compared scenarios||Set 1||Set 2|
|Wiener () and BM3D ()||0.0034|
|Wiener () and BM3D ()||0.0108|
|Wiener () and PM ()||0.2078|
|PM () and PM ()||0.077|
|Wiener () and Wiener ()||0.0661|
|BM3D () and BM3D ()||0.0404||0.2456|
|BM3D () and PM ()||0.039|
|BM3D () and Wiener ()||0.017|
To further quantify the gains obtained by applying the BM3D method on 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 the results can be extended to the application of CDI to tissues. Regardless of the image structure, a better denoising performance is expected from BM3D compared with the two other methods (PM and Wiener). BM3D compared with PM has a wider definition of image patterns. It deals with multiple groups of similar patches instead of only the edge or nonedge patterns dealt with in PM and on top of that it employs Wiener filtering. 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 as Gaussian has provided with the motivation in exploring appropriate denoising methods that would optimally remove the noise influence from the current measurements.
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 which simplifies the denoising operation in comparison to applying it on and images.
We would like to thank Mr. Eugen Hlasny and Mr. Neil Spiller from Toronto General Hospital Medical Imaging Department, who helped us in conducting MRI imaging and data collection. Their technical comments and kind support contributed to this work. We also thank Canadian Institute of Health Research (CIHR) grant number 111253 for funding this research.
Mohammadali Beheshti is a PhD candidate in the Department of Electrical and Computer Engineering at Ryerson University. He received his MSc degree in wireless communications from McMaster University in 2004 and his BSc degree in communications from Isfahan University of Technology in 2002. He has authored and coauthored publications in biomedical engineering and signal processing areas. His current research interest is medical imaging, biomedical system modeling, and signal and image processing.
Farbod H. Foomany is a researcher in signal processing and software/system security. He received his PhD from the University of Huddersfield, UK, in 2010, on analysis of biometrics for crime reduction. He has master's and bachelor's degrees in computer and electrical engineering, respectively. He has authored and coauthored several conference and journal papers in the area of security, crime reduction, and biomedical signal processing. His research interests include applications of biomedical signal processing, security of software, and control systems.
Karl Magtibay obtained his BEng and MASc degrees in biomedical engineering and electrical and computer engineering, in 2012 and 2014, respectively, from Ryerson University in Toronto, Canada. He is a research analyst at University Health Network, Toronto General Hospital. His research interest includes applications of medical imaging and digital signal and imaging processing techniques to cardiac electrophysiology.
David A. Jaffray is the head of the Radiation Physics Department at Princess Margaret Hospital as well as an associate professor in the Departments of Radiation Oncology and Medical Biophysics at the University of Toronto. He has pioneered the development of Cone Beam CT and is the recipient of many research awards. He completed his PhD in the Department of Medical Biophysics at the University of Western Ontario and has a BSc degree in physics from the University of Alberta in 1988.
Sridhar Krishnan received his BE degree in electronics and communication engineering from Anna University in 1993 and the MS and PhD degrees in electrical and computer engineering from the University of Calgary in 1996 and 1999, respectively. He is a professor in the Department of Electrical Engineering at Ryerson University. He has been a Canada research chair in biomedical signal analysis since 2007 and an associate dean for the Faculty of Engineering and Architectural Science since 2011.
Kumaraswamy Nanthakumar received his MD degree from the University of Toronto in 1995. He is a director of heart rhythm disorders at University Health Network and an associate professor of medicine at the University of Toronto. His research interests are the mechanisms of cardiac fibrillation, CPR (cardiopulmonary resuscitation), and ventricular tachycardia. He was the recipient of the Early Researcher Award from the Ministry of Innovation as well as the Clinician Scientist Award from the Canadian Institutes of Health Research.
Karthikeyan Umapathy is an associate professor in the Department of Electrical and Computer Engineering at Ryerson University. He received his BE degree in electronics and communication engineering from Bharathiar University in 1992. He received his MPhil degree in electronics, communication, and electrical engineering from the University of Hertfordshire in 2002 and his PhD in electrical and computer engineering from the University of Western Ontario in 2006. He holds an affiliate scientist position at St. Michael’s Hospital.