MRI intensity nonuniformity correction using simultaneously spatial and gray-level histogram information

Computerized Medical Imaging and Graphics 31 (2007) 81–90

MRI intensity nonuniformity correction using simultaneously spatial and gray-level histogram information Julien Milles a , Yue Min Zhu b,d,∗ , G´erard Gimenez b , Charles R.G. Guttmann c , Isabelle E. Magnin b a

Division of Image Processing, Department of Radiology, Leiden University Medical Center, P.O. Box 9600, 2300 RC Leiden, The Netherlands b CREATIS, CNRS UMR 5515 & INSERM Unit 630, INSA Lyon, Universit´ e Claude Bernard Lyon 1, Villeurbanne, France c Center for Neurological Imaging, Brigham and Women’s Hospital, Harvard Medical School, Boston, MA, USA d School of Electrical Engineering and Automation, Harbin Institute of Technology, China Received 13 July 2005; received in revised form 31 October 2006; accepted 9 November 2006

Abstract A novel approach for correcting intensity nonuniformity in magnetic resonance imaging (MRI) is presented. This approach is based on the simultaneous use of spatial and gray-level histogram information. Spatial information about intensity nonuniformity is obtained using cubic Bspline smoothing. Gray-level histogram information of the image corrupted by intensity nonuniformity is exploited from a frequential point of view. The proposed correction method is illustrated using both physical phantom and human brain images. The results are consistent with theoretical prediction, and demonstrate a new way of dealing with intensity nonuniformity problems. They are all the more significant as the ground truth on intensity nonuniformity is unknown in clinical images. © 2006 Published by Elsevier Ltd. Keywords: Magnetic resonance imaging; Brain; Intensity nonuniformity; Correction

1. Introduction Intensity nonuniformity, also called bias field, is an artifact often observed in magnetic resonance imaging (MRI). It consists in slow nonanatomical intensity variations across images. Its presence can significantly hamper the accuracy and reliability of quantitative analysis of magnetic resonance (MR) images. Various methods have previously been proposed for correcting intensity nonuniformity. They can be loosely classified into two types: pre-processing and post-processing methods. Pre-processing correction methods mostly consist in obtaining additional information about the nonuniformity (bias) field by either imaging a simple homogeneous object [1–4] or acquiring a second image of the same object using another coil [5,6] or sequence [7]. Another approach is to focus on modeling voxel intensity with regard to the sequence used in the imaging protocol to obtain additional information on the nonuniformity artefact [8–11].

∗

Corresponding author. Tel.: +33 4 72 43 89 74. E-mail address: [email protected], [email protected](Y.M. Zhu). 0895-6111/$ – see front matter © 2006 Published by Elsevier Ltd. doi:10.1016/j.compmedimag.2006.11.001

Post-processing correction methods represent a very large majority of correction methods. Expectation-maximization algorithms were used to iteratively classify and correct the images based upon some initial probability estimates [12,13]. Local and composite medians of specific tissue classes were compared to identify bias field [14]. Lee and Vannier [15] applied the fuzzy c-means clustering algorithm to estimate bias in MR images. Other solutions have been envisioned, such as extensive averaging [16] or histogram matching [17]. The N3 (nonparametric nonuniform intensity normalization) method [18] consists of an iterative deconvolution approach followed by a B-spline smoothing to estimate the distribution of the true tissue intensities. Other approaches have modeled the gain field by polynomial methods [19–22], interpolation between userselected points [23], sophistical statistical methods [24], iterative segmentation and B-spline fitting [25]. Homomorphic unsharp masking (HUM), which behaves as a sort of band notches filters, has often been used in practice [26–31]. Other approaches work in transformed domains, such as Fourier domain [32,33] or discrete wavelet domain [34]. For an acquired image or volume, the corresponding bias field should be the same whatever the intensity nonuniformity correction method used. However, in practice, different

82

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

correction methods do not lead to the exactly same results. This has been observed by several researchers [18,35]. This inconsistency is caused by the fact that the intensity nonuniformity artefact is the result of a complex process involving the MR scanner, the acquisition sequence and the imaged object [36,37] and that we lack ground truth regarding both intensity nonuniformity and the structures of objects being imaged. In this paper, we propose a correction approach that allows collaboration between two conceptually different algorithms. It thus allows to obtain a solution consisting in the best trade-off between two different sets of assumptions, partially overcoming the lack of ground truth on intensity nonuniformity. It is a post-processing method, but in comparison with existing methods, it has a salient aspect of making simultaneous use of complementary spatial and gray-level histogram information for nonuniformity correction. Spatial information of intensity nonuniformity is obtained using cubic B-spline smoothing. Gray-level histogram information of the image corrupted by intensity nonuniformity is exploited from a frequential point of view. The proposed method will be illustrated with the aid of both physical phantom and human brain images, and also compared with the N3 method [18]. 2. Methods 2.1. Bias field models Like most of post-processing correction methods [1– 3,5,6,12,13,18–20,23,21,34], we chose to use the following multiplicative model for an observed MR image v: v(x) = g(x).u(x) + n(x)

