Bio-Medical Materials and Engineering 24 (2014) 1239–1245 DOI 10.3233/BME-130925 IOS Press
1239
Bias Correction for Magnetic Resonance Images via Joint Entropy Regularization1 Shanshan Wang a,b, Yong Xia b,e, Pei Dong b, Jianhua Luo c*, Qiu Huang a,d, Dagan Feng b,d and Yuanxiang Lic* a
School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China School of Information Technologies, The University of Sydney, NSW 2006, Australia c School of Aeronautics and Astronautics, Shanghai Jiao Tong University, Shanghai 200240, China d Med-X Research Institute, Shanghai Jiao Tong University, Shanghai 200030, China e SAIIP, School of Computer Sciences, Northwestern Polytechnical University, Xi’an 710072, China b
Abstract. Due to the imperfections of the radio frequency (RF) coil or object-dependent electrodynamic interactions, magnetic resonance (MR) images often suffer from a smooth and biologically meaningless bias field, which causes severe troubles for subsequent processing and quantitative analysis. To effectively restore the original signal, this paper simultaneously exploits the spatial and gradient features of the corrupted MR images for bias correction via the joint entropy regularization. With both isotropic and anisotropic total variation (TV) considered, two nonparametric bias correction algorithms have been proposed, namely IsoTVBiasC and AniTVBiasC. These two methods have been applied to simulated images under various noise levels and bias field corruption and also tested on real MR data. The test results show that the proposed two methods can effectively remove the bias field and also present comparable performance compared to the state-of-the-art methods. Keywords: Bias correction, magnetic resonance (MR) images, joint entropy, total variation (TV)
1. Introduction Magnetic resonance images are often contaminated by a smooth and biologically meaningless bias field due to quite a few factors such as poor radio frequency coil uniformity, gradient-driven eddy currents and object-dependent electrodynamic interactions [1]. The bias field, which although is not a big problem for clinical visual inspection, can severely impede the quantitative analysis and automatic post-processing (like registration and segmentation) of MR images. Therefore, it is essential to develop an automatic and inexpensive algorithm for bias field removal. In the recent decades, there have been extensive literatures published and algorithms proposed, which can be roughly separated into two directions i.e. prospective and retrospective techniques. The prospective methods treat the bias artifact as a systematic error introduced during the MRI acquisition. 1 This work was supported in part by the ARC grants and in part by the China Scholarship Council under Grant 2011623084. *Corresponding author. E-mail: [email protected]. *Corresponding author. E-mail: [email protected] 0959-2989/14/$27.50 © 2014 – IOS Press and the authors. All rights reserved
1240
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
Therefore, this type of methods focuses on designing hardware acquisition sequences and calibrating the image system with additional information (from uniform phantom) or with different coils. On the other hand, the retrospective techniques attribute high credits to the information of the acquired image and the prior knowledge. This category is more general due to two points: 1) fewer assumptions are made about the acquisition process; and 2) patient dependant inhomogeneity can also be removed. Retrospective methods analyze the observed signal which is an integration of the imaged anatomy and the intensity inhomogeneity. To facilitate the extraction of the original signal, researchers commonly explore a prior knowledge, for instance, the spatial or intensity probability distribution of the anatomical information. Based on the prior knowledge and the operand, the retrospective approaches can be further divided into four sub-directions, namely, filtering methods, surface fitting, segmentation-based and Histogram-based methods [2]. Filtering methods assume that the bias field can be separated from the imaged anatomy via low-pass filtering [3]. This assumption is quite limited since most anatomy may possess low-frequency information which can be mistakenly removed by filters. Surface fitting method is quite popular due to its robustness to image noise. This kind of techniques fits a parametric surface to the image features. The bias is generally modelled as a polynomial or spline based multiplicative field [4; 5]. Segmentation based methods simultaneously remove the bias field and segment the MR images so that they can benefit from the result obtained by each other [6]. However, most of the above mentioned methods rely on the proper initialization and the priori knowledge on the imaged structures. Histogram based methods, on the other hand, analyse directly the image intensity histograms or minimization of the image entropy [7-9]. This makes the histogram based methods quite automatic and highly generally. In [7], a new method of automatically correcting the bias field is proposed, which we name it as IsoCeDBiasC, for short. This method combines the image intensity and gradient image features to formulate a new entropy-related cost function. However, this method only utilizes the central isotropic gradient features of the corrupted image, which hasn’t been proved to best suit the MR images. Therefore, this paper tries to explore the forward total variation (TV) of the MR images. With both isotropic and anisotropic TV features considered, two algorithms have been proposed. The original method proposed in [7] is also extended by considering the anisotropic central gradient features, AniCeDBiasC, for short. The proposed methods have been compared to three state-of-the-art algorithms, namely, IsoCeDBiasC, AniCeDBiasC, and the recent one named iteratively re-weighted least squares (IRLS) based algorithm [5]. Both simulated and real data have been tested on under a variety of noise levels and bias levels. 2. Methodology With ܷ denoting the original MR image, the corrupted data ܻ can be modelled as [2; 7] ܻሺݔǡ ݕሻ ൌ ܷሺݔǡ ݕሻܤሺݔǡ ݕሻ ݊ሺݔǡ ݕሻ
(1)
where ܤis the corresponding value of the smooth bias field, ݊ is the Rician distributed additive noise and ሺݔǡ ݕሻ is pixel index in the image. Let ܪand ܹ respectively denote the height and the width of the image, we have א ݔሾͳǡ ܪሿ and א ݕሾͳǡ ܹሿ. The goal is to restore ܷሺݔǡ ݕሻ from ܻሺݔǡ ݕሻ.
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
1241
2.1. Preprocessing The noise may introduce the dispersion of image intensity values in uniform areas. To avoid the impact of the noise, a simple low-pass filter is utilized to pre-process the image. Next, the foreground is identified according to the well-known Otsu method [10]. The background therefore can be removed prior to the bias correction since it mainly occurs at the imaged anatomy as instructed by [8]. 2.2. Bias removal After the preprocessing, the bias removal is conducted under the regularization of the joint entropy. The entropy should be used to measure the amount of the information inherited in the observation and the corrected image. Therefore the cost function of the entropy is expected to possess the ability to determine the image homogeneity. 2.2.1. Joint entropy formulation Shannon entropy is widely used for uniformity or homogeneity measurement [11], whose formulation can be given as follows
బ ௐ ܪଵ ൌ െ σୀଵ , where ൌ σு ௫ୀଵ σ௬ୀଵ ߜ ሼܻሺݔǡ ݕሻ െ ݅ ሽ ǡܮ א ݅
(2)
where ܮ is the grey levels of the corrupted image and ߜ is the Kronecker delta function. This simple entropy has been extended by Manjón [7] with the image local gradients, based on the assumption that uniform images will simultaneously present well-ordered intensities and wellclustered very low gradient values. The joint entropy ܪଶ is formulated as
బ భ σୀଵ ܪଶ ൌ െ σୀଵ where ൌ
ಽ
ೕ
(3)
ಽ
బ σ భ σసభ ೕసభ ೕ
ௐ ൌ σு ௫ୀଵ σ௬ୀଵ ߜ ሼܻሺݔǡ ݕሻ െ ݅ ሽߜ ሼܩሺݔǡ ݕሻ െ ݆ሽ ǡ ܮ א ݅ ǡ ܮ א ݆ଵ
(4)
ܩሺݔǡ ݕሻ is the associated magnitude image of the local gradient with ܮଵ grey levels. 2.2.2. The gradient features According to (4), it can be seen that the gradient space ܩሺݔǡ ݕሻ plays an essential role in the joint entropy and therefore its formulation is quite significant. Here, we consider four types of gradient features, the isotropic and anisotropic TV of the corrupted image, and the isotropic and anisotropic central difference of the corrupted image. Their respective mathematical equation is given in Eqs. (5-8). ଶ
ଶ
୍ܩୱ୭ ሺݔǡ ݕሻ ൌ ට൫ܻሺ ݔ ͳǡ ݕሻ െ ܻሺݔǡ ݕሻ൯ ൫ܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕሻ൯
(5)
ܩ୬୧ ሺݔǡ ݕሻ ൌ ȁܻሺ ݔ ͳǡ ݕሻ െ ܻሺݔǡ ݕሻȁ ȁܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕሻȁ
(6)
In the application of image denoising, isotropic TV settings have been shown better than anisotropic models in terms of rotational invariance. The anisotropic has a qualitative improvement at corners [12]. And it also has been shown that the introduction of anisotropy to TV regularization can improve the denoising results. Specifically, it reduces the staircasing effect and meanwhile creates less artifact. ଶ
ଶ
ଵ ୍ܩୱ୭େୣୈ ሺݔǡ ݕሻ ൌ ට൫ܻሺ ݔ ͳǡ ݕሻ െ ܻሺ ݔെ ͳǡ ݕሻ൯ ൫ܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕെ ͳሻ൯ ଶ
(7)
1242
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization ଵ
ଵ
ଶ
ଶ
ܩ୬୧େୣୈ ሺݔǡ ݕሻ ൌ ȁܻሺ ݔ ͳǡ ݕሻ െ ܻሺ ݔെ ͳǡ ݕሻȁ ȁܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕെ ͳሻȁ
(8)
The isotropic central gradient difference was originally used in the implementation of the algorithm presented in [7]. In Fig. 1, an example has been provided. As you can see, when the image is corrupted by the bias field, the image has some staircasing artifacts in all the four kinds of difference images compared to them of the bias free images. It can also be observed that the isotropic and anisotropic central difference images are smoother than both types of TV images.
AniTV(0%)
AniTV(40%)
AniCeD(0%)
AniCeD(40%)
IsoTV(0%)
IsoTV(40%)
IsoCeD(0%)
IsoCeD(40%)
Fig. 1. Example of the four types of gradient features for an axial PD-weighted image without bias and with 40% bias field.
2.2.3. Bias formulation In order to remove the bias field, like [7], we also model it as the linear sum of ݇ equidistant Bspline functions ݂ , which are low frequency basis functions. Mathematically, it can be denoted as ܤሺܽሻ ൌ σୀଵ ܽ ݂ , Ǥ Ǥ σୀଵ ݂ ሺݔ௩ ሻ ൌ ͳǡ א ݒሾͳǡ ܹܪሿ
(9)
where ݔis the vector of all image pixel locations of size ܹܪ. The bias field is determined by the value of ܽ, which can be obtained according to ܤൌ ܪଶ ൫൫ܻȀܤሺܽሻ൯൯
(10)
Following the common trend, in our implementation, the mean value of the corrected image was preserved as the same value of the signal so as to preserve the original contrast of the image [8].
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
1243
2.3. Image update The value of the coefficient matrixes are estimated according to (10) at each iteration and therefore the bias field can be calculated. Based on whether the change rate of the entropy formulation satisfies a previous threshold or not, the scale of the coefficients will be further adjusted in the iterations. Once it reaches the threshold, the bias will be formulated to the original dimension and the original image is obtained with the final estimated bias. 3. Experiments 3.1. Correcting simulated data The proposed algorithms have been compared to IsoCeDBiasC, AniCeDBiasC, and the recent IRLS based algorithm2, on the 3D MR volumes acquired from the BrainWeb simulated brain database at the McConnel Brain Imaging Center of the Montreal Neurological Institute, McGill University, Montreal, Canada [13]. As shown in Fig. 2, three imaging modalities have been considered, i.e. T1-, T2-, and PD-weighted volumes, each of which has the same dimension of ͳͺͳ ൈ ʹͳ ൈ ͳͺͳ and a voxel size of ͳ ൈ ͳ ൈ ͳଷ . The coronal, sagittal and axial planes were respectively extracted from these volumes with most anatomy. The test images were corrupted with two levels of bias field and five levels of noise. The corresponding image without bias field and noise were also downloaded and used as the ground truth to quantitatively evaluate the five algorithms.
(a) Sagittal T1
(b) Coronal T2
(c) Axial PD
Fig. 2. Examples of simulated test images
The results produced by our proposed two algorithms and three competing algorithms on bias correction of the simulated images are quantitatively evaluated by using the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) index (in brackets) [14], and presented in Table 1 respective෩ is defined as ൌ ͳͲ ଵ ሺʹͷͷଶ Ȁሻ where is the ly. PSNR of a bias corrected image ܷ ෩. The proposed IsoTVBiasC and AniTVBiasC algorithms outpermean squared error between ܷ and ܷ formed the other three algorithms in 8 out of 12 cases (shown in Table 2), regardless of the image quality metric used. The exploitation of isotropic TV yielded the best overall results among all five 2
http://www.mathworks.com/matlabcentral/fileexchange/27315. The default configuration was used in our experiments. No preprocessing is applied to the test image for IRLS such as background removal.
1244
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
algorithms, while our anisotropic TV-based algorithm is also capable of effectively removing the bias field compared to the IsoCeDBiasC and AniCeDBiasC. Table 1 PSNR/SSIM results of T1-weighted brain MR images after bias correction by all five algorithms. The best result in each case is highlighted in boldface. Bias
20%
40%
Noise 0 1 3 5 7 9 0 1 3 5 7 9
IsoCeDBiasC 30.138 (0.9925) 28.600 (0.9587) 25.589 (0.8569) 23.197 (0.7682) 22.259 (0.6890) 20.156 (0.6195) 25.981 (0.9826) 25.526 (0.9532) 23.554 (0.8490) 21.531 (0.7577) 20.008 (0.6858) 19.186 (0.6193)
AniCeDBiasC 29.401 (0.9919) 28.576 (0.9579) 25.408 (0.8569) 23.213 (0.7694) 22.311 (0.6908) 20.100 (0.6217) 25.920 (0.9823) 25.606 (0.9527) 23.741 (0.8539) 21.666 (0.7620) 20.145 (0.6888) 19.072 (0.6168)
IRLS-based 13.130 (0.6625) 14.448 (0.6740) 14.996 (0.6542) 15.007 (0.6165) 15.561 (0.6007) 15.485 (0.5564) 13.068 (0.6531) 14.186 (0.6541) 15.036 (0.6534) 15.016 (0.6212) 15.532 (0.6000) 14.802 (0.5419) Table 2
IsoTVBiasC 31.217 (0.9935) 29.915 (0.9613) 25.781 (0.8605) 23.191 (0.7665) 22.160 (0.6870) 20.271 (0.6228) 26.072 (0.9824) 25.512 (0.9540) 23.651 (0.8542) 21.716 (0.7649) 20.204 (0.6895) 18.935 (0.6137)
AniTVBiasC 31.647 (0.9942) 28.103 (0.9573) 25.271 (0.8572) 23.163 (0.7665) 22.096 (0.6886) 20.153 (0.6234) 25.925 (0.9822) 25.160 (0.9534) 23.758 (0.8534) 21.501 (0.7592) 20.212 (0.6892) 19.076 (0.6168)
Numbers of winning cases and win rates of five bias correction methods on simulated brain MR images. The statistics are obtained from all three test images, two bias levels and six noise levels. Metric PSNR SSIM
IsoCeDBiasC 10 (27.78%) 9 (25.00%)
AniCeDBiasC 8 (22.22%) 8 (22.22%)
IRLS-based 0 (0.00%) 1 (2.78%)
IsoTVBiasC 8 (22.22%) 10 (27.78%)
AniTVBiasC 10 (27.78%) 8 (22.22%)
3.2. Correcting real MR data Finally, the proposed algorithms were evaluated on the real MR data which were obtained from the American Radiology Services3. The bias field estimated by these algorithms have been presented in Fig. 3, which is an example of the ankle. It can be seen the bias field is quite smooth and most parts only cover the anatomical structure. We present the histogram of the bias corrupted and corrected images here in Fig. 3.
(a) Ankle
(b) IsoTVBiasC
(c) AniTVBiasC
(d) Histograms
Fig. 3. Correcting the bias-corrupted clinical image (a) Ankle: the bias fields estimated by the proposed (b) IsoTVBiasC and (c) AniTVBiasC algorithms, and the (d) histograms of the corrupted image and bias corrected.
3
http://www3.americanradiology.com/
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
1245
4. Conclusions The proposed two methods can be regarded as an extension to the nonparametric MRI inhomogeneity correction method [7], which only considers smooth local gradient information for the joint entropy regularization in its online implementation package. While in our method, both isotropic and anisotropic TV have been explored. The proposed algorithms have been tested on both simulated and clinical MR images. The experiments present results comparable to state-of-the-art bias correction approaches. References [1] J. Luo, Y. Zhu, P. Clarysse, and I.E. Magnin, Correction of bias field in MR images using singularity function analysis, IEEE Trans. Med. Imaging 24 (2005), 1067-1085. [2] U. Vovk, F. Pernus, and B. Likar, A review of methods for correction of intensity inhomogeneity in MRI, IEEE Trans. Med. Imaging 26 (2007), 405-421. [3] B.H. Brinkmann, A. Manduca, and R.A. Robb, Optimized homomorphic unsharp masking for MR grayscale inhomogeneity correction, IEEE Trans. Med. Imaging 17 (1998), 161-171. [4] Y. Zheng and J.C. Gee, Estimation of image bias field with sparsity constraints, in: Proceedings of the Twenty-Third IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2010, pp. 255-262. [5] Y. Zheng, M. Grossman, S.P. Awate, and J.C. Gee, Automatic correction of intensity nonuniformity from sparseness of gradient distribution in medical images, in: Medical Image Computing and Computer-Assisted Intervention (MICCAI), G.-Z. Yang, D. Hawkes, D. Rueckert, A. Noble, and C. Taylor, eds., Springer Berlin Heidelberg, 2009, pp. 852-859. [6] E.A. Vokurka, N.A. Watson, Y. Watson, N.A. Thacker, and A. Jackson, Improved high resolution MR imaging for surface coils using automated intensity non-uniformity correction: feasibility study in the orbit, Journal of Magnetic Resonance Imaging 14 (2001), 540-546. [7] J.V. Manjón, J.J. Lull, J. Carbonell-Caballero, G. García-Martí, L. Martí-Bonmatí, and M. Robles, A nonparametric MRI inhomogeneity correction method, Medical image analysis 11 (2007), 336-345. [8] B. Likar, M.A. Viergever, and F. Pernus, Retrospective correction of MR intensity inhomogeneity by information minimization, IEEE Trans. Med. Imaging 20 (2001), 1398-1410. [9] B. Likar, J.A. Maintz, M.A. Viergever, and F. Pernus, Retrospective shading correction based on entropy minimization, Journal of Microscopy 197 (2000), 285-295. [10] N. Otsu, A threshold selection method from gray-level histograms, Automatica 11 (1975), 23-27. [11] C.E. Shannon, A mathematical theory of communication, ACM SIGMOBILE Mobile Computing and Communications Review 5 (2001), 3-55. [12] H. Birkholz, A unifying approach to isotropic and anisotropic total variation denoising models, Journal of computational and applied mathematics 235 (2011), 2502-2514. [13] D.L. Collins, A.P. Zijdenbos, V. Kollokian, J.G. Sled, N.J. Kabani, C.J. Holmes, and A.C. Evans, Design and construction of a realistic digital brain phantom, IEEE Trans. Med. Imaging 17 (1998), 463-468. [14] Z. Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Trans. Image Process. 13 (2004), 600-612.
Copyright of Bio-Medical Materials & Engineering is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
1239
Bias Correction for Magnetic Resonance Images via Joint Entropy Regularization1 Shanshan Wang a,b, Yong Xia b,e, Pei Dong b, Jianhua Luo c*, Qiu Huang a,d, Dagan Feng b,d and Yuanxiang Lic* a
School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China School of Information Technologies, The University of Sydney, NSW 2006, Australia c School of Aeronautics and Astronautics, Shanghai Jiao Tong University, Shanghai 200240, China d Med-X Research Institute, Shanghai Jiao Tong University, Shanghai 200030, China e SAIIP, School of Computer Sciences, Northwestern Polytechnical University, Xi’an 710072, China b
Abstract. Due to the imperfections of the radio frequency (RF) coil or object-dependent electrodynamic interactions, magnetic resonance (MR) images often suffer from a smooth and biologically meaningless bias field, which causes severe troubles for subsequent processing and quantitative analysis. To effectively restore the original signal, this paper simultaneously exploits the spatial and gradient features of the corrupted MR images for bias correction via the joint entropy regularization. With both isotropic and anisotropic total variation (TV) considered, two nonparametric bias correction algorithms have been proposed, namely IsoTVBiasC and AniTVBiasC. These two methods have been applied to simulated images under various noise levels and bias field corruption and also tested on real MR data. The test results show that the proposed two methods can effectively remove the bias field and also present comparable performance compared to the state-of-the-art methods. Keywords: Bias correction, magnetic resonance (MR) images, joint entropy, total variation (TV)
1. Introduction Magnetic resonance images are often contaminated by a smooth and biologically meaningless bias field due to quite a few factors such as poor radio frequency coil uniformity, gradient-driven eddy currents and object-dependent electrodynamic interactions [1]. The bias field, which although is not a big problem for clinical visual inspection, can severely impede the quantitative analysis and automatic post-processing (like registration and segmentation) of MR images. Therefore, it is essential to develop an automatic and inexpensive algorithm for bias field removal. In the recent decades, there have been extensive literatures published and algorithms proposed, which can be roughly separated into two directions i.e. prospective and retrospective techniques. The prospective methods treat the bias artifact as a systematic error introduced during the MRI acquisition. 1 This work was supported in part by the ARC grants and in part by the China Scholarship Council under Grant 2011623084. *Corresponding author. E-mail: [email protected]. *Corresponding author. E-mail: [email protected] 0959-2989/14/$27.50 © 2014 – IOS Press and the authors. All rights reserved
1240
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
Therefore, this type of methods focuses on designing hardware acquisition sequences and calibrating the image system with additional information (from uniform phantom) or with different coils. On the other hand, the retrospective techniques attribute high credits to the information of the acquired image and the prior knowledge. This category is more general due to two points: 1) fewer assumptions are made about the acquisition process; and 2) patient dependant inhomogeneity can also be removed. Retrospective methods analyze the observed signal which is an integration of the imaged anatomy and the intensity inhomogeneity. To facilitate the extraction of the original signal, researchers commonly explore a prior knowledge, for instance, the spatial or intensity probability distribution of the anatomical information. Based on the prior knowledge and the operand, the retrospective approaches can be further divided into four sub-directions, namely, filtering methods, surface fitting, segmentation-based and Histogram-based methods [2]. Filtering methods assume that the bias field can be separated from the imaged anatomy via low-pass filtering [3]. This assumption is quite limited since most anatomy may possess low-frequency information which can be mistakenly removed by filters. Surface fitting method is quite popular due to its robustness to image noise. This kind of techniques fits a parametric surface to the image features. The bias is generally modelled as a polynomial or spline based multiplicative field [4; 5]. Segmentation based methods simultaneously remove the bias field and segment the MR images so that they can benefit from the result obtained by each other [6]. However, most of the above mentioned methods rely on the proper initialization and the priori knowledge on the imaged structures. Histogram based methods, on the other hand, analyse directly the image intensity histograms or minimization of the image entropy [7-9]. This makes the histogram based methods quite automatic and highly generally. In [7], a new method of automatically correcting the bias field is proposed, which we name it as IsoCeDBiasC, for short. This method combines the image intensity and gradient image features to formulate a new entropy-related cost function. However, this method only utilizes the central isotropic gradient features of the corrupted image, which hasn’t been proved to best suit the MR images. Therefore, this paper tries to explore the forward total variation (TV) of the MR images. With both isotropic and anisotropic TV features considered, two algorithms have been proposed. The original method proposed in [7] is also extended by considering the anisotropic central gradient features, AniCeDBiasC, for short. The proposed methods have been compared to three state-of-the-art algorithms, namely, IsoCeDBiasC, AniCeDBiasC, and the recent one named iteratively re-weighted least squares (IRLS) based algorithm [5]. Both simulated and real data have been tested on under a variety of noise levels and bias levels. 2. Methodology With ܷ denoting the original MR image, the corrupted data ܻ can be modelled as [2; 7] ܻሺݔǡ ݕሻ ൌ ܷሺݔǡ ݕሻܤሺݔǡ ݕሻ ݊ሺݔǡ ݕሻ
(1)
where ܤis the corresponding value of the smooth bias field, ݊ is the Rician distributed additive noise and ሺݔǡ ݕሻ is pixel index in the image. Let ܪand ܹ respectively denote the height and the width of the image, we have א ݔሾͳǡ ܪሿ and א ݕሾͳǡ ܹሿ. The goal is to restore ܷሺݔǡ ݕሻ from ܻሺݔǡ ݕሻ.
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
1241
2.1. Preprocessing The noise may introduce the dispersion of image intensity values in uniform areas. To avoid the impact of the noise, a simple low-pass filter is utilized to pre-process the image. Next, the foreground is identified according to the well-known Otsu method [10]. The background therefore can be removed prior to the bias correction since it mainly occurs at the imaged anatomy as instructed by [8]. 2.2. Bias removal After the preprocessing, the bias removal is conducted under the regularization of the joint entropy. The entropy should be used to measure the amount of the information inherited in the observation and the corrected image. Therefore the cost function of the entropy is expected to possess the ability to determine the image homogeneity. 2.2.1. Joint entropy formulation Shannon entropy is widely used for uniformity or homogeneity measurement [11], whose formulation can be given as follows
బ ௐ ܪଵ ൌ െ σୀଵ , where ൌ σு ௫ୀଵ σ௬ୀଵ ߜ ሼܻሺݔǡ ݕሻ െ ݅ ሽ ǡܮ א ݅
(2)
where ܮ is the grey levels of the corrupted image and ߜ is the Kronecker delta function. This simple entropy has been extended by Manjón [7] with the image local gradients, based on the assumption that uniform images will simultaneously present well-ordered intensities and wellclustered very low gradient values. The joint entropy ܪଶ is formulated as
బ భ σୀଵ ܪଶ ൌ െ σୀଵ where ൌ
ಽ
ೕ
(3)
ಽ
బ σ భ σసభ ೕసభ ೕ
ௐ ൌ σு ௫ୀଵ σ௬ୀଵ ߜ ሼܻሺݔǡ ݕሻ െ ݅ ሽߜ ሼܩሺݔǡ ݕሻ െ ݆ሽ ǡ ܮ א ݅ ǡ ܮ א ݆ଵ
(4)
ܩሺݔǡ ݕሻ is the associated magnitude image of the local gradient with ܮଵ grey levels. 2.2.2. The gradient features According to (4), it can be seen that the gradient space ܩሺݔǡ ݕሻ plays an essential role in the joint entropy and therefore its formulation is quite significant. Here, we consider four types of gradient features, the isotropic and anisotropic TV of the corrupted image, and the isotropic and anisotropic central difference of the corrupted image. Their respective mathematical equation is given in Eqs. (5-8). ଶ
ଶ
୍ܩୱ୭ ሺݔǡ ݕሻ ൌ ට൫ܻሺ ݔ ͳǡ ݕሻ െ ܻሺݔǡ ݕሻ൯ ൫ܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕሻ൯
(5)
ܩ୬୧ ሺݔǡ ݕሻ ൌ ȁܻሺ ݔ ͳǡ ݕሻ െ ܻሺݔǡ ݕሻȁ ȁܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕሻȁ
(6)
In the application of image denoising, isotropic TV settings have been shown better than anisotropic models in terms of rotational invariance. The anisotropic has a qualitative improvement at corners [12]. And it also has been shown that the introduction of anisotropy to TV regularization can improve the denoising results. Specifically, it reduces the staircasing effect and meanwhile creates less artifact. ଶ
ଶ
ଵ ୍ܩୱ୭େୣୈ ሺݔǡ ݕሻ ൌ ට൫ܻሺ ݔ ͳǡ ݕሻ െ ܻሺ ݔെ ͳǡ ݕሻ൯ ൫ܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕെ ͳሻ൯ ଶ
(7)
1242
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization ଵ
ଵ
ଶ
ଶ
ܩ୬୧େୣୈ ሺݔǡ ݕሻ ൌ ȁܻሺ ݔ ͳǡ ݕሻ െ ܻሺ ݔെ ͳǡ ݕሻȁ ȁܻሺݔǡ ݕ ͳሻ െ ܻሺݔǡ ݕെ ͳሻȁ
(8)
The isotropic central gradient difference was originally used in the implementation of the algorithm presented in [7]. In Fig. 1, an example has been provided. As you can see, when the image is corrupted by the bias field, the image has some staircasing artifacts in all the four kinds of difference images compared to them of the bias free images. It can also be observed that the isotropic and anisotropic central difference images are smoother than both types of TV images.
AniTV(0%)
AniTV(40%)
AniCeD(0%)
AniCeD(40%)
IsoTV(0%)
IsoTV(40%)
IsoCeD(0%)
IsoCeD(40%)
Fig. 1. Example of the four types of gradient features for an axial PD-weighted image without bias and with 40% bias field.
2.2.3. Bias formulation In order to remove the bias field, like [7], we also model it as the linear sum of ݇ equidistant Bspline functions ݂ , which are low frequency basis functions. Mathematically, it can be denoted as ܤሺܽሻ ൌ σୀଵ ܽ ݂ , Ǥ Ǥ σୀଵ ݂ ሺݔ௩ ሻ ൌ ͳǡ א ݒሾͳǡ ܹܪሿ
(9)
where ݔis the vector of all image pixel locations of size ܹܪ. The bias field is determined by the value of ܽ, which can be obtained according to ܤൌ ܪଶ ൫൫ܻȀܤሺܽሻ൯൯
(10)
Following the common trend, in our implementation, the mean value of the corrected image was preserved as the same value of the signal so as to preserve the original contrast of the image [8].
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
1243
2.3. Image update The value of the coefficient matrixes are estimated according to (10) at each iteration and therefore the bias field can be calculated. Based on whether the change rate of the entropy formulation satisfies a previous threshold or not, the scale of the coefficients will be further adjusted in the iterations. Once it reaches the threshold, the bias will be formulated to the original dimension and the original image is obtained with the final estimated bias. 3. Experiments 3.1. Correcting simulated data The proposed algorithms have been compared to IsoCeDBiasC, AniCeDBiasC, and the recent IRLS based algorithm2, on the 3D MR volumes acquired from the BrainWeb simulated brain database at the McConnel Brain Imaging Center of the Montreal Neurological Institute, McGill University, Montreal, Canada [13]. As shown in Fig. 2, three imaging modalities have been considered, i.e. T1-, T2-, and PD-weighted volumes, each of which has the same dimension of ͳͺͳ ൈ ʹͳ ൈ ͳͺͳ and a voxel size of ͳ ൈ ͳ ൈ ͳଷ . The coronal, sagittal and axial planes were respectively extracted from these volumes with most anatomy. The test images were corrupted with two levels of bias field and five levels of noise. The corresponding image without bias field and noise were also downloaded and used as the ground truth to quantitatively evaluate the five algorithms.
(a) Sagittal T1
(b) Coronal T2
(c) Axial PD
Fig. 2. Examples of simulated test images
The results produced by our proposed two algorithms and three competing algorithms on bias correction of the simulated images are quantitatively evaluated by using the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) index (in brackets) [14], and presented in Table 1 respective෩ is defined as ൌ ͳͲ ଵ ሺʹͷͷଶ Ȁሻ where is the ly. PSNR of a bias corrected image ܷ ෩. The proposed IsoTVBiasC and AniTVBiasC algorithms outpermean squared error between ܷ and ܷ formed the other three algorithms in 8 out of 12 cases (shown in Table 2), regardless of the image quality metric used. The exploitation of isotropic TV yielded the best overall results among all five 2
http://www.mathworks.com/matlabcentral/fileexchange/27315. The default configuration was used in our experiments. No preprocessing is applied to the test image for IRLS such as background removal.
1244
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
algorithms, while our anisotropic TV-based algorithm is also capable of effectively removing the bias field compared to the IsoCeDBiasC and AniCeDBiasC. Table 1 PSNR/SSIM results of T1-weighted brain MR images after bias correction by all five algorithms. The best result in each case is highlighted in boldface. Bias
20%
40%
Noise 0 1 3 5 7 9 0 1 3 5 7 9
IsoCeDBiasC 30.138 (0.9925) 28.600 (0.9587) 25.589 (0.8569) 23.197 (0.7682) 22.259 (0.6890) 20.156 (0.6195) 25.981 (0.9826) 25.526 (0.9532) 23.554 (0.8490) 21.531 (0.7577) 20.008 (0.6858) 19.186 (0.6193)
AniCeDBiasC 29.401 (0.9919) 28.576 (0.9579) 25.408 (0.8569) 23.213 (0.7694) 22.311 (0.6908) 20.100 (0.6217) 25.920 (0.9823) 25.606 (0.9527) 23.741 (0.8539) 21.666 (0.7620) 20.145 (0.6888) 19.072 (0.6168)
IRLS-based 13.130 (0.6625) 14.448 (0.6740) 14.996 (0.6542) 15.007 (0.6165) 15.561 (0.6007) 15.485 (0.5564) 13.068 (0.6531) 14.186 (0.6541) 15.036 (0.6534) 15.016 (0.6212) 15.532 (0.6000) 14.802 (0.5419) Table 2
IsoTVBiasC 31.217 (0.9935) 29.915 (0.9613) 25.781 (0.8605) 23.191 (0.7665) 22.160 (0.6870) 20.271 (0.6228) 26.072 (0.9824) 25.512 (0.9540) 23.651 (0.8542) 21.716 (0.7649) 20.204 (0.6895) 18.935 (0.6137)
AniTVBiasC 31.647 (0.9942) 28.103 (0.9573) 25.271 (0.8572) 23.163 (0.7665) 22.096 (0.6886) 20.153 (0.6234) 25.925 (0.9822) 25.160 (0.9534) 23.758 (0.8534) 21.501 (0.7592) 20.212 (0.6892) 19.076 (0.6168)
Numbers of winning cases and win rates of five bias correction methods on simulated brain MR images. The statistics are obtained from all three test images, two bias levels and six noise levels. Metric PSNR SSIM
IsoCeDBiasC 10 (27.78%) 9 (25.00%)
AniCeDBiasC 8 (22.22%) 8 (22.22%)
IRLS-based 0 (0.00%) 1 (2.78%)
IsoTVBiasC 8 (22.22%) 10 (27.78%)
AniTVBiasC 10 (27.78%) 8 (22.22%)
3.2. Correcting real MR data Finally, the proposed algorithms were evaluated on the real MR data which were obtained from the American Radiology Services3. The bias field estimated by these algorithms have been presented in Fig. 3, which is an example of the ankle. It can be seen the bias field is quite smooth and most parts only cover the anatomical structure. We present the histogram of the bias corrupted and corrected images here in Fig. 3.
(a) Ankle
(b) IsoTVBiasC
(c) AniTVBiasC
(d) Histograms
Fig. 3. Correcting the bias-corrupted clinical image (a) Ankle: the bias fields estimated by the proposed (b) IsoTVBiasC and (c) AniTVBiasC algorithms, and the (d) histograms of the corrupted image and bias corrected.
3
http://www3.americanradiology.com/
S. Wang et al. / Bias correction for magnetic resonance images via joint entropy regularization
1245
4. Conclusions The proposed two methods can be regarded as an extension to the nonparametric MRI inhomogeneity correction method [7], which only considers smooth local gradient information for the joint entropy regularization in its online implementation package. While in our method, both isotropic and anisotropic TV have been explored. The proposed algorithms have been tested on both simulated and clinical MR images. The experiments present results comparable to state-of-the-art bias correction approaches. References [1] J. Luo, Y. Zhu, P. Clarysse, and I.E. Magnin, Correction of bias field in MR images using singularity function analysis, IEEE Trans. Med. Imaging 24 (2005), 1067-1085. [2] U. Vovk, F. Pernus, and B. Likar, A review of methods for correction of intensity inhomogeneity in MRI, IEEE Trans. Med. Imaging 26 (2007), 405-421. [3] B.H. Brinkmann, A. Manduca, and R.A. Robb, Optimized homomorphic unsharp masking for MR grayscale inhomogeneity correction, IEEE Trans. Med. Imaging 17 (1998), 161-171. [4] Y. Zheng and J.C. Gee, Estimation of image bias field with sparsity constraints, in: Proceedings of the Twenty-Third IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2010, pp. 255-262. [5] Y. Zheng, M. Grossman, S.P. Awate, and J.C. Gee, Automatic correction of intensity nonuniformity from sparseness of gradient distribution in medical images, in: Medical Image Computing and Computer-Assisted Intervention (MICCAI), G.-Z. Yang, D. Hawkes, D. Rueckert, A. Noble, and C. Taylor, eds., Springer Berlin Heidelberg, 2009, pp. 852-859. [6] E.A. Vokurka, N.A. Watson, Y. Watson, N.A. Thacker, and A. Jackson, Improved high resolution MR imaging for surface coils using automated intensity non-uniformity correction: feasibility study in the orbit, Journal of Magnetic Resonance Imaging 14 (2001), 540-546. [7] J.V. Manjón, J.J. Lull, J. Carbonell-Caballero, G. García-Martí, L. Martí-Bonmatí, and M. Robles, A nonparametric MRI inhomogeneity correction method, Medical image analysis 11 (2007), 336-345. [8] B. Likar, M.A. Viergever, and F. Pernus, Retrospective correction of MR intensity inhomogeneity by information minimization, IEEE Trans. Med. Imaging 20 (2001), 1398-1410. [9] B. Likar, J.A. Maintz, M.A. Viergever, and F. Pernus, Retrospective shading correction based on entropy minimization, Journal of Microscopy 197 (2000), 285-295. [10] N. Otsu, A threshold selection method from gray-level histograms, Automatica 11 (1975), 23-27. [11] C.E. Shannon, A mathematical theory of communication, ACM SIGMOBILE Mobile Computing and Communications Review 5 (2001), 3-55. [12] H. Birkholz, A unifying approach to isotropic and anisotropic total variation denoising models, Journal of computational and applied mathematics 235 (2011), 2502-2514. [13] D.L. Collins, A.P. Zijdenbos, V. Kollokian, J.G. Sled, N.J. Kabani, C.J. Holmes, and A.C. Evans, Design and construction of a realistic digital brain phantom, IEEE Trans. Med. Imaging 17 (1998), 463-468. [14] Z. Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Trans. Image Process. 13 (2004), 600-612.
Copyright of Bio-Medical Materials & Engineering is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
