Efficient MR Inhomogeneity Correction by Regularized Entropy Minimization and Proximal Alternations Bo Zhang1 and Hans Peeters2 1


Medisys, Philips Research, France MRI Development, Philips, Netherlands

ABSTRACT Magnetic Resonance (MR) images usually exhibit intensity inhomogeneity (bias field) due to non-uniformity generated from RF coils or radiation-patient interactions. This inhomogeneity could lead to undesired perception and difficult diagnosis. In this work, we propose a nonparametric retrospective bias corrector by minimizing a regularized-entropy criterion through a proximal alternating scheme. Numerically, the proposed solver is particularly efficient, scalable and parallelizable compared to existing approaches. Additionally, the effectiveness and the robustness of the corrector have been validated by simulations and vast clinical tests. Index Terms—inhomogeneity correction, bias entropy minimization, proximal iterative solver 1.



Magnetic Resonance (MR) image intensity inhomogeneity (aka bias field, shading) arises from imperfections of image acquisition process, and manifests itself as slowly varying intensities across the image. Many factors may give rise to the inhomogeneity artifact, such as magnetic field nonuniformity and radiation-patient interactions [1]. This shading phenomenon reduces the general image quality and is unfavorable for image perception and diagnosis. Both prospective and retrospective bias-field reduction has been proposed in the literature (see [1]). Prospective correction aims at system calibration, usually requiring a specific protocol prior to patient scans. Retrospective approaches have more flexibility of being data-driven, such that the bias-field estimate mainly relies on the information of scanned data. In this paper, we propose a retrospective bias corrector, which is based on the minimization of a regularized Shannon-entropy functional. As a matter of fact, the slowly varying bias shading acts as a multiplicative artifact on the image of tissues ideally modeled to be piecewise constant. On the image histogram, the shading introduces a blurring nature: the histogram originally containing a limited number of peaks (indicating different tissue classes) tends to be smoothed away, with those peaks spread out as the bias modulates the intensity inside the tissue regions. From an

978-1-4799-2374-8/15/$31.00 ©2015 IEEE

information-theoretical perspective, this process increases the image entropy. This rationale underlies both our method, and existing works such as [2] [3] [4] [5] which exploit entropy criterion minimization. It is common to have an explicit parametric bias modeling (e.g. polynomials) which embodies some shading regularity assumptions [3] [4] [5]. The regularity parameters, e.g. the order of polynomials, are closely related to the degrees of freedom of the model, and therefore, they strongly impact the computational complexity of the corrector. However in reality, MR users are facing tens or even more MR acquisition protocols each being able to produce a different inhomogeneity field. With the same bias corrector, a varying delay is undesirable, as it increases the workflow complexity while reducing the predictability of the system performance. As a consequence, we opt for a nonparametric bias model, in which the shading regularity is penalized through a fractional Laplacian operator. In order to further improve the computational efficiency, we propose to approximate the entropy by the image likelihood, and to apply a proximal-relaxation based schema that alternatively minimizes the entropy and the regularization penalty. It turns out that our numerical solver only involves parallelizable pointwise operations in the image domain and in a transform domain (in our case the discrete cosine transform or DCT). This reveals to improve significantly the shading-correction performance by reducing the processing time from order of minutes [2] [3] [5] to seconds. Finally, the effectiveness and the robustness of our corrector have been benchmarked on simulated data, and validated by vast clinical tests. 2. METHOD 2.1 Problem formulation Our bias-field to estimate is modeled as a point-wise multiplicative field in an image volume: (1) = ⋅ , where is the observed image contaminated by the bias, the unobserved bias-free image, and = ; ; the location of any voxel. By applying the logarithmic scale on both sides, (1) is transformed into an additive model:



+ = log


(2) ,

= log Above, we have set , and = log . In order to estimate the bias , we minimize the following regularized-entropy functional: (3) argmin ≔ℰ + + . The term ℰ denotes the Shannon entropy of the corrected image; stands for a regularization term imposing some smoothness penalty on the bias field; represents an indicating function of some constraints ! that the bias field is assumed to satisfy a priori, i.e., ∶= #

0 +∞

∈! . '(ℎ*+,-.*

Details of these terms are given in the following sections. 2.1.1 The entropy Let the probability of observing an intensity level / in the − corrected image to be 01 ≔ Pr − = / ≈ #6 : = /8/:, with # denoting a set cardinality and : the total number of image voxels. Then Shannon entropy, by definition, is expressed in the domain of the image-intensity range ℒ, i.e., ℰ ≔ − ∑1∈ℒ 01 log 01 . It turns out to be especially useful to approximate the entropy by the image log-likelihood. This allows us to deal with the entropy in the image space, in a pointwise manner: ℰ =−