(1)

where x is the spatial position, v the observed image, u the ideal or true image, g the bias field whose probability density function G is assumed to be a centered gaussian, and n the noise, which presents a Rician probability density function in MR images [38]. If noise is neglected, a logarithmic transform can be applied to Eq. (1) to transform the multiplicative model into an additive one. If we note aˆ (x) = log(a(x)), Eq. (1) becomes vˆ (x) = gˆ (x) + u( ˆ x)

(2)

Eq. (2) corresponds to an under-determined problem as we lack information about g and u. That is why the assumptions that g is spatially smooth, hence slowly varying, and that u is piecewise constant are widely used. The logarithmic transform does not change the nature of these hypotheses and they can directly be applied to gˆ and u. ˆ Eq. (2) can be equivalently expressed in terms of corresponding probability density functions as V (ˆv) = G(ˆv) × U(ˆv)

(3)

where V, G and U are the probability density functions of respectively vˆ , gˆ and u, ˆ under the assumption that gˆ and uˆ are independent random variables. This equation constitutes the mathematical basis of the N3 method. Eqs. (2) and (3) provide two different representations of the same bias field. The first represents it in the spatial domain,

whereas the second characterizes it in the probability density function domain. Since probability density function corresponds to gray-level histogram when considering discrete images, Eq. (3) encodes the gray-level histogram information about the bias field. Due to the complementary informations provided by those two representations, the idea is to estimate the bias field by associating Eqs. (2) and (3) in a cooperative framework according to the flowchart of Fig. 1. For a given initial volume, bias field correction consists in four parts: pre-processing of the data, spatial information exploitation, gray-level histogram information exploitation and combination of corrected volumes. The cooperative aspect relies in the use of the volume corrected using model (2) as an initial volume for model (3) and vice-versa until convergence. The essential aspects of the proposed correction method are described in the following. 2.2. Pre-processing of the data Pre-processing is composed of two steps. First we determine a region of interest (ROI) by thresholding the original data. The threshold is set to the first minimum value of the whole volumes’ gray-level histogram. This step is used to reduce volume of the data to be computed and also to avoid low intensity values, which cause numerical problems when using the logarithmic transform. The second step is the logarithmic transform. Applying this transform implies that we have previously removed regions of low intensities. Logarithmic transform is used to convert the initial multiplicative model (1) into an additive model so that models (2) and (3) can be conjointly used. Another effect of this transform is that it makes transitions between tissue compartments smoother since for any value x strictly positive, x ≥ log(x). 2.3. Spatial information exploitation This constitutes the spatial part of the correction algorithm that aims to estimate the bias field and correct the initial volume through exploiting spatial information. Starting from the hypotheses about the spatial smoothness of gˆ and u, ˆ a cubic Bspline smoothing can be used to model gˆ in a nonparametric framework. To set the smoothing parameter, we used the assumption that the ideal image being piecewise constant, its entropy should be minimal. Entropy is a measurement of information, which is computed from the probability density function of the volume using  H(ˆv) = − p(ˆvi ) log(p(ˆvi )) (4) i

where p(ˆvi ) designates the probability density of gray-level vˆ i . Once smoothing has been done, initial data are corrected using vˆ corrected = vˆ initial − vˆ smoothed

(5)

Smoothing and correction steps are iteratively alternated. The smoothing kernel extent varies for each iteration. We use two

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

83

Fig. 1. Flowchart of intensity nonuniformity correction method.

criteria to stop those iterations. The first one is an entropy-based criterion, iteration process being stopped when the entropy of the corrected volume becomes less than the one of the initial dataset (H(ˆvcorrected ) < H(ˆvinitial )). The other stopping criteron is based on the fact that when correction is complete no smooth information can be further extracted from corrected data. In this case entropy does not vary significantly between two iterations. For this reason we also check for the variation of entropy between two consecutive iterations, and stop the process if no significant variation in entropy is observed. Spatially corrected volume is then saved for further computations. 2.4. Gray-level histogram information exploitation This constitutes the statistical part of the correction algorithm. Its aim is to estimate the bias field through exploiting gray-level histogram informations. The model (3) expresses the intensity nonuniformity artefact in a convolutive framework based on probability density functions. Since for digital data probability

