A Review on MR Image Intensity Inhomogeneity Correction

Intensity inhomogeneity (IIH) is often encountered in MR imaging, and a number of techniques have been devised to correct this artifact. This paper attempts to review some of the recent developments in the mathematical modeling of IIH field. Low-frequency models are widely used, but they tend to corrupt the low-frequency components of the tissue. Hypersurface models and statistical models can be adaptive to the image and generally more stable, but they are also generally more complex and consume more computer memory and CPU time. They are often formulated together with image segmentation within one framework and the overall performance is highly dependent on the segmentation process. Beside these three popular models, this paper also summarizes other techniques based on different principles. In addition, the issue of quantitative evaluation and comparative study are discussed.


INTRODUCTION
With the frequent application of the magnetic resonance (MR) imaging method to clinical diagnosis, automatic analysis of the acquired images using techniques from computer vision and pattern recognition has received considerable attention. In developing such computer-aided diagnosis tools, a commonly encountered problem is to correct the intensity inhomogeneity (IIH) in MR images.
The IIH (also termed as the intensity nonuniformity, the bias field, or the gain field in the literature) usually refers to the slow, nonanatomic intensity variations of the same tissue over the image domain. It can be due to imaging instrumentation (such as radio-frequency nonuniformity, static field inhomogeneity, etc.) or the patient movement [1][2][3][4][5]. This artifact is particularly severe in MR images captured by surface coils. Two real MR images with severe IIH artifact are shown in Figure 1(a), where one can see that the intensity varies significantly for the pixels of the same tissue and the intensity values overlap markedly between the pixels of the different tissues. For comparison, the IIH corrected images by a surface fitting technique [6] are given in Figure 1(b), from which the improvement in image quality is clearly visible. The estimated IIH maps are given in Figure 1(c).
Let x denote the measured intensity and x the true intensity. Then the most popular model in describing the IIH effect is where α denotes the IIH effect and ξ the noise. Notation of bold letters refers to 2D or 3D MR data. Figure 2 displays the widely used BrainWeb [7] simulated MR images, where on To simplify the computation, one often ignores the noise and takes the logarithmic transform of intensity where x i is the intensity at voxel i (i = 1, . . . , n). Here, to avoid numerical problems, care should be taken for those pixels/voxels with low intensities, which are usually excluded from computation. In general, the presence of IIH can significantly reduce the accuracy of image segmentation and registration, hence decreasing the reliability of subsequent quantitative measurement. A number of techniques have been proposed to deal with this issue. In general, if a map of the IIH in the image domain (Figure 1(c) for instance) is known or can be estimated, then it is simple to correct the IIH by division in (1) or subtraction in the log-domain (2). One can obtain the IIH map from measurement in vivo [8][9][10][11][12][13][14][15], typically of a uniform phantom, [16][17][18][19][20], which often requires extra measurement (and increases the scanning time) or needs additional hardware which may not be readily available in some clinical departments. Also there are theoretical modeling approaches 2 International Journal of Biomedical Imaging  [ [21][22][23][24][25][26][27][28] to approximate the IIH map. However, due to the complexity that causes the IIH, it is difficult to model the IIH under a variety of imaging conditions. In particular, the object-induced IIH is hard to be accounted for by phantom study or theoretical modeling.
More often, the IIH map is derived retrospectively from the image data alone. A number of research efforts have been put in this direction and many techniques have been proposed. Popular mathematical models for IIH description can be classified as follows: (1) low-frequency model, which assumes the IIH to constitute low-frequency components in frequency domain and the IIH map can be recovered by lowpass filtering; (2) hypersurface model, which fits the IIH map by a smooth functional, whose parameters are usually obtained using regression; (3) statistical model, which assumes the IIH to be a random variable or a random process and the IIH map can be derived through statistical estimation; (4) others, which are based on different principles, and sometimes without explicit assumptions on the IIH field.
With this in mind, the IIH correction methods are categorized into lowpass filtering, statistical modeling, surface fitting and others, which are detailed, respectively, in the following sections. For an early literature review, interested readers are referred to [29], where an evaluation of the IIH correction effect for brain tumor segmentation is also reported. This paper attempts to summarize the recent progress and focus will be on mathematical modeling for IIH removal. Nevertheless, it is by no means an exhaustive summary. For simplicity, the description will be on single-channel data only.

LOWPASS FILTERING
Since the IIH is slowly varying in the image domain, its spectrum in frequency domain will be concentrated in the lowfrequency end. Therefore the IIH could be separated from the true image by a lowpass filter, L. After lowpass filtering in log-domain, one would approximately have This procedure to correct the IIH is similar to the homomorphic filtering in digital image processing for the correction of illumination inhomogeneity [30]. In fact, (1) easily reminds Zujun Hou 3 one of the illumination-reflectance model in optical imaging [30], where the artifact from the illumination inhomogeneity is often termed "shading" in the literature. Thus, techniques for shading correction, such as the homomorphic filtering, can be adopted for IIH removal and the converse also holds. An investigation of applying IIH correction methods to deal with the shading problem in microscopic images has been carried out in [31]. Due to their simplicity and efficiency in implementation, lowpass filtering methods have been widely used [30][31][32][33][34][35][36][37][38][39][40][41][42][43]. For a summary, interested readers are referred to [42]. Also in [42], the impact of filter width on IIH correction was investigated and it was found that these methods should be used with care to avoid intensity distortion and artificial artifacts in the corrected images. Basically, for MR images, due to the overlapping spectrum between the patient data and the IIH, the effectiveness of most conventional lowpass filtering in removing the IIH is generally quite limited.
Luo et al. [44] presented a technique to recover lowfrequency components which correspond to anatomic structure and are lost during the lowpass filtering. The method expresses the signal with a linear combination of singularity functions. The higher-frequency components are assumed to be less affected by the IIH and are used to reconstruct the true image, after which the ratio between the observed and estimated image is used for IIH map approximation.
Recently, lowpass filtering methods have been extended using the wavelet transform [45,46] and were shown to be effective in removing IIH in images acquired by surface coils and phase array coils. Compared with usual lowpass filtering methods, the multiresolution analysis allows one to select an optimal scale from which the approximate band in the wavelet transform domain is used for estimating the IIH map.
In [47], an improvement of a lowpass filtering method [43] was presented. The method varies the filter kernel size to minimize the segmentation error. The idea is generally similar to [45,46] in addressing IIH correction from the scale space, but differs in the criterion to determine the optimal scale.

SURFACE FITTING
Since the inhomogeneity field is slowly varying, it is natural to approximate the IIH by a parametric smooth functional [6,[48][49][50][51][52][53][54][55][56][57][58]. Very often, the parameter estimation is linked to image segmentation. In this way, the two different problems, IIH correction and image segmentation, are formulated in one framework and solved simultaneously. Alternatively, the parameter searching can be guided through the variation of some global image feature in an iterative process. A typical example is to minimize the entropy of gray-level histogram.

Segmentation
A large category of surface fitting approaches search the parameters by fitting with respect to a set of tissue points encoding information about the IIH. Let I = {1, . . . , n} in-dex the voxel coordinates of brain tissue. Then, in order to determine the parameters of this functional, one needs to find/segment a set of voxels S I ⊆ I which convey information about the IIH map. Among these methods, the essential difference lies in the identification of S I , and hence decoding the IIH information from S I .
Dawant et al. [48] proposed to manually select S I such that they belong to the same type of tissue. As a result, the intensity variation among these voxels can largely be attributed to IIH. However, expert supervision to select points is time consuming and error prone, especially for volume data. Wang et al. [59] has presented an automated method for generating the reference points.
Meyer et al. [49] employed the LCJ method [60] to preliminarily segment the image and then fit a smooth functional over the segmented image. The LCJ segmentation method assumes the image to be piecewise smooth and requires that the different objects are well separated at the boundaries, which is quite stringent in practice, particularly when image quality is poor due to perturbation such as noise, partial volume artifact, or IIH. Beside that, not every clinical department can afford the computer cost to run the parallel LCJ algorithm.
Liew and Yan [58] approximated the IIH as a stack of Bspline surfaces with continuity constraints across slices. The estimation of IIH interwines with a fuzzy c-means clustering process. In [61,62], segmentation that utilizes local scale as homogeneous criteria has been presented and applied to IIH correction.
When a statistical classifier, such as Gaussian classifier or random field modeling, is exploited [50,53,63], the process is similar to the parameter estimation in Section 4, where the parameters are associated with a probability distribution and can be estimated with common statistical estimation methods such as maximum-likelihood estimation.

Entropy minimization
As a frequently used criterion to characterize the intensity distribution of an image, entropy has been employed to design algorithms for image restoration, thresholding, or classification [64,65]. Also, it has been utilized to quantify the image property with IIH present and guide the parameter searching for IIH removal [54][55][56][57].
It is assumed that the intensity distribution of the original image is multimodal, and the presence of IIH causes the intensity overlapping between objects. Figure 3 shows the histograms of brain tissue in a BrainWeb-simulated image (Figure 3(a)) and the head image in Figure 1 (Figure 3(b)). On Figure 3 correction can be achieved through searching the parameter space of the IIH model such that the entropy of the image is reduced. It should be pointed out that direct minimization on the entropy would lead to the null field [55,66]. To avoid this pitfall, constraints over the solution space are necessary. Mangin [55] constrained the solution to minimize the distance between the mean values of the restored and the original image. In [57], the restored image was constrained to have the same mean value as the original one.
Evidently, other quantities relevant to image features, variance for example, can also be applicable in a similar fashion. Again, constraints upon the solution space are necessary.

STATISTICAL MODELING
The statistical methods [67][68][69][70][71] may assume that the IIH follows a distribution, the Gaussian distribution for example, or model the IIH as a random process, such as the Markov random field.

Bayesian framework
The Bayes' rule has frequently been employed to estimate the IIH map when the IIH is modeled by a distribution. Let β be a random vector (β 1 , . . . , β n ) with probability density p(β).
To estimate β, one can maximize the conditional probability of β given y (the log-transform of x) as follows: This is called the maximum a posterior (MAP) estimate and, by the Bayes rule, is equivalent to Wells et al. [67] used the Gaussian distribution to model the entire log-transformed bias field and the observed intensity at voxel i: where Γ i is the tissue class at voxel i with mean value μ(Γ i ), and with ψ x as the covariance matrix. By assuming the statistical independence of voxel intensities and from (7), one can derive When the image data is not polluted by IIH, the above method is simply the tissue classification using a mixture Gaussian model. Hence, this method essentially interleaves the IIH correction with a Gaussian classifier. Guillemaud and Brady [68] observed that the effect of IIH correction by Wells et al. algorithm is substantially affected by the Gaussian classifier. In real images, it is very possible for the histogram to deviate from the mixture Gaussian distribution. A modification was then proposed by introducing a tissue class Γ other with a non-Gaussian distribution Zujun Hou 5 With this modification, the IIH is only estimated with respect to the Gaussian classes. Including the non-Gaussian component makes the Gaussian classifier less influenced by possible outliers arising in images.

Spatial modeling
In the method by Wells et al. the IIH correction has not explicitly considered the context of IIH map. Since the IIH field is slowly varying in the image domain, the values of the IIH map in neighboring pixels/voxels would be close. When this spatial relation is taken into account in the IIH modeling, one would likely arrive at a smoother approximation of the IIH map. A useful tool for spatial modeling is the Markov random field (MRF), which is first employed by Geman and Geman [72] for image segmentation. Held et al. [69] have employed the MRF to model IIH. According to the Hamersly-Clifford theorem [73], the prior probability p(y) is given by the Gibbs distribution with the Gibbs energy where i, j sums over every voxel i and its neighbours j.

The N3 method
Different from most IIH correction method which involves a classification step, Sled et al. [74] proposed a nonparametric nonuniform intensity normalization (N3) method which searches for the IIH field to maximize the frequency content of the image intensity distribution. The method simplified the problem in log-domain as a deconvolution problem by realizing that if v 1 and v 2 are two independent random variables with distributions V 1 and V 2 , respectively, then the distribution of their sum is the convolution of V 1 and V 2 [75].
To constrain the solution space, the IIH field is modeled as a Gaussian distribution with small variance. Code for this method is publicly available. 1

Comparison between local and global statistics
There are efforts [76,77] to estimate the IIH by comparing a local statistic with the global one. The two statistics are assumed to characterize the same population. These methods essentially relate the IIH correction to tissue segmentation, and the implicit assumptions are (1) the constant intensity 1 http://www.bic.mni.mcgill.ca/software/N3 for a tissue and (2) that the intensity variation within a tissue is solely due to IIH. Not surprisingly, these methods are sensitive to the estimation of tissue statistics, which is nontrivial in practice. When GM and WM are combined as one class and the local statistic is estimated from a sample in a local region as done in [76,77], the method can be regarded as a generalized white matter method by Dawant et al. [48], where the reference tissue is the combination of GM and WM. Also, it can be taken as a lowpass filtering estimation. Although it is usually much easier to identify GM and WM together than to identify WM alone, a potential problem is that the local sample could fail to adequately characterize the feature of the combined tissue class even though there is no artifact like IIH.
A possible solution is to carry out more detail tissue classification in each local region. For example, the technique proposed in [78,79] 2 estimated the global tissue mean values by empirical thresholding, while the local statistics are derived through fitting the local histogram with a theoretical distribution. After obtaining estimations of local correction factors, a smooth function is fitted among these data and applied to the whole brain volume for IIH removal.

Image feature-based methods
An image feature-based technique was reported in [80], where the IIH correction was decomposed into row and column correction. The correction factor at a voxel is firstly related to first-order difference at other voxels in the same row/column, then combined with those calculated in the rows/columns from the initial to the current one. The rationale underlying the computation is obscure from the description. However, it is very similar to that in [81]. In the latter technique, a smooth "variation" image was firstly derived from normalized intensity gradient field, where pixels with low intensity or high gradient magnitude are excluded. Then numerical integration was applied to the "variation" (first-order derivative) image to obtain an image, which only contains small variations and was used to determine the IIH map.
Vovk et al. [82,83] proposed to use the probability distribution of image features for IIH correction, where the image feature includes the intensity and the second spatial derivative of the image. Similar to the usual intensity histogram, the joint probability distribution also contains information for classifying tissues and such information was encoded by entropy. The correction factor was derived so that the entropy would decrease, similar to [57].

Estimate without explicit modeling
There are methods [84][85][86][87][88] which consider the IIH as model parameters formulated in a segmentation framework. Suppose the segmentation is to optimize a functional Φ(y, θ, β), where y denotes the observed data, β the IIH term, and θ 6 International Journal of Biomedical Imaging other parameters. Then one way to estimate β can be obtained by Rajapakse and Kruggel [84] used an MRF formulation, whereas Farag and his group exploited the fuzzy c-means clustering framework [85][86][87][88]. It is noted that in these methods no explicit assumption has been made on the IIH field, which can be advantageous over model-based methods, since the assumptions with a model could be violated in practice. On the other hand, the absence of constraints on the IIH solution could result in erroneous IIH maps that deviate far away from the truth. In addition, voxel-wise updating the IIH in an iterative process is time-consuming, hence techniques such as multigrid computing may help reduce the computation load.

Registration against template
Image registration has also been utilized to aid the IIH correction [89,90]. In [90], the patient data was registered against a tissue reference template, which allows to estimate the IIH map by direct comparison between two images. Here human intervention was employed to ensure the correct correspondence between the distorted and the reference image.

Shape recovery
Lai and Fang [91] transformed the IIH correction into the problem of shape recovery with orientation constraint and solved the latter using regularization theory. The approach may result in solving a linear equation with a large matrix.

Deformed thin plate model
Bansal et al. [92] modelled the IIH field as a thin plate deforming elastically under a body force: where μ and λ are the elasticity constants. The body force b(β) is evaluated to minimize the entropy of the observed image.

Integrated approaches
As seen from the description above, many IIH correction methods relate the problem with image segmentation and solve these two problems alternatively through an iteration process. Evidently, accurate segmentation would significantly ease the burden of IIH correction. Conversely, if the IIH has been precisely removed, the segmentation accuracy will in general be improved. Thus, it is not surprising to see overwhelming techniques addressing these two problems within a common framework.
Actually, it has been a trend in computer vision to imitate the human intelligent system and solve the different problems simultaneously. A typical computer vision system may consist of several individual processes, which can be solved sequentially. However, the solution of one process could be beneficial to the solution of another one and the converse may also hold. As an example, image denoising and edge detection are two closely related problems. And it is common to require a denoising algorithm able to preserve image edge structures, and an edge detection method robust against noise. For the problem of IIH correction, beside the connection with image segmentation as mostly noted, there are efforts that relate the solution to image registration, because the image quality can impose an impact on the accuracy of image registration and conversely a good registration against a template could greatly help derive a high-quality image. In addition, there even exist attempts to address these three processes, IIH correction, segmentation, and registration together [93,94].
Among some processes, their relationship may not be very intimate, but different implementation order could result in different performance. In [68], it was found that denoising after IIH correction is more preferable. Madabhushi and Udupa [95] investigated the interplay between IIH correction and intensity standardization, and concluded that the better sequence is IIH correction followed by intensity standardization.

Validation and comparative study
For end users, it is natural to ask questions such as how to assess the performance of an IIH correction method, which method should be recommended when a practical medical image processing system encounters the problem of IIH correction, or if there is a method which exclusively outperforms others. To answer these questions, we need to do extensive comparisons under a variety of data sets. However, this turns out to be a difficult task, because the true amount of IIH is unknown for real data. The lack of ground truth is a common problem in evaluating a computer vision algorithm. There are two possible ways to circumvent this difficulty. One is to approximate the golden standard by experts' estimation, which is often a tedious task. Alternatively, we can use synthetic data. In IIH correction, the simulated brain images [7,96] from the Montreal Neurological Institute 3 have been widely employed for validation.
(1) Criteria The criteria that have been frequently used are listed in the following.
Cr-I The variance of the fully or a partially segmented image, which is supposed to decrease after the IIH correction. When this criterion is used for comparing different methods, the result could be misleading because the variance is scale-variant. Usually, the mean-preserving condition is utilized to avoid this problem. Cr-II The coefficient of variation, cv, of class Γ i : It can be shown that this quantity overcomes the limitation of the image variance. But the cv alone only characterizes the within-class scatter and a criterion that also takes into account the between-class scatter is as follows. Cr-III The coefficient of joint variations between two classes Moreover, one can also use the relative change of c jv as defined below: where the subscripts a and b denote after and before IIH correction. Cr-IV Mean-square error, which directly measures the distance between the derived and the true IIH map. Cr-V Segmentation accuracy, which indirectly reflects the effect of IIH correction. Care should be taken in interpreting the segmentation result since the latter could be complicated by other factors, like subject, scanner, noise, segmentation method, and so forth. Cr-VI Stability, which means that an IIH correction algorithm is recursively applied to the corrected volume. For a good algorithm, the extracted IIH map is assumed to converge rapidly. Cr-VII Computer requirement and CPU time.
From the list, one can observe that most criteria have their own limitations and some are applicable to simulated data only. However, simulated data might not adequately characterize real ones. For example, in [47], the proposed method was reported to be inferior to methods such as the N3 in terms of cv or c jv when tested on simulated data, but the order is reversed when tested against real volumes. Furthermore, the adaptivity of an algorithm is also important. For a method with good adaptivity, the approximated IIH map would approach a constant when the real IIH approaches zero.
(2) Comparative study Compared to the numerous techniques for IIH correction, only a few studies have been carried out towards the comparative evaluation of existing algorithms. Sled et al. [97] have compared three IIH correction methods, the expectation maximization (EM) [67], the white matter (WM) [48], and the N3 method [74] using simulated T1, T2, and PD weighted data. It was shown that the WM method performs better than the other two methods for T1 weighted volumes, which might be due to the high contrast between the WM and other tissues in T1 weighted images. The EM method made excessively large corrections to voxels that fall outside the classifier's tissue model, as is consistent with that pointed out in [68]. Overall, the N3 method performs the most stable for all simulated images.
Velthuizen et al. [29] have evaluated four IIH correction methods (a phantom method [17], two lowpass filtering methods [36,39], and a surface fitting method with reference points selected from white matter [48]) in brain tumor segmentation. The surface fitting method was found to be inferior to others, which could be due to the way the reference points were generated. As mentioned in Section 3.1, the latter is crucial to the performance of the surface fitting method. An automatic method to generate such reference points has been presented in [59]. Hou and Huang [98] have also developed a similar technique based on order statistics, which is very well comparable with the state-of-the-art IIH correction methods.
Although it turned out no improvement in tumor assessment after the IIH correction [29], it does not mean that IIH correction is not an obstacle to automatic medical image processing in general, since the tumor segmentation is characterized by the localization of the tumor region as well as the intensity contrast with surrounded tissue. Thus, the tumor segmentation could be less affected by the IIH artifact.
A more comprehensive study was presented in [66], where six algorithms, n3 [74], hum [42], eq [43], bfc [78], spm (statistical parametric mapping) 4 [52], and cma 5 were compared against BrainWeb-simulated data as well as real volumes including repeated scans of the same subject, scans under different magnetic fields and different scanners. Three of the methods (hum, eq, and cma) are lowpass filtering based. The spm method is based on surface-fitting, and its parameters are estimated through integration with a tissue mixture model. It was found that the IIH maps obtained by filtering based methods can exhibit higher-frequency structures pertaining to brain anatomy. The spm method could be unstable when operating on relatively uniform image volumes and could lead to spurious solution for some volume. Overall, the n3 and the bfc methods are superior to the other four methods. At lower bias levels, the estimated bias by bfc is more accurate than that by n3, and at higher bias levels, the order reverses. Nevertheless, none of the six methods performs ideally under all the circumstances investigated.
The problem of the spm might be similar to that of the EM method by Wells et al. [67]. Both methods utilized the mixture Gaussian classifier, which may be inadequate to model the image intensity distribution arising in practice. It should be pointed out that the spm method used in [66] is the SPM99 version, which has been updated to version SPM2 in 2003 with substantial improvement in theoretical modeling 8 International Journal of Biomedical Imaging or algorithmic design, and the latest version is SPM5. As to the three filtering-based methods, they lack a scheme to adapt the filtering strength to data quality, which may explain the inefficiency compared with the n3 and the bfc methods. As mentioned in Section 2, filtering methods [45][46][47] with data adaptivity have been developed recently, which might outperform their conventional counterparts.
Although further comparative study using more extensive MR images is necessary, it might be inappropriate for end users to expect an algorithm superior to others and exclusively applicable. In general, each method has its underlying assumptions and limitations and the choice of which method to use is intimately interwined with the problem to solve, the source, and quality of the data. Many methods have attempted to correct the IIH artifact in brain MR images, some of which require the removal of the scalp/skull before the correction process, while others do not. Although sophisticated methods may be able to correct for IIH more accurately, one would also have to consider the expense of computer cost as well as the final segmentation error. Among the publicly available softwares, the N3 method has been widely used and its performance has been well demonstrated, while the BFC method can be advantageous when the image is also contaminated by severe noise [98].

CONCLUSION
This paper presented a summary of the recent progress on MR image IIH correction. The most popular models to describe the IIH field are the low frequency, the hypersurface, and the statistical model. Filtering methods are fast, easy to code and widely used. With optimization in scale space, the filtering method can also be adaptive to image data. Surface fitting and statistical methods are easy to integrate with other knowledge such as segmentation, registration, or some image feature, thus could in principle provide more reliable solution, which have been and will be the trend in the field. Some techniques based on other IIH correction principles were also reviewed in the paper. In future, it might be of interest to have more extensive investigations on evaluation of existing methods.