≈ −=



#6 :

1 = log Pr? : @E


− :

= /8



log 01 BCD

1 = log Pr : @E




where we have set Pr @ ∶= Pr? @ − @ A. The leading factor 1/: in (4) can be deemed as a confidence weighting attributing to each voxel. More generally, one may introduce a weighting function F (see (5)) encoding some customized confidence levels. A typical example is the binary-weighting case, where we associate zero to voxels belong to the background, and one for those of object of interest. ℰ

= −= F @


log Pr



2.1.2 Regularization Regarding the smoothness penalty of the bias-field, we K

J J choose ≔ GHI G , where ‖∙‖ is the NK -norm, and HI is a (negative) fractional Laplacian operator such that: J


∶= O−PKQ





Parameter P = VPQ ; PS ; PT W balances the amount of the regularization imposed on the bias field compared to the other terms in (3); parameter X > 0 controls the degree of the smoothness. The negative signs in (6) ensure the semipositive-definiteness of the operator. Operator HIJ should be understood in the distribution sense. In finite-dimensional discrete space, HIJ represents a

matrix defined in its discrete spectral domain. As a matter of fact, HIJ is diagonalizable by an orthogonal basis Z (Z ∗ Z = \ with Z ∗ being the conjugate transpose of Z): J (7) HI = Z ∗ ]J Z, ] = −PKQ ]QQ − PKS ]SS − PKT ]TT ]QQ , ]SS and ]TT are diagonal matrices encoding, on their main diagonals, eigen-values of some 2nd-order finitedifference matrices on x, y and z directions, respectively. By adopting an even-symmetric image boundary assumption, our basis Z is given by the DCT matrix [6]. This fact is of significant importance for designing computationally efficient numerical schemes. 2.1.3 Constraints We assume our bias field to be lower and upper bounded and to be of zero average, i.e., ∈ ! where ! ≔ 6 |


_bQ , = Q

= 08.


The zero-average condition preserves the global average of the uncorrected image . ! being a convex set, we can define 0+'c as the projection of a bias field onto the constraint set !. 2.2 Numerical scheme For solving (3), we propose to apply the Alternating Direction Method of Multipliers (ADMM) [7] [8], which employs a divide-and-conquer strategy. We summarize the scheme in (9)-(11) where de > 0 is non-decreasing (f = 0,1,2, …): eiD



= arg min jGHI


G +

de ‖ − ke − le ‖K m 2

CD 2 J∗ J HI HI + \o ke + le de de = arg min qℰ k + k + ‖ 2 p de = arg min = V eiD @ − k 2 p∈



− F @ log Pr @ = le − eiD + keiD

(9) eiD @

− k − le ‖K r

− le




(10) (11) Briefly speaking, the algorithm splits (3) into two groups consisting of the regularity term and of the constrained entropy ℰ + , respectively. Both groups are relaxed using proximal operators before being minimized alternatively (i.e., (9) and (10)). Following the primal minimization, the dual variables are updated in (11). The entire process is iterated till convergence with the last obtained e being our solution. More specifically, the minimizer to (9) can be efficiently solved in the transform domain, thanks to the decomposition in (7), i.e., leiD


= Z∗ n

CD 2 KJ ] + \o Z ke + le de


The applications of Z and Z ∗ in (12) are DCT and inverse DCT, respectively. The diagonal inverse matrix in (12)


encodes the DCT spectrum, which can be pre-computed once for all. The solution to (10) is approximated by a step of projected gradient-descent, i.e.: ∗ (13) keiD = 0+'c VseiD − .eiD ∇u seiD W, K de u k = = V eiD @ − k @ − le @ W 2 @