density functions can be approximated by gray-level histograms, model (3) shows that the histogram of the observed data is the result of convolution between the histogram of true data and that of the bias field G. Correcting for bias field then amounts to deconvolving G, the bias field’s histogram, from the histogram of the initial volume. This provides a robust approach for correcting intensity nonuniformity since only gray levels of the images are used and no spatial information is required. Moreover, deconvolution techniques available in the literature can be readily used to recover U from V. For example, making the assumption that the probability density of the bias field can be seen as a gaussian function, the probability density of the true image can then be restored by a conventional deconvolution technique. This constitutes the core of the N3 method [18]. However, as this correction model (3) only takes into account first order statistics of the processed volumes, spatial information encoded in the original dataset is lost. The deconvolution of the histogram has been processed using ˜ the Fourier transform of G a Wiener filter [39]. If we note G ∗ ˜ and G the associated complex conjugate, the Wiener filter is

84

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

given by F˜ =

˜∗ G ˜ 2 + Z2 |G|

(6)

where Z is an empirically set normalization constant. The Fourier transform of the histogram of the original volume can then be approximated by ˜ = F˜ .V˜ U

(7)

The strategy for determining the width of the gaussian kernel is the same as that used to set the distance between control points for the spatial part. The width is first set to a small initial value and increased until the entropy of the corrected volume is smaller than that of the initial volume, as done for the extraction of spatial information. Once U is estimated, in order to map U to u, ˆ a lookup table is computed using E[u|ˆ ˆ v] =

G(ˆv) × (uU( ˆ u)) ˆ G(ˆv) × U(u) ˆ

(8)

and is applied to vˆ to obtain the corrected volume [18]. 2.5. Combination of the corrected volumes The cubic B-spline smoothing correction and the gray level histogram correction model (3) may not lead to the exactly same correction results. The fact that many bias correction algorithms disagree with each other has already been observed and studied by other researchers, as mentioned before. The goal here is to look for a solution, which corresponds to the best trade-off between spatial and probability density domain corrections. Once the corrections corresponding to the spatial and statistical parts have been performed, their results are compared. A normalized mean square error (NMSE) is used as a criterion to compute the global difference between the two estimates. If this difference is smaller than a given threshold, the algorithm goes on to next fusion step. Otherwise, the two estimates being too different, the volume corrected by the spatial part is used as the initial volume for the statistical part and the volume corrected by the statistical part is used as the initial value for the spatial part. In case the convergence criterion is met, that is, the NMSE between the two corrected volumes is smaller than the threshold set by the user, the results are fused. In the present study, fusion of the results is simply done by averaging the two corrected volumes that are supposed to be similar. An exponential transform is then applied to the fused volume, yielding the final corrected volume.

modified spin echo sequences with a varying flip angle and heavily proton density-weighted (TR = 4000 ms, TE = minimum required to ensure full sampling of the k-space). The second dataset contains different volumes corresponding to T1- and T2weighted images of a human head, alongside with segmented volumes for intra-cranial tissues. The correction method was first evaluated qualitatively by visually assessing the shape and intensity variations of a socalled “concatenated histogram pattern”, which is obtained by juxtaposing slice histograms along the z-direction of the MR volume. The vertical direction of the pattern indicates the slice number. For each given slice number, the obtained line profile corresponds to the gray level histogram of that slice. With such a concatenated histogram representation, the in-slice nonuniformity is characterized by the width of the main peak of each line profile, and the slice-to-slice nonuniformity by the shape of the pattern’s ridge. In the ideal case, for a uniform object, the pattern is a vertical line with one pixel thickness and constant intensity. Second, we assessed quantitatively for the quality of the correction results by using the criteria defined in [3] for in-slice and slice-to-slice nonuniformity. Those criteria are defined as σ gin-slice (%) = 100 × (9) μ and gslice-to-slice (%) = 100 ×

σμ μμ

(10)

