Image background inhomogeneity correction in MRI via intensity standardization

1. Introduction 

Background image intensity variations caused by imperfections in imaging devices pose major challenges for further processing, segmentation, and analysis of acquired images. In magnetic resonance imaging (MRI), such effects originating from imperfections in the radio frequency field and from the human body are well known [1]. A variety of techniques have been developed over the past 15 years to address this problem [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13]. While many of these methods have provided an effective solution in MRI, there is room for improvement in the sense of developing a general method that fulfills the following requirements: (R1) no need for user help especially on a per-scene basis; (R2) no need for accurate prior segmentation of objects in the image; and (R3) no need for prior knowledge of intensity distribution.

Early methods for intensity inhomogeneity estimation and correction used phantoms to empirically measure the inhomogeneity [6]. Wells et al. [8] proposed an expectation maximization (EM)-based algorithm consisting of interleaved voxel classification and inhomogeneity correction steps. Guillemaud and Brady [9] introduced a modified EM algorithm that replaced the distribution of the class other, which included all tissues not explicitly modeled, by a uniform probability density function. The correction was claimed to be more robust and to overcome some limitations of Wells’ original method. Similar to these works [8], [9], Van Leemput et al. [10] proposed a model-based method that also employed an iterative EM strategy that interleaved voxel classification with estimation of class distribution and inhomogeneity parameters. However, they used a priori information from a digital brain atlas to initialize the expected location and parameters of tissue classes, and claimed yielding more objective and reproducible results. Sled et al. [11] described a non-parametric non-uniform intensity normalization () method. The method is independent of MRI pulse sequences and is claimed to be insensitive to pathological data that may otherwise violate model assumptions. An iterative approach was employed to estimate both the multiplicative inhomogeneity and the distribution of the true tissue intensities. Styner et al. [12] developed a correction method called PABIC based on a simplified method of the imaging process, a parametric model of tissue class statistics, and a polynomial model of the inhomogeneity field. Likar et al. [13] described a model-based correction method that made the assumption that an image corrupted by intensity inhomogeneity contains more information than the corresponding uncorrupted image. The image degradation process was described by a linear model consisting of a multiplicative and an additive component which were modeled by a combination of smoothly varying basis functions.

It is clear that some object information is vital in order to estimate (and, therefore, to correct) the background variation component of the intensities in the image. This, however, we argue, need not come in the form of explicit object segmentation or by way of prior tissue intensity distribution. To liberate a solution to inhomogeneity correction from the vagaries of image segmentation will be a good idea, since the latter itself is a difficult problem which has defied a solution fulfilling (R1)–(R3) up till now. What is needed is sample (not necessarily complete) object regions constituting the same object material (tissue in medical imaging). In this paper, we present a significant body of evidence that a solution fulfilling (R1)–(R3) is indeed feasible through the use of intensity standardization techniques, introduced in [14], [15].

A major problem in MR images is the non-standardness of the MR image intensity gray scale, which implies the lack of a tissue-specific numeric intensity meaning, even within the same MRI protocol, for the same body region, for images obtained on the same scanner, and for the same patient. The intensity standardization method described in [14] offers a simple way of transforming the MR images nonlinearly so that there is a significant gain in similarity of the resulting images. We have demonstrated that this technique is very useful in segmentation and tissue characterization [16], and that the same segmentation methods perform better after standardization than before [17]. The purpose of this work is to demonstrate how standardization can be employed to reduce background inhomogeneity in MR images. Our method, described in Section 2, consists of first standardizing the MR images to the standard intensity scale, then obtaining different tissue sample regions from simply thresholding by using fixed intensity intervals, and then estimating and correcting for intensity variations based on the thresholded regions, which constitute the sample regions mentioned earlier but not complete or accurate segmentations. The process is iteratively repeated on the resulting images. We demonstrate in Section 3 both qualitatively and quantitatively the effectiveness of the method in correcting even severe inhomogeneities coming from surface coils. Our concluding remarks are stated in Section 4. An early version of this paper was presented at the SPIE Medical Imaging 2006 conference whose proceedings contained an abbreviated version of this paper [18].

2. The standardization-based method 