(14) @ log Pr @ , ∗ Ru Ru Ru (15) ∇u k = j , ,…, m, Rk Rk D Rk BCD Ru = de Vk @ + le @ − eiD @ W Rk @ Pr′V @ −k @ W + F @ . PrV (16) @ −k @ W In (16), the intensity probability and its derivative are computed from the (smoothed) histogram of − k . The gradient descent in (13) is computed independently on each summing term of (14), which boils down to 1D descent problem on every voxel. Every component of the steplength vector .eiD may also be computed element-wise (through, e.g., line searching), such that u k does not increase after the gradient descent. The result of the gradient descent is then projected onto the constraint set !. In words, the algorithm at each iteration alternates between (i) a filtering through a point-wise DCT spectral modulation, and (ii) a voxel-wise projected gradientdescent. Both steps are highly scalable and parallelizable as well as the DCT. 3. RESULTS 3.1 Results on simulated brains To quantify the bias reduction, we benchmarked our shading corrector with a histogram-deconvolution method [9] (N3) and a classification-based method [10] (FC) on healthy brain phantoms of McGill's BrainWeb1. Acquisition protocols T1, T2 and PD were involved. Bias fields were measured from separate scans, and were linearly scaled to obtain 20%, 40% and 60% intensity variations. Three levels of Rician noise were added on each image to three SNR levels. For each SNR, ten replications were carried out for generating different noise realizations. In all, we built a database of 270 simulated volumes each of size 181 × 217 × 181 and of resolution 1zz{ . Three tissue areas viz white matter (WM), gray matter (GM) and cerebrospinal fluid (CSF) were segmented and provided by BrainWeb. The intensity uniformity of a given tissue area | was quantified by the coefficient of variation (i.e. ratio between the standard deviation and the average of intensity): (17) !}~ ≔ •~ /d ~ . 1

Our bias fields were slowly varying which showed a low-pass nature. We found that the images could be subsampled to a much coarser resolution (6zz{ ) without degrading the bias estimate. Further, we used a binary function F (see (5)) for foreground and background separation. The background intensity was detected as the first local minimum of the image histogram. We also used no more than 15 iterations, with PQ = PS = PT = 4, X = 1 and de = 2. The reduction of CV was evaluated for all three shading correctors. Regarding FC in particular, we also included the hypotheses of 3 (FC-3), 4 (FC-4) and 5 (FC-5) tissue classes. The ranking results are shown in TABLE I. It can be seen that our approach outperforms the other correctors in most of the cases. TABLE I. Reduction of CV of WM, GM and CSF according to different shading correctors. Each table cell shows (i) the percentage of a given method ranked first for a given criterion; (ii) in the parenthesis, the average rank of the method for that criterion. The best results among the different correctors are shown in bold. N3 FC-3 FC-4 FC-5 Crit. Proposed !}‚ƒ 59.3% (1.6) 29.6% (2.0) 7.4% (4.0) 3.7% (3.7) 0.0% (3.7) !}„ƒ 63.0% (1.6) 25.9% (2.1) 7.4% (3.8) 3.7% (3.7) 0.0% (3.7) !} …† 55.6% (1.7) 25.9% (2.3) 11.1% (3.6) 3.7% (3.6) 3.7% (3.8)

3.2 Clinical validations Our methods were clinically tested on a large range of breast and brain datasets, where radiologists were asked to blindly assess the quality of images before and after correction presented side by side. For breast application, 21 patients containing 43 volumes of T1, T2 and DCE scans (both 1.5T and 3T) were evaluated. All corrected images were ranked higher than original ones, with no loss of diagnostic information. Regarding the brain application, 50 patients containing 242 volumes of T1, post-contrast T1, T2 and FLAIR scans (both 1.5T and 3T) were evaluated. About 90% corrections were preferred over uncorrected images, with not a single report on diagnostic quality loss. Fig. 1 shows several cases of uniformity improvement due to the correction. In the brain examples, the bias-field is compensated either in the anterior-posterior direction or a in the central area. In the breast cases, unbalanced intensities mainly on the left and right breasts are corrected. This shows the effectiveness and the robustness of the algorithm in real applications. 3.3 Results on computational load The algorithm converges fast in practice: in all our tests, we used no more than 15 iterations which seemed enough to produce satisfactory corrections. This requires around 100 evaluations of the entropy term (14) per image, about an order of magnitude less than that reported in [4].



80% computation parallelized are shown in blue. Green curve shows our measurements using up to 6 cores.

4. CONCLUSION  In this paper, we have proposed an entropy-based nonparametric retrospective bias-field corrector with a very efficient numerical implementation. The approach seems promising on various anatomical applications, and has shown excellent robustness on brain and breast clinical tests.



The authors would like to thank Ad Moerland, Niccolo Stefani, Yuxi Pang and Trevor Andrew of Philips for their organization of the clinical evaluations of the shading corrector. 6.  Fig. 1 Inhomogeneity corrections: blue arrows show the transition from original to corrected images; green markers highlight the inhomogeneous areas of interest. From left to right and top to bottom: brain-T1, brain-FLAIR, breast-T2, and breast-T1 cases.