where μ and σ are respectively the average and standard deviation of a given tissue intensity in a slice and μμ and σμ the average and standard deviation of the averaged intensities of a given tissue throughout the volume (along the z direction). Those criteria are suitable for measuring uniform objects, we have used them to evaluate each brain tissue compartment separately. 3. Results and discussion 3.1. Results 3.1.1. Correction of homogeneous physical phantoms Each initial dataset presents a different intensity nonuniformity degree varying from 9 to 25% for in-slice nonuniformity and 8–28% for slice-to-slice nonuniformity, as shown in Table 1. Intensity nonuniformities for each dataset have been computed Table 1 In-slice and slice-to-slice nonuniformity for a homogeneous physical phantom (all numerical results are percentages) Dataset

2.6. Evaluation protocol The proposed algorithm has been tested on different physical datasets alongside with the N3 method [18] whose parameters were set to the default ones used with version 1.0 (fwhm = 0.15, stop = 0.001, iterations = 50, shrink = 4 and distance = 200). The first dataset contains 6 volumes corresponding to 60 slices acquired from a spherical phantom filled with doped water. The two-dimensional (2D) multislice acquisitions were based on

1

2

In-slice Initial data N3 corrected Proposed method

3

10.5 9.1 8.7

10.1 9.0 8.5

9.7 9.0 8.5

Slice-to-slice Initial data N3 corrected Proposed method

16.6 13.5 0.2

15.6 12.7 0.2

13.7 11.5 0.2

4

5

6

9.4 9.0 8.4

11.0 9.0 8.3

25.1 25.0 10.1

10.6 9.9 0.4

8.5 8.5 0.6

28.8 24.9 2.4

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

85

Fig. 2. Correction of homogeneous physical phantom datasets. From top to bottom: datasets 1, 3, 5 and 6.

using Eqs. (9) and (10). The more the intensity nonuniformity degree is, the greater the values of gin-slice and gslice-to-slice are. We compared our results with those obtained using the N3 method, which is solely based on histogram deconvolution, as stated in Section 2. Fig. 2 shows the central slice of initial and corrected volumes for dataset numbers 1, 3, 5 and 6. Both our method and the N3 method improve image quality for the first three displayed datasets, the last one, which presents the highest intensity nonuniformity level, being only slightly corrected by the N3 method but more significantly by the proposed method. To assess for the slice-to-slice nonuniformity correction, we used the concatenated histograms described in Section 2.6. Fig.

3 shows those histograms for dataset numbers 1, 3 and 5. Initial datasets 1 and 3 present a crescent-like pattern with a varying width, which means that both strong in-slice and slice-to-slice intensity nonuniformities are present. Dataset 5 has a bow-like pattern, this means that slice-to-slice nonuniformity evolves differently for this dataset. However, it presents almost the same in-slice and slice-to-slice nonuniformities values as the other two datasets. With the N3 method, the width of those patterns has been thinned, showing that data are corrected for in-slice nonuniformity, but the crescent- and bow-like shapes are almost unchanged, meaning that slice-to-slice nonuniformity has not been totally corrected. Our method presents the

86

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

Fig. 3. Concatenated histograms of initial and corrected homogeneous physical phantom datasets. From top to bottom: datasets 1, 3 and 5.

same decrease in in-slice nonuniformity as the N3 method, but the concatenated histogram patterns appear to be straighter, indicating that slice-to-slice nonuniformity has been better corrected. Those results are confirmed quantitatively when computing in-slice and slice-to-slice nonuniformities, as shown in Table 1. When considering in-slice nonuniformity, relative improvement in data quality provided by the N3 method ranges from 0.4% for dataset number 6 to 18.2% for dataset number 5, with an average improvement of 9%. The same measurement applied to data corrected with our method ranges from 10.6% for dataset number 4–40.2% for dataset number 6, with an average improvement of 20.1%, so the N3 method fails when applied to data corrupted by strong intensity nonuniformity but gives good results on weaker artefacts. Our method is efficient whatever the intensity nonuniformity is. This difference in behavior becomes more important when considering slice-to-slice nonuniformity. While improvement provided by the N3 method for this type of intensity nonuniformity ranges from 0% for dataset 5–18.7% for dataset 1, with an average improvement of 12.2%, our method nearly eliminates slice-to-slice nonuniformity, improvement going from 91.7 to 98.8%, with an average improvement of 96.1%. That demonstrates that combining spatial information with gray-level histogram information provides a powerful mean for correcting both types of nonuniformity.