2.3. Description of steps in SBC 

2.3.1. S1: standardizing MRI scene intensities 

As demonstrated in [14], MRI scenes in do not have scene intensity gray scales that can be considered to be nearly the same, even for scenes obtained on the same scanner, for the same subject. The method described in [14], and enhanced in [15], offers a way of transforming the scenes in so that there is a significant gain in the similarity of the resulting scenes in terms of their intensity-wise tissue characterizability. This is achieved in two steps—a training step that is executed only once for a given P and R, and a transformation step that is executed for each given scene in . In the training step, certain landmarks of a standard histogram are estimated from a given subset of . In the transformation step, the actual intensity transformation from the intensity scale of the input scene to the standard scale is computed by mapping the landmarks determined from the histogram of the given scene to those of the standard scale. The standardization process is illustrated in Fig. 1 wherein histograms of the white matter (WM) region in a set of 10 PD-weighted MRI scenes are shown before and after standardization. The effect of intensity standardization is readily seen by comparing Fig. 1(a) and (c), and (b) and (d). We use the method of [15] to produce the standardized scene corresponding to given scene to be corrected.

View Large Image Download to PowerPoint

Fig. 1. (a) Histograms of WM regions in 10 original PD-weighted MRI scenes. (b) A histogram obtained by combining the 10 histograms of (a) which represents WM intensity distribution in PD scenes. (c) Histograms of WM regions in the 10 corresponding standardized scenes. (d) Similar to (b) but obtained from the 10 histograms in (c).

2.3.2. S2: determining object regions 

After MRI scene intensity standardization, the scene intensity intervals for different tissues for the specific P and R are determined and fixed once for all, and used throughout the subsequent background inhomogeneity correction steps for all scenes in . Fig. 2 (a) shows a slice of a given standardized scene , which is a simulated PD-weighted MRI scene of the head [19], and Fig. 2(b) and (c) shows the gray matter (GM) and WM regions obtained by thresholding via fixed intensity intervals. The largest connected component of each thresholded result is extracted, and produced as a binary scene . This binary scene gives us the sample object region (although not with a perfect segmentation by any means) containing roughly the same object material for each object/tissue. Let be m binary scenes produced in this manner for m object/tissue regions in .

View Large Image Download to PowerPoint

Fig. 2. (a) One slice of the original PD-weighted MRI. (b) GM region and (c) WM region resulting from thresholding this slice with fixed intensity intervals; (d) corresponding slice of ; and (e) corresponding slice of corrected scene .

2.3.3. S3: stopping criterion 

The basic premise of this method is that, as correction is effected, the object region represented in each will change gradually. When this change is insignificant (), the iterative procedure stops. From our numerous experiments, we have observed that the iterative procedure stops within 20 iterations. The variations of and with iteration within SBC are illustrated in Fig. 3 for the scene example used in Fig. 2. The effect of inhomogeneity correction with iteration is shown via intensity histogram of the foreground only of the corrected scene in Fig. 3(d), where the overlap between the intensity distributions of GM and WM caused by background variations is reduced gradually with increasing number of iterations.

View Large Image Download to PowerPoint

Fig. 3. The improvement in observed with increasing number of iterations (left to right: 3, 5, 10, and 20) illustrated through the examples of Fig. 2. (a, b) thresholded scenes corresponding to GM and WM threshold intervals. (c) . (d) The intensity histogram of . The improved definition of the two modes corresponding to WM (larger mode) and GM with iteration is clear. For clarity of display, the histograms shown do not include the background.

2.3.4. S4: estimation of inhomogeneity map 

In this step, the background variation is estimated first and then corrected. We treat every tissue region from Step S2 as an independent entity with its own local, best fit polynomial to the background variation component. We then combine these independent polynomials into a single unified function and use it for correcting inhomogeneities in the whole scene. We explain our implementation below in detail.

For any tissue region , let be the set of all 1-valued voxels in , and be the mean intensity within . Without loss of generality, let . The background variation in each is then estimated as a discrete function defined as follows. For any ,

(1)

