2

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.

field,

INTRODUCTION

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:

1069

=

+ = 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: ℰ =−

≈ −=

1∈ℒ

BCD

#6 :

1 = log Pr? : @E

@

− :

−

= /8

@

A=−

log 01 BCD

1 = log Pr : @E

@

,

(4)

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

@

(5)

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

HI

∶= O−PKQ

RK RK RK K K − P − P U S T RTK RQK RSK

J

.

(6)

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.

(8)

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

keiD

J

= arg min jGHI

K

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∈

=n

@

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

(9) eiD @

− k − le ‖K r

− le

@

W

K

(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

eiD

= Z∗ n

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

(12)

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

1070

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 @

−F

(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].

http://www.brainweb.bic.mni.mcgill.ca/brainweb/

1071

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.

5.

ACKNOWLEDGMENT

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

2

Measured accelerating factors

1.5 1

2

4

6 8 10 12 Number of threads

14

16

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

REFERENCES

[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.

1072