3.1.2. Correction of brain images To assess for the efficiency of our method, we compared it with the N3 method for T1- and T2-weighted datasets. A visual comparison of the results for T1 -weighted images is shown in Fig. 4. The original image shows a slight drop in intensity when considering the rear-left part of the brain. This nonanatomical variation is detected by both methods. Estimated bias fields are globally similar. Some local differences do appear. There are more local variations in the bias field provided by our method than by the N3 method. This can also be observed on quantitative results shown in Table 2. We observe that the N3 method, contrary to the results observed on phantoms, reduces more significantly slice-to-slice nonuniformity comparing to in-slice nonuniformity. The latter shows a decrease ranging from 0.6 to 2.9%, depending on considered tissue, with a global average decrease of 1.6%, whereas slice-to-slice nonuniformity of tissues decreases from 4.9 to 8.8% with a global average decrease of 6.9%. It is worthy to underline that our method always gives substantially better results for slice-to-slice nonuniformity than for in-slice nonuniformity. The decrease in in-slice nonuniformity varies between 4.9 and 8.8% with an average decrease of 6.9%. Results for slice-to-slice nonuniformity are more constrasted, as there is a strong relative decrease, between 42.3 and 68.6% with an average of 55.7% for white and gray matter, and also for lesions, but a small increase of 6.1% for cerebrospinal

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

87

Fig. 4. Correction results for a T1 -weighted brain image. (a) Initial image, (b) difference image between the corrected images (e) and (c), (c) image corrected using the N3 method, (d) bias field estimated using the N3 method, (e) image corrected using our method, (f) bias field estimated using the proposed method. Table 2 In-slice and slice-to-slice nonuniformity for a T1 -weighted brain dataset (all numerical results are percentages) Gray matter

White matter

CSF

Lesions

In-slice Initial data N3 corrected Proposed method

10.2 10.0 8.3

14.0 13.6 11.3

20.5 20.3 18.1

17.3 17.2 14.3

Slice-to-slice Initial data N3 corrected Proposed method

13.7 12.5 7.9

22.0 20.7 6.9

16.3 15.5 17.3

19.0 17.5 8.3

fluid (CSF). This misbehavior can be explained by the fact the contrast between ventricles, filled with CSF, and surrounding tissue compartments was very high, which led to extracting wrong spatial informations when correcting data using spatial smoothing. In all cases, when comparing our method with the N3 method, the correction results obtained with the T1-weighted images are consistent with those obtained in the case of phantoms. Our method provides significant correction improvements for both in-slice and slice-to-slice nonuniformity. Results obtained on T2-weighted data confirm the trend observed on T1-weighted images. Fig. 5 illustrates correction results for a T2-weighted brain image. Visually, the correc-

88

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

Fig. 5. Correction results for a T2 -weighted brain image. (a) Initial image, (b) difference image between the corrected images (e) and (c), (c) image corrected using the N3 method, (d) bias field estimated using the N3 method, (e) image corrected using our method, (f) bias field estimated using the proposed method.

tion seems to be stronger, in view of more contrast between tissues, with our method than with the N3 method. This overcorrection can be understood when observing the estimated bias fields shown in Fig. 5. Stronger local variations of the bias field estimates occur in the same spatial location as that where the corrected images exhibit higher contrast. The N3 method reduces more slice-to-slice nonuniformity than in-slice nonuniformity, relative decrease varying between respectively 3.7 and 18.3% with an average decrease of 9.4% for in-slice nonuniformity and 8.2 and 16.1% with an average decrease of 11.0% for slice-to-slice nonuniformity as shown in Table 3. Our method

provides greater decreases varying between 12.9 and 72.9% with an average decrease of 37.9% for in-slice nonuniformity and between 15.7 and 82.8% with an average decrease of 47.8% for slice-to-slice nonuniformity. The most important nonuniformity reduction was observed for CSF, even if bias fields show stronger variations around those regions. 3.2. Discussion The proposed method appears to be more efficient for slice-to-slice nonuniformity correction than for in-slice nonuni-

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90 Table 3 In-slice and slice-to-slice nonuniformity for a T2 -weighted brain dataset (all numerical results are percentages) Gray matter

White matter

CSF

Lesions

In-slice Initial data N3 corrected Proposed method

15.8 14.4 13.4

14.7 13.7 12.8

10.7 10.3 2.9

11.5 9.4 5.7

Slice-to-slice Initial data N3 corrected Proposed method

11.5 10.2 9.7

17.0 15.6 9.4

8.7 7.3 1.5

11.7 10.7 6.1