where is the operation (mentioned in assumption (A2)) between and that yields . We denote this operation by , which equivalently means, for any , . In MRI, is usually considered to be multiplication. A 2nd order polynomial is then estimated from the discrete function by using a least-squares-error fit.

As per assumption (A1), the background variation in a scene is slow varying. However, the inhomogeneity functions , are separately estimated for different homogeneous regions (corresponding to different tissues/objects) in an independent manner. Consequently, the continuity between the inhomogeneity map estimated over and over , cannot be guaranteed, as illustrated in Fig. 4. This discontinuity is caused by the fact that the expected intensity for region is corrupted independent of and . At this stage, it is impossible to determine the bias in the mean intensity in the context of global inhomogeneity. The situation is exacerbated by the presence of intensity non-standardness itself. Since MRI intensities do not have absolute tissue-specific meaning, it becomes impossible to make assumptions about the relative shifts (biases) among the estimated . This points out the arbitrariness that exists in determining the overall level of the correcting function. This also clarifies as to how inhomogeneity correction itself can introduce non-standardness [20]. Our approach to fit an overall correcting function to estimated maps is as follows, and also illustrated in Fig. 4.

View Large Image Download to PowerPoint

Fig. 4. Illustration of the possibility of existence of discontinuity between inhomogeneity maps estimated independently from different tissue regions and , and how to combine them into a single unified inhomogeneity map. For the convenience of presentation, the scene is considered to be 1D in this illustration, and only two objects are considered. (a) Original discrete inhomogeneities are estimated independently from different tissue regions, and two smooth inhomogeneity maps and are fitted for two objects and in an independent manner. (b) is shifted by a weight factor so that the difference between and is minimized. (c) A new discrete inhomogeneity map is estimated by combining and . (d) A single smooth unified inhomogeneity map is fitted to .

At first we combine the inhomogeneity maps and over the largest and the second largest regions and , respectively, as follows:

(1)Find a weight factor to minimize . Essentially, this step computes the bias to be applied to so that the difference between and is minimized.

(2)Combine the two inhomogeneity maps and to obtain a new discrete inhomogeneity map such that, for any

(2)

where is the minimum distance between c and a voxel in .

(3)Determine a 2nd degree polynomial that constitutes a least-squares-error fit to .

After merging and as above, we now have distinct regions and inhomogeneity maps . The above steps are then repeated until we have only one region and a single unified inhomogeneity map . This unified inhomogeneity map is represented as a scene , where .

Since is invertible, , and therefore, by (A2),

(3)

In MRI, is generally considered to be a multiplication operation, and therefore, becomes division. In general, if is known and invertible, (i.e., is known), then from Eq. (3), we can determine . In Fig. 2(d) and (e), we display and for the running example after the first iteration.

After scene correction, we go to the standardization step and repeat the procedure. We note that, the parameters of the standardizing transform are determined by training only once and then fixed for a given P and R (and so also the threshold intervals). The standardization step (actually the transformation step) S1 in each iteration uses these same parameters for effecting standardization. Therefore our strategy of using fixed threshold intervals in the standardized scene is perfectly legitimate. Since inhomogeneity correction methods introduce some non-standardness of their own [20], as described above, in each iteration of SBC, the standardization step S1 becomes necessary.

3. Results and evaluation 

In Section 3.1, to demonstrate qualitatively how effective SBC is in correcting background variation, we present a set of examples showing scenes before and after applying the SBC correction procedure. We present the results of a quantitative evaluation procedure in Section 3.2 based on 12 simulated MRI scenes obtained from the BrainWeb database [19] and on 10 PD- and 10 T2-weighted clinical MRI scenes with Multiple Sclerosis (MS) lesions obtained from our own database. Information about the scene data sets used in our evaluation is listed in Table 1. We compare the SBC method with the N3 method [11]. The main reasons for selecting N3 for comparison are: (i) it seems to be commonly used, (ii) it comes closest to SBC in its spirit of not requiring any user input or a priori information, and (iii) following the results of the comparative study among six algorithms for correcting intensity inhomogeneity reported by Arnold et al. [21], the N3 method generally outperformed other methods. Briefly, the basic idea of the N3 method is to find a smooth, slowly varying, multiplicative field that maximizes the frequency content of the probability density function of the logarithm of scene intensities. It is an iterative process in which, in addition to the intensity inhomogeneity, the non-parametric model of the scene intensities is also derived directly from the data. The implementation of the N3 method available at http://ftp.bic.mni.mcgill.ca/pub/mni_n3 is utilized in our comparison. In our experiments, the parameters of the N3 method were set to the default values.

