Wellcome Department of Imaging Neuroscience, London, UK.

john@fil.ion.ucl.ac.uk

Introduction

MR images are usually corrupted by a smoothly varying artefact that modulates the intensity of the image. This poster presents a method for removing most of this artefact.

Most current methods are either:

Parametric: Bias correction is incorporated into a mixture of Gaussians type approach, possibly as a refinement of a tissue classification algorithm.

Bias correction component is normally applied to log-transformed intensities (e.g. Wells III et al, 1996; Van Leemput et al, 1999).

Non-parametric: Bias correction applied to histograms of intensities in order to maximise entropy.

Most widely used approach is applied to histograms of log-transformed intensities (Sled et al, 1998).

Another approach minimises entropy of the histogram of the original intensities, with a modification to preserve the average intensity in the image (Mangin, 2000).

This poster will present a non-parametric approach, based on optimising an objective function similar to the entropy of the log-transformed intensity distribution, but using histograms of non-transformed intensities.

Mixture of Gaussians Based Derivation - I

A distribution can be modelled by a mixture of Gaussians (MOG). For univariate data, the kth Gaussian is modelled by a mean (k), variance (k2) and mixing proportions (k, where kk=1 and k>=0). Fitting a MOG involves maximising the likelihood of the data (y), given the parameterisations. The likelihood of a datum with intensity yi, given that it belongs to the kth Gaussian is:

P(yi|k,k,k) = (2k2)-1/2 exp(-(yi-k)2/(2k2))

To incorporate bias correction, extra parameters are used to model smooth intensity variations. A field modelling the variation at element i is denoted i(), where is a vector of unknown parameters. Intensities from the ith cluster are assumed normally distributed, with mean k/i() and variance (k/i())2:

Probability of a voxel, irrespective of intensity, belonging to kth Gaussian, given the proportion of voxels belonging to that Gaussian is P(k|k)=k .

Therefore, the joint probability of cluster k and intensity yi is:

P(yi,k|k,k,k,) = P(yi|k,k,k,) P(k|k)

By integrating over k, we obtain the likelihood of yi:

Consider a simple mixture model, where the means and variances are known, and therefore fixed by the model. Also assume the means are equally spaced, and the variances are identical. The only unknown parameters are the mixing proportions (). By representing the density of the kth Gaussian at yi by k(yi), the negative log-likelihood of the model is:

E = -log{P(y|)} = -i log{k k k(yi)}

In fact, k(yi) does not need to be a Gaussian density. It could be a B-spline with local support. In the limiting case of zeroeth and first degree B-splines, maximising this cost function is equivalent to generating a simple histogram of the data. E has a simple relationship with the entropy.

The bias correction model will also fit within this framework:

E = -log{P(y|,)} = -i log{i() k k k(i()yi)}

B-spline Density Representations

The algorithm is iterative, and involves alternating between:

A Newton-Raphson optimisation scheme can be derived from this:

First and second derivatives of E can be efficiently computed by parameterising the bias field in terms of separable basis functions.

A Quick Evaluation

Discussion

Optimising the current objective function is equivalent to optimising the entropy of the probability distribution of log-transformed intensities.

Entropy of log-transformed data is:

H = -I-1i log{i() yi k k k(i()yi)} = I-1 E - I-1 i log{yi}

where I is the number of voxels.

Histograms are produced from intensities that are not log-transformed.

Same method can be applied to images containing low or negative image intensities.

Resolves a problem pointed out by Arnold et al (2001), with the SPM99 bias correction approach (Ashburner, 2000):

The entropy of non-transformed intensity distribution is maximised when the estimated bias field is uniformly zero.

Resulted in side effects because average bias (rather than average corrected image intensity) was constrained to one.

Accuracy of the results depends on form and magnitude of regularisation.

References

J.B. Arnold, J.S. Liow, K.A. Schaper, J.J. Stern, J.G. Sled, D.W. Shattuck, A.J. Worth, M.S. Cohen, R.M. Leahy, J.C. Mazziotta & D.A. Rottenberg. Qualitative and quantitative evaluation of six algorithms for correcting intensity nonuniformity effect. NeuroImage 3(15):931-943, 2001.

J. Ashburner & K.J. Friston. Voxel-based morphometry - the methods. NeuroImage 11:805-821, 2000.

R.K.-S. Kwan, A.C. Evans & G.B. Pike. An Extensible MRI Simulator for Post-Processing Evaluation. Visualization in Biomedical Computing (VBC'96). Lecture Notes in Computer Science, vol. 1131. Springer-Verlag, 1996. 135-140

J.-F. Mangin. Entropy minimisation for automatic correction of intensity nonuniformity. In Proc. IEEE Workshop on Mathematical Methods in Biomedical Image Analysis, 2000.

J.G. Sled, A.P. Zijdenbos & A.C. Evans. A non-parametric method for automatic correction of intensity nonuniformity in MRI data. IEEE Transactions on Medical Imaging 17(1):87-97, 1998.

K. Van Leemput, F. Maes, D. Vandermeulen & P. Seutens. Automated model-based bias field correction of MR images of the brain. IEEE Transactions on Medical Imaging 18(10):885-896, 1999.

W.M Wells III, W.E.L Grimson, R. Kikinis & F.A. Jolesz. Adaptive segmentation of MRI data. IEEE Transactions on Medical Imaging 15(4):429-442, 1996.