formity correction. This can be explained by the fact that data used in the present were 2D multislice acquisitions instead of 3D isotropic acquisitions. This anisotropy in the z-direction leads to stronger and more frequent intensity variations. The method being programmed as a full 3D algorithm, no predominance is given for any direction and thus does not imply that the method is less suitable for correcting in-slice nonuniformity. The results obtained on both phantom and real brain images consistenly show that the proposed method is at least as efficient as the N3 method, whose basic assumptions are not satisfied for 2D multislice datasets. The structure of the algorithm, with two parallel branches, yields questioning regarding computation time. This concern has been integrated in the algorithm design from the start by choosing to use simple, proven and computationally efficient sub-algorithms instead of more elaborated ones. As a result, our experiments, carried out using a standard PC with a single Intel Pentium III processor, showed that computation time is of course increased but is still lower than applying successively the two sub-algorithms. This can be explained by the fact that convergence is accelerated by using the estimate provided by one of the branch for initialisation of the other branch. Moreover, in case a multi-processor computer is used, parallelisation of the code should result only in a minor, if any, increase in computation time. Due to the simple methods used as sub-algorithms, contradictions can occur and as a result the outcome may not be as good as other specifically designed methods. This may be the case for highly contrasted images, but those might not need any correction. However it has to be pointed out that our method, by alternating between two different methods, avoids the pitfalls for which the original sub-algorithms may fail. 4. Conclusion A novel approach to correction of intensity nonuniformity in MRI was proposed. The method was based on the cooperation of two algorithms exploiting respectively spatial information from cubic B-spline smoothing and the frequential information encoded in gray-level histogram of the image. This approach offers an interesting conceptual framework for intensity correction problems when the ground truth on intensity nonuniformity is unknown, which is the case for the human brain. It also allows great flexibility with respect to the assumptions made

89

on the intensity nonuniformity and the a priori information injected during correction processes. The results obtained with both phantoms and human brain images demonstrated that the proposed method allows for significant reduction of intra-slice intensity nonuniformity, and more particularly slice-to-slice intensity nonuniformity, when compared to a state of the art method for both 2D multislice and 3D isotropic data. References [1] Axel L, Constantini J, Listerud J. Intensity correction in surface-coil MR imaging. Am J Roentgenol 1987;148:418–20. [2] Condon BR, Patterson J, Wyper D, Jenkins A, Hadley DM. Image nonuniformity in magnetic resonance imaging: its magnitude and methods for its correction. Brit J Radiol 1987;60(709):83–7. [3] Wicks DAG, Barker GJ, Tofts PS. Correction of intensity nonuniformity in MR images of any orientation. Magn Reson Imag 1993;11(2):183–96. [4] Collewet G, Davenel A, Toussaint C, Akoka S. Correction of intensity nonuniformity in spin-echo T1 -weighted images. Magn Reson Imag 2002;20(4):365–73. [5] Brey WW, Narayana PA. Correction for intensity falloff in surface coil magnetic resonance imaging. Med Phys 1988;15(2):241–45. [6] Murakami JW, Hayes CE, Weinberger E. Intensity correction of phasedarray surface coil images. Magn Reson Med 1996;35(4):585–90. [7] Liney GP, Turnbull LW, Knowles AJ. A simple method for the correction of endorectal surface coil inhomogeneity in prostate imaging. J Magn Reson Imag 1998;8(4):994–7. [8] Mihara H, Iriguchi N, Ueno S. A method of RF inhomogeneity correction in MR imaging. Magn Reson Mater Phys 1998;7:115–20. [9] Thulborn KR, Boada FE, Shen GX, Christensen JD, Resse TG. Correction of B1 inhomogeneities using echo-planar imaging of water. Magn Reson Med 1998;39(3):369–75. [10] Clare S, Alecci M, Jezzard P. Compensating for B1 inhomogeneity using active transmit power modulation. Magn Reson Imag 2001;19(10):1349– 52. [11] Deichmann R, Good CD, Turner R. RF inhomogeneity compensation in structural brain imaging. Magn Reson Med 2002;47(2):398–402. [12] Wells WM, Grimson WEL, Kikinis R, Jolesz FA. Adaptative segmentation of MRI data. IEEE Trans Med Imag 1996;15(4):429–42. [13] Guillemaud R, Brady M. Estimating the bias field of MR images. IEEE Trans Med Imag 1997;16(3):238–51. [14] DeCarli C, Murphy DGM, Teichberg D, Campbell G, Sobering GS. Local histogram correction of MRI spatially dependent image pixel intensity nonuniformity. J Magn Reson Imag 1996;6(3):519–28. [15] Lee SK, Vannier MW. Post-acquisition correction of MR inhomogeneities. Magn Reson Med 1996;36(2):275–86. [16] Koivula A, Alakuijala J, Tervonen O. Image feature based automatic correction of low-frequency spatial intensity variations in MR images. Magn Reson Imag 1997;15(10):1167–75. [17] Wang L, Lai H-M, Barker GJ, Miller DH, Tofts PS. Correction for variations in MRI scanner sensitivity in brain studies with histogram matching. Magn Reson Med 1998;39(2):322–27. [18] Sled JG, Zijdenbos AP, Evans AC. A nonparametric method for automatic correction of intensity nonuniformity in MRI data. IEEE Trans Med Imag 1998;17(1):87–97. [19] Tincher M, Meyer CR, Gupta R, Williams WM. Polynomial modeling and reduction of RF body coil spatial inhomogeneity in MRI. IEEE Trans Med Imag 1993;12(2):361–65. [20] Meyer CR, Bland PH, Pipe J. Retrospective correction of intensity inhomogeneities in MRI. IEEE Trans Med Imag 1995;14(1):36–41. [21] Van Leemput K, Maes F, Vandermeulen D, Suetens P. Automated modelbased bias field correction of MR images of the brain. IEEE Trans Med Imag 1999;18(10):885–96. [22] Styner M, Brechb¨uhler C, Sz´ekely G, Gerig G. Parametric estimate of intensity inhomogeneities applied to MRI. IEEE Trans Med Imag 2000;19(3):153–65.

