Another Bias Correction Method

1. Another Bias Correction Method John Ashburner. Wellcome Department of Imaging Neuroscience, London, UK. john@fil.ion.ucl.ac.uk

2. 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. Another Bias Correction Method, J Ashburner, 2002

3. 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 (mk), variance (sk2) and mixing proportions (gk, where Skgk=1 and gk>=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,mk,sk) = (2psk2)-1/2 exp(-(yi-mk)2/(2sk2)) • To incorporate bias correction, extra parameters are used to model smooth intensity variations. A field modelling the variation at element i is denoted ri(b), where b is a vector of unknown parameters. Intensities from the ith cluster are assumed normally distributed, with mean mk/ri(b) and variance (sk/ri(b))2: P(yi|k,mk,sk,b) = (2p(sk /ri(b))2)-1/2 exp {-(yi-mk /ri(b))2/(2(sk/ri(b))2)} = ri(b) (2psk2)-1/2 exp {-(yi-mk /ri(b))2/(2sk2)} Another Bias Correction Method, J Ashburner, 2002

4. Mixture of Gaussians Based Derivation - II • Probability of a voxel, irrespective of intensity, belonging to kth Gaussian, given the proportion of voxels belonging to that Gaussian is P(k|gk)=gk . • Therefore, the joint probability of cluster k and intensity yiis: P(yi,k|mk,sk,gk,b) = P(yi|k,mk,sk,b)P(k|gk) • By integrating over k, we obtain the likelihood of yi: P(yi|m,s,g,b) =SkP(yi,k|mk,sk,gk,b) =SkP(yi|k,mk,sk,b)P(k|gk) • Likelihood of the entire dataset is derived by assuming independence: P(y|m,s,g,b) = Pi P(yi|m,s,g,b) = Pi{Sk ri(b) (2psk2)-1/2 exp {-(yi-mk /ri(b))2/(2sk2)} } • Likelihood is maximised when the following is minimised: E = -log{P(y|m,s,g,b)} =-Si log{Sk ri(b) (2psk2)-1/2 exp {-(yi-mk /ri(b))2/(2sk2)} } Another Bias Correction Method, J Ashburner, 2002

5. Non-parametric Generalisation • 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 (g). By representing the density of the kth Gaussian at yi by fk(yi), the negative log-likelihood of the model is: E = -log{P(y|g)} = -Si log{Skgkfk(yi)} • In fact, fk(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|g,b)} = -Si log{ri(b) Skgkfk(ri(b)yi)} Another Bias Correction Method, J Ashburner, 2002

6. The algorithm is iterative, and involves alternating between: Re-estimating g, while holding b fixed. Re-estimating b, while holding g fixed. Re-estimating g involves building a probability density representation of the data, which have been corrected with the current bias estimate. Simplest case involves generating a simple histogram. An iterative method is needed for B-splines of degree greater than 1. Similar to algorithm for ML-EM reconstruction of PET images. Re-estimating b makes use of the histogram and its gradient w.r.t. intensity - therefore zeroeth degree splines can not be used. B-spline Density Representations Zeroeth (nearest neighbour), first (linearly interpolated), second and third degree B-spline density representations. Another Bias Correction Method, J Ashburner, 2002

7. Regularisation + Optimisation • Distinction between intensity variations due to: • bias artefact due to MR physics • different tissue properties. • Don’t want to over-fit the data, so regularisation is used. • Maximise the joint probability of y and b, assuming independence ofb and g: P(y,b|g) = P(y|b,g) P(b) • where P(b) is assumed to be a d dimensional multi-normal distribution: P(b) = N(b0,C0) = ((2p)d|C0|)-1/2 exp{ -1/2 (b-b0)TC0 -1 (b-b0)} • 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. Another Bias Correction Method, J Ashburner, 2002

8. A Quick Evaluation The method was applied to the 40% non-uniformity, 3% noise, 1mm slice thickness simulated T1 BrainWeb (Kwan et al, 1996) image from: http://www.bic.mni.mcgill.ca/brainweb/ . When applied only to a region containing brain and CSF, the simulated and recovered biases have a correlation coefficient of 0.9964. Applied Bias Estimated Bias Note that intensities are uniformly scaled to match the estimated bias. Histograms of data before and after recovery - original intensities - recovered intensities Another Bias Correction Method, J Ashburner, 2002

9. 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-1Si log{ri(b) yiSkgkfk(ri(b)yi)} = I-1 E - I-1 Si 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. Another Bias Correction Method, J Ashburner, 2002

10. 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. Another Bias Correction Method, J Ashburner, 2002