Table 1. Information about scene data utilized in evaluation.

	Data set	Body region	Number of scenes	Scene domain	Voxel size (mm3)
	D1	Head-coronal	1
	D2	Head-transverse	1
	D3	Abdomen	1
	D4	Hips	1
	D5	Foot	1
	D6	Simulated brain MRI	12
	D7	Clinical brain MRI	10 T2+10 PD

3.1. Qualitative 

In Fig. 5, we present five examples constituting data sets D1–D5 of Table 1 involving different P and R. One slice from the original scene, and the same slice of the corrected scene produced by N3 and SBC methods are shown for each data set. In all displays, the same gray level window settings were used. The corresponding intensity histograms (for the whole 3D scene) are displayed under each scene. These figures give a qualitative indication of the effectiveness of SBC in suppressing background variation. In all examples, SBC worked automatically requiring no specific setting or adjustment of parameters, and so did N3. Since SBC has no parameters whose values are to be adjusted or ascertained for different situations by user interaction or by using priors, this automatic mode became feasible on scenes listed in Table 1 and on a variety of others. From the histograms displayed, SBC seems to have performed better than N3 in all cases yielding a better definition of the modes in the histogram.

View Large Image Download to PowerPoint View Large Image Download to PowerPoint

Fig. 5. A slice of MRI scenes of different body regions acquired by using different MRI protocols is shown on the left in (a)–(e). The corresponding slices of the scenes after inhomogeneity correction by using the N3 (middle) and SBC procedures (right). For each of (a)–(e), the corresponding intensity histogram (of the whole 3D scene) is shown under each scene. The data sets in (a)–(e) correspond to D1–D5 in Table 1.

3.2. Quantitative 

We assess the performance of the correction methods by comparing the percent coefficients of variations () of intensity within individual tissue regions in the original and corrected scenes. After the scene is segmented, of tissue class K is computed as

(4)

where and are the standard deviation and mean intensity of class K, respectively. The quantity , which represents the normalized standard deviation in a given tissue class, is invariant to uniform multiplicative intensity transformation.

3.2.1. Simulated data 

Twelve (four T1, four T2, and four PD-weighted) MRI scenes from the BrainWeb database were used [19](D6 in Table 1). All scenes had noise. Six scenes each with and inhomogeneity were used, three of each were normal and three with MS lesions. Table 2 lists the results. The corrected scenes indicate that both N3 and SBC can effectively minimize the intensity variations within GM and WM. Generally, the values of all corrected scenes for inhomogeneity are similar to the of the scenes for inhomogeneity for both methods. One slice of a T2-weighted MRI scene with MS lesions before and after applying the N3 and SBC methods and the corresponding (3D) intensity histograms are shown in Fig. 6. Generally the values for SBC seem to be as good as those for N3. Since multiple data sets with multiple anatomies are not available in this database under each protocol, we cannot conduct a statistical comparison.

Table 2. Percent coefficient of variation of tissue intensities in segmented GM and WM regions for 12 simulated MRI scenes before correction, and after correction by N3 and SBC methods.

	Date set	Protocol	Inhomogeneity	(GM)			(WM)
				Original			Original
	D6 (normal)	T1	20%	11.0	9.9	9.9	6.7	5.1	5.1
			40%	13.5	10.0	9.9	9.2	5.2	5.2
		T2	20%	18.4	18.0	17.9	12.0	11.8	11.7
			40%	20.3	20.0	17.9	13.3	13.0	11.8
		PD	20%	6.3	4.6	4.5	5.5	4.7	4.7
			40%	9.7	4.6	4.5	7.5	4.6	4.6
		T1	20%	11.2	10.1	10.1	6.9	5.3	5.3
			40%	13.7	10.2	10.1	9.3	5.3	5.4
	D6 (lesions)	T2	20%	10.9	10.0	9.8	8.7	8.3	8.1
			40%	13.7	10.1	9.8	10.6	8.2	8.2
		PD	20%	5.8	3.9	3.8	5.3	4.4	4.3
			40%	9.4	4.0	3.9	7.4	4.3	4.3