90

J. Milles et al. / Computerized Medical Imaging and Graphics 31 (2007) 81–90

[23] Dawant BM, Zijdenbos AP, Margolin RA. Correction of intensity variations in MR images for computer-aided tissue classification. IEEE Trans Med Imag 1993;12(4):770–81. [24] Wells WM, Kikinis R, Grimson WEL, Jolesz FA. Statistical intensity correction and segmentation of magnetic resonance image data. In: Proceedings of visualisation in biomedical computation (VBC ’94); 1994. p. 21–4. [25] Gilles S, Brady M, Declerck J, Thirion J-P, Ayache N. Bias field correction of breast MR images. In: Proceedings of visualisation in biomedical computation (VBC ’96); 1996. p. 153–58. [26] Lim KO, Pfefferbaum A. Segmentation of MR brain images into cerebrospinal fluid spaces, white and gray matter. J Comput Assist Tomogr 1989;13:588–93. [27] Harris GJ, Barta PE, Peng LW, Lee S, Brettschneider PD, Shah A, Henderer JD, Schlaepfer TE, Pearlson GD. MR volume segmentation of gray matter and white matter using manual thresholding: dependence on image brightness. Am J Neuroradiol 1994;15(2):225–30. [28] Narayana PA, Borthakur A. Effect of radio frequency inhomogeneity correction on the reproducibility of intra-cranial volumes using MR image data. Magn Reson Med 1995;33(3):396–400. [29] Johnston B, Atkins MS, Mackiewich B, Anderson M. Segmentation of multiple sclerosis lesions in intensity corrected multispectral MRI. IEEE Trans Med Imag 1996;15(2):154–69. [30] Reiss AL, Abrams MT, Singer HS, Ross JL, Denckla MB. Brain development, gender and IQ in children. A volumetric imaging study. Brain 1996;119(5):1763–74. [31] Brinkmann BH, Manduca A, Robb RA. Optimized homomorphic unsharp masking for MR grayscale inhomogeneity correction. IEEE Trans Med Imag 1998;17(2):161–71. [32] Wald LL, Carvajal L, Moyher SE, Nelson SJ, Grant PE, Barkovich AJ, Vigneron DB. Phased array detectors and an automated intensity-correction algorithm for high-resolution MR imaging of the human brain. Magn Reson Med 1995;34(3):433–9. [33] Cohen MS, DuBois RM, Zeineh MM. Rapid and effective correction of RF inhomogeneity for high field magnetic resonance imaging. Hum Brain Mapp 2000;10(4):204–11. [34] Han C, Hatsukami TS, Yuan C. A multi-scale method for automatic correction of intensity nonuniformity in MR images. J Magn Reson Im 2001;13(3):428–36. [35] Arnold JB, Liow J-S, Schaper KA, Stern JJ, Sled JG, Shattuck DW, Worth AJ, Cohen MS, Leahy RM, Mazziotta JC, Rottenberg DA. Qualitative and quantitative evaluation of six algorithms for correcting intensity nonuniformity effects. NeuroImage 2001;13(5):931–43. [36] Alecci M, Collins CM, Smith MB, Jezzard P. Radio frequency magnetic field mapping of a 3 tesla birdcage coil: experimental and theoretical dependence on sample properties. Magn Reson Med 2001;46(2):379–85. [37] Milles J, Zhu YM, Chen N, Panych L, Gimenez G, Guttmann CRG. Computation of transmitted and received b1 fields in magnetic resonance imaging. In: Proceedings of computer assisted radiology and surgery (CARS ’03); 2003. p. 1323. [38] Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson Med 1995;34(6):910–14. [39] Pratt WK. Digital image processing. 3rd ed.. New York (USA): John Wiley & Sons, Inc.; 2001. p. 722.