Further, to estimate the time performance and the degree of parallelism of our algorithm (C++ implemented), we tested it using a breast image of size 512 × 512 × 250 voxels. The voxel resolution was 0.78 × 0.78 × 2.4 zz{ which was subsampled to 6 × 6 × 2.4 zz{ for bias estimation. Fig. 2 shows the computational time as a function of number of threads used. Even on a single core, the correction lasts less than 10 seconds, in contrast to order of minutes reported in [2] [3] [5] on images of comparable or even smaller sizes. In the same figure, we show the time accelerating factors (ratio of single-thread computational time over multi-thread time) as a function of the number of threads. These statistics confirm the high parallelizability of the algorithm which is estimated to be at least 80%. 4

Accelerating factor

# Threads Comput. time 1 8.34 secs 2 4.83 secs 3 3.76 secs 4 3.12 secs 5 2.80 secs 6 2.65 secs

3.5 3 2.5 Theoretical accelerating factors


Measured accelerating factors

1.5 1



6 8 10 12 Number of threads



Fig. 2. Left: Execution time on a breast image (512 × 512 × 250 voxels, resolution: 0.78 × 0.78 × 2.4 zz{ subsampled to 6 × 6 × 2.4 zz{ for bias estimation) with 15 iterations at most on a 6-core 3.33 GHz Intel® Xeon CPU. Right: Acceleration factors as a function of number of threads. Theoretical factors assuming


[1] U. Vovk, F. Pernus and B. Likar, "A Review of Methods for Correction of Intensity Inhomogeneity in MRI," IEEE T Med Imaging, vol. 26, no. 3, pp. 405-421, 2007. [2] J.-F. Mangin, "Entropy minimization for automatic correction of intensity nonuniformity," in Math Method Biomed Image Anal, Hilton Head Island, SC, 2000. [3] B. Likar, M. A. Viergever and F. Pernus, "Retrospective Correction of MR Intensity Inhomogeneity by Information Minimization," IEEE T Med Imaging, vol. 20, no. 12, pp. 1398-1410, 2001. [4] Q. Ji, J. O. Glass and W. E. Reddicka, "A novel, fast entropyminimization algorithm for bias field correction in MR images," Magn Reson Imaging, vol. 25, no. 2, pp. 259-264, 2007. [5] J. V. Manjon, J. J. Lull, J. Carbonell-Caballero, G. GarciaMarti, L. Marti-Bonmati and M. Robles, "A nonparametric MRI inhomogeneity correction method," Med Image Anal, vol. 11, no. 4, pp. 336-345, 2007. [6] G. Strang, "The discrete cosine transform," SIAM review, vol. 41, no. 1, pp. 135-147, 1999. [7] J. Eckstein and D. Bertsekas, "On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators," Math Program, vol. 55, pp. 293-318, 1992. [8] M. Afonso, J. Bioucas-Dias and M. Figueiredo, "Fast Image Recovery Using Variable Splitting and Constrained Optimization," IEEE T. Image Process., vol. 19, no. 9, pp. 2345-2356, 2010. [9] J. G. Sled, A. P. Zijdenbos and A. C. Evans, "A nonparametric method for automatic correction of intensity nonuniformity in MRI data," IEEE T Med Imaging, vol. 17, no. 1, pp. 87-97, 1998. [10] D. L. Pham and J. L. Prince, "An adaptive fuzzy C-means algorithm for image segmentation in the presence of intensity inhomogeneities," Pattern Recognit. Lett., vol. 20, no. 1, pp. 57-68, 1999.


Efficient MR Inhomogeneity Correction by Regularized ... - IEEE Xplore

Magnetic Resonance (MR) images usually exhibit intensity inhomogeneity (bias field) due to non-uniformity generated from RF coils or radiation-patient ...

624KB Sizes 2 Downloads 190 Views

Recommend Documents

Transparent Error Correction for Communication ... - IEEE Xplore
Jun 15, 2011 - TCP/IP throughput by an order of magnitude on a 1-Gb/s link with 50-ms ... protocols, aggregating traffic for high-speed encoding and using a.

An Efficient MR Inhomogeneity Corrector Using ...
flexibility of being data-driven compared to prospective correction which usually requires a specific protocol ... Fast Image Recovery Using Variable Splitting and.

Energy-Efficient Opportunistic Interference Alignment - IEEE Xplore
Abstract—We introduce an energy-efficient distributed op- portunistic interference alignment (OIA) scheme that greatly improves the sum-rates in multiple-cell ...

Computationally Efficient Template-Based Face ... - IEEE Xplore
head poses, illuminations, ages and facial expressions. Template images could come from still images or video frames. Therefore, measuring the similarity ...

Efficient Multiple Hypothesis Tracking by Track Segment ... - IEEE Xplore
Burlington, MA, USA. {chee.chong, greg.castanon, nathan.cooprider, shozo.mori balasubramaniam.ravichandran}@baesystems.com. Robert Macior. Air Force ...

IEEE Photonics Technology - IEEE Xplore
Abstract—Due to the high beam divergence of standard laser diodes (LDs), these are not suitable for wavelength-selective feed- back without extra optical ...

Spatial-Temporal Enhancement Pattern, for MR-based ... - IEEE Xplore
features of tumor, but also the spatial variations of pixel-wise temporal enhancement of tumor. Although the latter has been widely used by radiologist during ...

Energy Efficient Content Distribution in an ISP Network - IEEE Xplore
The content data is delivered towards the clients following a path on the tree from the root, i.e., the Internet peering point. A storage cache can be located at each node of the network, providing a potential facility for storing data. Moreover, cac

Bandwidth-Efficient WDM Channel Allocation for Four ... - IEEE Xplore
52, NO. 12, DECEMBER 2004. Bandwidth-Efficient WDM Channel Allocation for. Four-Wave Mixing-Effect Minimization. Vrizlynn L. L. Thing, P. Shum, and M. K. ...

An Energy Efficient Multi-channel MAC Protocol for ... - IEEE Xplore
Department of Computer Engineering, Kyung Hee University, 449-701, ... Department of Electronics and Communications Engineering, Kwangwoon University, ...

A Gram-Based String Paradigm for Efficient Video ... - IEEE Xplore
Mar 13, 2013 - semination of video data has created an urgent demand for the large-scale ... video sequence, retrieval is performed by transforming the query.

An Efficient Geometric Algorithm to Compute Time ... - IEEE Xplore
An Efficient Geometric Algorithm to Compute Time-optimal trajectories for a Car-like Robot. Huifang Wang, Yangzhou Chen and Philippe Sou`eres.

Efficient Estimation of Critical Load Levels Using ... - IEEE Xplore
4, NOVEMBER 2011. Efficient Estimation of Critical Load Levels. Using Variable Substitution Method. Rui Bo, Member, IEEE, and Fangxing Li, Senior Member, ...

Efficient Pattern Matching Algorithm for Memory ... - IEEE Xplore
intrusion detection system must have a memory-efficient pat- tern-matching algorithm and hardware design. In this paper, we propose a memory-efficient ...

wright layout - IEEE Xplore
tive specifications for voice over asynchronous transfer mode (VoATM) [2], voice over IP. (VoIP), and voice over frame relay (VoFR) [3]. Much has been written ...

Device Ensembles - IEEE Xplore
Dec 2, 2004 - time, the computer and consumer electronics indus- tries are defining ... tered on data synchronization between desktops and personal digital ...

wright layout - IEEE Xplore
ACCEPTED FROM OPEN CALL. INTRODUCTION. Two trends motivate this article: first, the growth of telecommunications industry interest in the implementation ...

A Peak Power Efficient Cooperative Diversity Using Star ... - IEEE Xplore
Abstract—In this paper, we propose a new simple relaying strategy with bit-interleaved convolutionally coded star quadra- ture amplitude modulation (QAM).

Tree Filtering: Efficient Structure-Preserving Smoothing ... - IEEE Xplore
Dec 17, 2013 - Linchao Bao, Yibing Song, Qingxiong Yang, Member, IEEE, Hao Yuan, and Gang ... joint filtering, minimum spanning tree, structure-preserving.

Efficient Computation of Transfer Function Dominant ... - IEEE Xplore
Jan 12, 2006 - and produce good modal equivalents automatically, without any ... modal equivalents, model reduction, poorly-damped oscillations,.

Efficient Computation of Transfer Function Dominant ... - IEEE Xplore
Jan 12, 2006 - dominant poles of a high-order scalar transfer function. The al- gorithm ... power system dynamics, small-signal stability, sparse eigenanal-.

Study and Correction of Field Inhomogeneity in ...
5.4 Linear mapping from data to displayable values. . . . . . . . . . . 54 ...... soft pulse (because smoother) while a non-selective pulse is called hard pulse. 13This method is not ..... and recover I(x) ? Not yet. ...... 5.12, we see that the line

Evolutionary Computation, IEEE Transactions on - IEEE Xplore
search strategy to a great number of habitats and prey distributions. We propose to synthesize a similar search strategy for the massively multimodal problems of ...