View Large Image Download to PowerPoint

Fig. 6. (a) A slice of a simulated T2-weighted MRI scene with MS lesions and with inhomogeneity. The corresponding slice of the scene corrected by using (b) N3 method, and (c) SBC method. The scene intensity histograms are shown in the second row.

3.2.2. Clinical data 

Clinical MRI scene data (D7 in Table 1) used in our evaluation experiments came from 10 different MS patients and each had T2- and PD-weighted scenes. These data sets were randomly selected from our database, and were acquired on the same GE Signa 1.5T scanner with a FOV of 22 cm, and TR of 2500 ms and TE of 18ms/90ms. Again the tissues considered were GM and WM. Fig. 7 shows one slice of one PD-weighted scene and its corrected scenes created by using the N3 and SBC methods. The GM and WM regions were obtained by using the fuzzy connectedness-based segmentation method [22] implemented in 3DVIEWNIX [23], and by correcting the results when necessary under the guidance of an expert. Table 3 lists the mean and standard deviation (in parentheses) of values for both T2- and PD-weighted scenes, and for GM and WM regions. We observe that, for both tissue regions and for both protocols, the two methods achieve smaller values than before correction. Further, the values for the SBC method are smaller than those for the N3 method, and this difference has been found to be statistically significant under a paired -test for each protocol at a level of .

View Large Image Download to PowerPoint

Fig. 7. (a) A slice of a MRI PD scene of the head of a patient with Multiple Sclerosis. The corresponding slice of the corrected scene output by (b) N3, and (c) SBC methods, and the respective histograms (bottom).

Table 3. The mean and standard deviation (in parentheses) of percent coefficient of variation of tissue intensities in segmented GM and WM regions for 10 clinical T2- and PD-weighted MRI scenes of MS patients before correction, and after applying the N3 and SBC correction methods.

	Protocol	(GM)			(WM)
		Original			Original
	T2	16.7 (1.61)	14.9 (1.18)	14.7 (1.23)	12.9 (1.12)	11.5 (1.07)	11.2 (0.98)
	PD	7.1 (0.32)	5.9 (0.20)	5.6 (0.19)	7.8 (0.72)	6.6 (0.62)	6.2 (0.51)

4. Concluding remarks 

Background variation in MRI scenes due to imperfections in imaging devices is a problem commonly encountered in image analysis in many areas. Although many solutions have been proposed in the past, they all had varying degrees of dependency on the particular protocol, body region, application domain, and the type of variation. In an attempt to liberate this dependency, we have devised a method based on MRI intensity standardization. The method has no parameters whose values are to be determined or to be provided by the user on a per scene basis. Our qualitative and quantitative evaluation experiments indicate that its performance and behavior match our intuitive expectation. Its performance seems, at the least, to be comparable to that of the state-of-the-art N3 method, and our experiments indicate a slight edge for the proposed SBC method over N3. The N3 method assumes that the bias field is normally distributed, which may not be true for some clinical MRI data. In addition, the N3 method treats all voxels alike for inhomogeneity estimation. This is somewhat unnatural, since it is obvious that the WM voxels in a brain MRI scene, which have a narrow intensity histogram, are much more suited for inhomogeneity estimation than, for instance, the tissues surrounding the brain or ventricular CSF [10]. Other advantages of SBC are its simplicity, intuitive nature, and standardization of the scene as a byproduct of the method.

As related to assumption (A2), some comments are in order. There are basically the following three types of operations commonly done on MRI scenes before further carrying out other image processing and analysis operations on them: correction of background variation, filtering to suppress noise [24], and intensity scale standardization. There is perhaps some interaction among the effects of these operations, which, to our knowledge, has not been studied. To what extent the order in which these operations are carried out influences the results is also unknown. Given the importance of these operations, it is worthwhile to study the above issues. We have done some preliminary work along these lines [20], [24].