Julien Milles received the MSc degree in electrical engineering from the Institut National Polytechnique of Grenoble, France, in 1999 and the PhD degree in image processing from the Institut National des Sciences Appliquées of Lyon, France, in 2002. From 2002 to 2004, he was a Research Fellow at the Department

of Biophysics of the Maastricht University, The Netherlands. Since 2004, he has worked as a researcher in the Division of Image Processing, Department of Radiology, of the Leiden University Medical Center, The Netherlands. His research interests are in magnetic resonance images processing for both neurological and cardiological purposes. Yuemin Zhu received a MS degree (1984), a PhD degree (1988) and the “Habilitation a` Diriger des Recherches” (1993), all from the INSA (Institut National des Sciences Appliquées) of Lyon, France. He is a permanent researcher (Directeur de recherche) of the CNRS (Centre National de la Recherche Scientifique) of France. He has conducted over 20 research projects including those related to European Union, in the field of industrial nondestructive evaluation and testing, and medical imaging. His research interests include medical MRI (modeling, reconstruction, bias correction, segmentation, visualization), data fusion, timefrequency analysis, multiresolution analysis, and local spectral analysis. He has published over 120 articles. He serves as China Affairs Manager in the INSA Lyon. He is also Honorary Professor of Chongqing Three Gorges University and Harbin Institute of Technology (HIT). G´erard Gimenez received his MSc degree from the University of Lyon in 1970, his Doctorat from the University of Marseille in 1972, and his PhD degree from the University of Paris-Sud Orsay in 1978. He has been with CREATIS, a joint laboratory of the Centre National de la Recherche Scientifique (CNRS), the Institut National des Sciences Appliqu´ees de Lyon (INSA Lyon), and the University Claude Bernard Lyon 1 since 1980. His research interests are in the field of image and signal processing, particularly involving interaction between waves and biological tissues. Involved in the management of CREATIS since 1992, he was the director of the lab between 1995 and 2001. Since 2001 he is director of “CERMEP, imagerie du vivant”, an imaging center in Lyon currently evolving toward a multi-modalities imaging facility for Man and Animal (PET, ␮-PET, MRI, ␮-MRI, MEG/EEG, ␮-X-rays scanner, ␮-echographs). He is currently a professor with the Electrical Engineering Department at INSA Lyon. Charles R.G. Guttmann is the Director of the Center for Neurological Imaging at Brigham and Women’s Hospital and an Assistant Professor in Radiology at Harvard Medical School. His main interest is the quantitative evaluation of normal and pathological states of the brain using MRI. Among Dr. Guttmann’s specific goals are the understanding of the natural course of MS and of white matter disorders in the elderly. Using MRI findings as phenotypic descriptors and elucidating the relationship between brain morphological changes and functional deficits is particularly emphasized. He graduated from the University of Zurich with a Doctorate in Medicine and completed the Post-Graduate Course in Experimental Medicine and Biology at the same Institution. His MD thesis on the invasivity of human glioma cell lines in the central nervous system was done under the direction of Dr. Martin E. Schwab at the Brain Research Institute of the University of Zurich. Following his Post-Doctoral training under Dr. Ferenc A. Jolesz at Brigham and Women’s Hospital (Harvard Medical School) he went on to found the CNI in 2000. Isabelle Magnin is the Head of CREATIS (Centre of Research and Applications in Signal and Image Processing) in Lyon, France. At a French level, she is the Director of CNRS-INSERM Research group in Health and Information Technology. She is an expert for the European commission IST e-Health and Honorary Professor of Harbin University in China. Her research interest mainly focuses on Medical Image Processing including tomography, motion estimation, analysis and modeling of 2D and 3D deformable objects from multimodality medical image sequences (CT, MRI, US, PET). Since 2000, she promoted biomedical image processing applications (Cardiac imaging and MRI simulation) in European Healthgrid Projects. In 1994, she got the Gold medal for innovation from the European Society of Non Destructive Testing for original work in ultrasonic image processing. She is the author of about 230 publications.