Biostatistics

Biostatistics Advance Access originally published online on April 11, 2007
Biostatistics 2007 8(4):784-799; doi:10.1093/biostatistics/kxm005

This Article Abstract FREE Full Text (PDF) Supplementary Material All Versions of this Article: 8/4/784    most recent kxm005v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Request Permissions Disclaimer Google Scholar Articles by Frandsen, J. Articles by Vedel Jensen, E. B. Search for Related Content PubMed PubMed Citation Articles by Frandsen, J. Articles by Vedel Jensen, E. B. Social Bookmarking          What's this?

© The Author 2007. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oxfordjournals.org.

Bayesian regularization of diffusion tensor images

Jesper Frandsen

Department of Neuroradiology, Centre for Functionally Integrative Neuroscience, Århus University Hospital, Århus, Denmark

Asger Hobolth

Bioinformatics Research Center, Aarhus University, Århus, Denmark

Leif Østergaard and Peter Vestergaard-Poulsen

Department of Neuroradiology, Centre for Functionally Integrative Neuroscience, Århus University Hospital, Århus, Denmark

Eva B. Vedel Jensen^*

The T.N. Thiele Centre of Applied Mathematics in Natural Science Department of Mathematical Sciences, Aarhus University Ny Munkegade, DK-8000, Århus C, Denmark eva{at}imf.au.dk

^* To whom correspondence should be addressed.

	SUMMARY

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
Diffusion tensor imaging (DTI) is a powerful tool in the study of the course of nerve fiber bundles in the human brain. Using DTI, the local fiber orientation in each image voxel can be described by a diffusion tensor which is constructed from local measurements of diffusion coefficients along several directions. The measured diffusion coefficients and thereby the diffusion tensors are subject to noise, leading to possibly flawed representations of the 3-dimensional (3D) fiber bundles. In this paper, we develop a Bayesian procedure for regularizing the diffusion tensor field, fully utilizing the available 3D information of fiber orientation. The use of the procedure is exemplified on synthetic and in vivo data.

Keywords: Bayesian regularization; Diffusion tensor imaging; Fiber tracking; Markov chain Monte Carlo; Prior model; Rice distributions; Tensors

	1. INTRODUCTION

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
The human brain may be divided into 2 main components, gray matter and white matter. The gray matter contains the functional centers of the brain and has been studied for the last 2 decades by magnetic resonance imaging. This type of imaging technique is essential for understanding how we perform cognitive tasks and may also be used to study a wide range of diseases. The white matter connects the functional centers of the brain. The integrity and course of white matter fiber bundles are of key importance in understanding the structural basis of the functional integration of cortical centers in cognitive tasks, of the origin of functional impairment in focal brain lesions, and finally brain plasticity.

With the development of diffusion tensor imaging (DTI), white matter microstructures can be indirectly probed by measuring the directionality of Brownian movements of free water (Basser and others, 1994; Basser and Pierpaoli, 1996; Pajevic and Basser, 2003). By the hindrance of water diffusion across the cell membranes of white matter fiber tracts, the local fiber direction is believed to be inferable from the preferred free water diffusion direction, as derived from diffusion coefficients measured locally along several directions. The local fiber orientation in each image voxel is often described in terms of a diffusion tensor which is constructed from the measured diffusion coefficients. The diffusion tensor field forms the basis for determining the course of entire white matter fiber bundles.

Fiber tracking based on the diffusion tensor field is now used as a tool for noninvasive assessment of white matter fiber bundles (Conturo and others, 1999; Mori and others, 1999), both as a means of studying functional connectivity in the normal brain (Lehericy and others, 2004; Wakana and others, 2004) and in the study and treatment of disease (Dubey and others, 2005), for example, presurgical planning in cerebral neoplasms (Berman and others, 2004).

Due to the inherent noise of DTI measurements, the diffusion tensors are inaccurate, leading to possibly erroneous results in the derived white matter fiber paths. There is an extensive literature on different regularization techniques, see Anderson (2001), Wang and others (2003, 2004), McGrew and others (2004), Chen and Hsu (2005), and references therein. In particular, a sophisticated regularization technique based on local smoothing of the diffusion tensors has been suggested in Wang and others (2003). The latter paper also contains a review of DTI regularization techniques.

In the present paper, we develop a Bayesian framework for regularization of the tensor field. Prior information about the fiber architecture can thereby be exploited. Bayesian regularization techniques are well-known in the statistical literature, initiated in the 1980s by Geman and Geman (1984) and Besag (1986), and further developed in Besag and others (1991), Ripley (1991), and Titterington (1997). See also the recent review of Hurn and others (2003). Bayesian regularization of the primary directions of the diffusion tensors has been developed in Poupon and others (2000). Our paper is a natural continuation of this work. By regularizing the field of diffusion tensors, we fully utilize the 3-dimensional (3D) information on fiber orientation available in the tensors. It should be emphasized that the full tensor field contains important information other than the primary diffusion direction such as information about the local degree of anisotropy and crossing fiber bundles.

Our paper is organized as follows: In Section 2, we describe the field of diffusion tensors, more generally the field of diffusion functions. In Sections 3 and 4, the prior and likelihood models are discussed. The Bayesian regularization procedure is described in detail in Section 5. In Section 6, we study the performance of the regularization procedure on synthetic data, and in Section 7 we consider in vivo data. An extension of the prior model is described in Section 8, and perspectives are summarized in Section 9. Technical aspects of the regularization are collected in 3 appendices as supplementary material (available at Biostatistics online, http://www.biostatistics.oxfordjournals.org).

	2. THE FIELD OF DIFFUSION FUNCTIONS

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
Let be the finite set of voxels representing the white matter. For each voxel w, the true diffusion coefficient in a direction uS² on the unit sphere is denoted by f_w(u). Since the diffusion coefficient is the same in opposite directions, we have

We denote the field of diffusion functions by

It is believed that the diffusion coefficient f_w(u) is large if u is close to the main direction of the nerve fibers passing through w. One possibility for graphical representation of a diffusion function is the diffusion function profile

which is the spatial surface having the distance f_w(u) to the origin in the direction uS².

The diffusion function is commonly modeled as a quadratic form (cf. Basser and others, 1994). In this case,

(2.1)

where _w is a 3x3 positive semi-definite matrix, referred to as a tensor, and u^* is the transpose of u. The field of tensors will also be denoted by , for convenience. Let _{w1 w2 w3} 0 be the eigenvalues of _w with corresponding orthonormal eigenvectors u_w1,u_w2, and u_w3. In Figure 1, diffusion function profiles of quadratic diffusion functions (tensors) are shown. If _w1 > _w2, a diffusion function of the form (2.1) attains its maximal value _w1 in the unique unorientated direction u = ±u_w1. This direction is called the primary diffusion direction. In the DTI literature, the local variability in fiber orientation is often specified by the fractional anisotropy (FA) index

(2.2)

where

cf., for example, Basser and Pierpaoli (1996). Note that 0 FA_w 1. The lower limit is obtained if all eigenvalues are identical in which case the local fiber orientation distribution is isotropic. The upper limit is obtained if only the largest eigenvalue _w1 is different from 0, corresponding to perfectly aligned fibers at voxel w.

View larger version (21K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 1. Diffusion function profiles for quadratic diffusion functions (tensors) with eigenvalues 1 = 2 = 3 (left), 1 = 2 > 3 (middle), and 1 > 2 = 3 (right).

 
Figure 2 shows a 3D image of diffusion tensors. More precisely, we show the projection of an image consisting of 10x10x3 diffusion tensor profiles. The shape of an individual diffusion tensor profile is here typically like an ellipsoid or a "peanut shell." The colors are determined by the FA index. The actual values of the FA index can be looked up in the color bar to the right. To the left, the slice of the brain shows the location of the 3D tensor image. In the slice, the yellow fiber bundle with the shape of a horse shoe is the splenium of the corpus callosum which is the largest fiber bundle connecting the left and right hemispheres of the human brain.

View larger version (48K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 2. An example of a 3D diffusion tensor field of size 10x10x3, delineated by green. Each tensor is represented by its diffusion tensor profile, colored by its FA index, see the color bar to the right. To the left, the slice of the brain shows the position of the 3D tensor image. For more details, see the text.

 

	3. PRIOR MODEL

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
In order to perform Bayesian regularization, a prior model must be defined on the field of diffusion functions. In Bayesian analysis, Gibbs-type distributions (also sometimes referred to as Markov random fields) have been successfully applied as prior models for low-level tasks such as image restoration (see, e.g. Geman and Geman, 1984; Besag, 1986; Guttorp, 1995; Hurn and others, 2003). In the present paper, the same type of approach is suggested for the field of diffusion functions. Under the prior model, diffusion functions in neighbor voxels tend to be similar.

Let be a neighborhood relation on indicating the voxels from which local information is taken into account during regularization. In 3D, the common neighborhood structures are the 6, 26, or 32 nearest neighbors. In our application, we use the 26 nearest neighbors. A general Gibbs-type prior density of is

(3.1)

where Z is a normalizing constant and > 0. The summation involves all pairs of neighbor voxels. The function d(·,·) measures the distance between diffusion functions and ||w – w^'|| is the Euclidean distance between voxel w and w^'. In principle, any distance between diffusion functions can be used with the constraint that (3.1) should specify a proper prior, that is, the integral of the exponential term should be finite. Note that under the prior model, high probability fields have the property that diffusion functions in neighbor voxels are similar.

As the diffusion function contains directional as well as size information, the field of diffusion functions is commonly normalized before any further analysis (Basser and Pierpaoli, 1996). The function f_w is hence replaced by

(3.2)

Under the tensor model (2.1), we have

(3.3)

cf. Appendix A of the supplementary material (available at Biostatistics online), and (3.2) reduces to

In the following, we consider Gibbs-type prior densities on the normalized diffusion function field

(3.4)

Without specific assumptions on the form of the diffusion functions, we may choose

(3.5)

where g:[0,)[0,) is an increasing function with g(0) = 0. A simple choice would be

If the tensor model (2.1) holds, the integral of (3.5) reduces to

cf. Appendix A of the supplementary material (available at Biostatistics online). Here, the Frobenius norm

cf. Basser and Pierpaoli (1996), is used on the space of 3x3 matrices = {_ij}. Thus, (3.5) reduces to

(3.6)

where the function g in (3.6) has the same properties as in (3.5).

We now address the special case where only one eigenvalue is positive. Then,

(3.7)

where · indicates inner product, cf. Appendix A of the supplementary material (available at Biostatistics online). Here, u_w1·u_w^'1 is cosine of the angle between the 2 primary diffusion directions u_w1 and u_w^'1. In this case,

(3.8)

where g in (3.8) has the same properties as in (3.5). Poupon and others (2000) use a model similar to (3.8) as prior for the primary diffusion direction.

	4. LIKELIHOOD MODEL

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
We assume that the DTI data set consists of the following measured diffusion coefficients

Here, u₁,...,u_k are the directions in which the diffusion coefficients are measured. The minimum value of the number of directions k is 6. In the applications described in Sections 6 and 7, we used the algorithm described in Jones and others (1999) based on the principle of electrostatic repulsion to generate k evenly distributed directions.

The variability of a measured diffusion coefficient F_w(u) at voxel w in direction u depends on the true value of the diffusion coefficient f_w(u). Expressed more formally, we have

(4.1)

Note that (4.1) involves the usual conditional independence assumption. If we assume that the measured diffusion coefficients are independent and normally distributed, then the likelihood becomes

(4.2)

Because of detailed insight in the DTI scanning procedure, the likelihood model can be given theoretical support for a particular choice of h. In DTI, the diffusion coefficient F_w(u) in direction u is determined by the equation

where S_w0 is the signal intensity without any gradient, S_w(u) is the measured signal intensity, and b is the diffusion-encoding strength factor. The signal intensities S_w(u) and S_w0 follow Rice distributions, S_w(u)R(a,²) and S_w0R(a₀,²) (Sijbers and others, 1998). The signal-to-noise ratio (SNR) is defined for S_w(u) and S_w0 as SNR = a/ and SNR₀ = a₀/, respectively. The SNRs are known characteristics of the scanning procedure. If SNR is not too small (SNR 8), then according to Appendix B of the supplementary material (available at Biostatistics online) the approximation ²

(4.3)

holds, that is, the noise is Gaussian. Under the tensor model, we get

(4.4)

Note that the resulting model is scale invariant in the sense that, under (4.4), the likelihood (4.2) is proportional to the likelihood based on a scaled version of data

where . The scaled data follows a model of the type (4.4) with replaced by while is unchanged.

Recall from (3.3) that

and therefore can be unbiasedly estimated by

if the experimental directions u₁,...,u_k are uniformly distributed on S².

	5. THE BAYESIAN REGULARIZATION PROCEDURE

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
We have used the Metropolis–Hastings algorithm (see, for instance, Gilks and others, 1996) to simulate a regularized field of diffusion functions from the posterior distribution

We have worked under the tensor model in which case the likelihood p(F|) given by (4.2) is determined by

The prior p() takes the form, cf. (3.4) and (3.6),

We have used in the examples presented in Sections 6 and 7.

We start from some suitable initial field

The transition from the field _t at time t to the field _{t + 1} at time t + 1 is performed as follows:

1) Choose a normalized diffusion tensor for updating.

2) Sample a candidate normalized diffusion tensor from a proposal distribution q(·|_t), which may or may not depend on the current state _t of the chain.

3) Accept the candidate tensor with probability ß, where

Here, ^' is identical to _t except for the diffusion tensor at , which is .

4) If is accepted, then _{t + 1} = ^', else _{t + 1} = _t.

The choice of diffusion tensor to update may be done at random or according to a sequential ordering of the voxels. Another possibility is to use a random permutation scheme where each diffusion tensor is visited once during a sweep of the algorithm, but in random order (see, e.g. Guttorp, 1995, p. 55). In our implementation, we chose to visit the voxels in sequential order.

In theory, the proposal distribution q may have any form (for regularity conditions, see Roberts, 1996), but in practice it is important to select q carefully so that the Markov chain {_t} moves rapidly around the tensor space, yet has a reasonable proportion of candidate tensors accepted. This can be achieved by ensuring that the candidate tensor is close to the actual tensor of the chain, yet not so close that it takes many steps to move around in the tensor space.

The tensor space consists of 3x3 normalized positive definite matrices. The Wishart distribution generates positive definite matrices and a natural choice for the proposal distribution is therefore a normalized Wishart distribution. It is easy to simulate from the Wishart distribution, and flexibility is provided through the choice of degrees of freedom. Let X be distributed according to the Wishart distribution with mean and degrees of freedom n. We then propose to use the distribution of the normalized matrix as proposal distribution. In order to calculate the acceptance probability, we need the density of which is derived in Appendix C of the supplementary material (available at Biostatistics online).

Let _i,i = 1,2,3, denote the eigenvalues of the product . Then, we get from Appendix C

For the calculation of the posterior ratio p(^'|F)/p(_t|F), it is important to note that the chain only updates the diffusion tensor in one voxel at a time and therefore only terms in (3.4) and (4.2) involving this particular voxel contribute to the ratio.

In the direct method of simulating from the Wishart distribution, it is used that X = Z_iZ is if Z_i,i = 1,...,n(n 3), are independent, each distributed according to . For more efficient methods of simulating from the Wishart distribution than the direct one, see Liu (2001) and references therein. The value of the degrees of freedom n determines how close the candidate tensor is to the current tensor; as n increases the density of the proposal becomes more concentrated around the mean .

	6. SYNTHETIC DATA

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
To make an initial test of the regularization procedure, we analyzed the torus model presented in Tournier and others (2002). The fiber bundle is modeled as a horizontal torus of radius R and circular cross-section of radius r. The side length of a voxel is taken to be one unit length, and each voxel is represented by its center point.

At a voxel w inside the torus, a tensor _w is positioned with primary diffusion direction horizontal and perpendicular to the ray from O to w. The eigenvalues of _w satisfy

and are normalized such that and . The degree of anisotropy is the same for all tensors and is determined by = _w2/_w1. Since , we furthermore have

In terms of , the FA index (2.2) can be expressed as

If a voxel w is outside the torus, we let . There is thereby no preferred direction, = 1 and FA = 0. If the voxel is on the boundary of the torus, the diffusion tensors are determined as a volume-weighted average of the tensors outside and inside the torus. Finally, noise is added in all k directions for each voxel, according to the observation model (4.4) with .

The effect of regularization was assessed as follows. The normalized tensor field

was calculated and the noisy data

were simulated. Next, a least-squares method was used to fit a 3x3 positive definite tensor matrix _w(0) to the data in each voxel w. The corresponding normalized tensor field is denoted by

This tensor field was used as a starting value for the Markov chain Monte Carlo (MCMC) regularization algorithm described in Section 5. The voxels were visited in sequential order. In one sweep, each tensor was visited exactly once. The total number of tensors in the torus was 1335.

The regularized tensor field after s sweeps of the algorithm is denoted by

The average gain obtained by regularization was assessed by comparing the regularized tensor field after s sweeps to the true tensor field, that is, by considering the average tensor difference

as a function of the number s of sweeps.

Figure 3 shows the results obtained in an example where the parameters were chosen as R = 7.5, r = 3 (in voxel units), FA = 0.6, = 7.5, k = 17, n = 200, and SNR₀ = 25. The regularization is based on normalized tensors and therefore it does not depend on the value of b, as described in the discussion of scale invariance in Section 4. SNR₀ = 25 is the value obtained in the in vivo experiment. Tensors positioned outside the torus were not regularized.

View larger version (27K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 3. Left: Image of the projection of the primary diffusion direction and the corresponding fiber tracking before regularization. Right: Image of the projection of the primary diffusion direction and the corresponding fiber tracking after regularization. Tensors positioned outside the torus were not regularized.

 
In Figure 3, a central section of the torus is shown. At each voxel, the projection onto the section of the primary diffusion direction is shown before (left) and after (right) regularization. An example of fiber tracking in the raw and regularized fields of tensors is also shown in Figure 3, using an algorithm presented in Mori and others (1999). Regularization seemingly results in a more well-defined fiber orientation distribution. A plot of the average difference between the true diffusion tensors and the regularized tensors as a function of the number of sweeps through the tensor field showed convergence after approximately 400 sweeps.

The regularized image in Figure 3 (right) as well as those shown in Figures 6 (middle) and 7 (middle) represent one state of the Markov chain after convergence. We also determined estimates of posterior means based on 400 sweeps after convergence, and the result was very close to what is shown in Figure 3 (right). We advise future users of our regularization procedure to establish empirically the number of states needed in their application for obtaining a stable estimate of the posterior mean.

View larger version (18K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 6. Left: Tensors based on an individual data set before regularization. Middle: Tensors after regularization (400 sweeps through the full tensor field). Right: Average tensor field. Normalized tensors are shown. The colors indicate the value of the FA index. The same color coding as in Figure 2.

 

	7. IN VIVO DATA

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
All the scanning experiments were performed on a normal adult male volunteer in a 1.5-T GE Signa system (GE Medical Systems, Milwaukee, WI). The diffusion-weighted scanning was performed in k = 14 directions isotropically distributed in space. The diffusion-encoding strength was characterized by a b factor of 1000 bad hbox while SNR₀ = 25. The scanning resulted in 55 consecutive sections of the brain, each of thickness 2.5 mm. Each section consisted of 128 by 128 pixels, covering an area of 22 by 22 cm. The dimensions of a voxel are therefore 1.7x1.7x2.5 mm. In addition, 2 scannings without diffusion encoding (b = 0s/mm²) were acquired. The total scanning time was 306 s. This design was repeated 10 times providing 10 individual data sets. Finally, a T1-weighted sequence provided the high resolution structural image shown in Figure 4. During the data acquisition, the head of the volunteer was immobilized by a special foam pillow.

View larger version (56K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 4. Left: The position of the representative 10x10x9 cube '. Right: The projection of the 10x10 diffusion tensor profiles from the central horizontal slice of ' in one of the individual data sets. Normalized tensors are shown. The colors indicate the value of the FA index. The same color coding as in Figure 2 is used.

 
The observed diffusion coefficients will be denoted

where F_wij is the observed coefficient at voxel w, direction i, and replication j. In order to get an impression of the individual data available, we show in Figure 4 the central 10x10 slice of a representative 10x10x9 cube ^' and its position in the brain. In order to find an appropriate value of the prior parameter , we fitted a normalized tensor to each individual data set j at each voxel w. The averages

were regarded as the "true" tensors in a subsequent estimation of . The total prior tensor difference

(7.1)

is a sufficient statistic of the prior model which belongs to the family of exponential models. According to standard exponential family theory, the maximum likelihood estimate of is the value at which the mean of the sufficient statistic equals the observed value (Barndorff-Nielsen, 1978). The mean of the sufficient statistic as a function of was determined by MCMC simulation for voxels in the 10x10x9 cube ^'. The value of was estimated to be 7.5.

We applied our regularization procedure to each of the 10 individual data sets. All individual data sets gave similar results. In Figure 5, the average tensor difference

is shown for one of the data sets, as a function of the number s of sweeps through the full tensor field together with the corresponding plot for the voxels in ^'. In one sweep, each tensor was visited exactly once. The total number of voxels in the white matter was 156 000.

View larger version (8K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 5. Average tensor difference between tensors from an individual data set and average tensors. The upper curve is for all voxels in white matter (all WM voxels), while the lower curve is for voxels in the 10x10x9 subsample (selected WM voxels). The average tensor difference is shown as a function of the number s of sweeps through the full tensor field. In one sweep, each tensor is visited exactly once. The total number of voxels in white matter is 156 000.

 
In Figure 6, a central slice of tensors is shown before and after regularization together with the average tensor field. The overall impression is that the tensors after regularization resemble more closely the average tensors. The changes are moderate but important from a practical point of view. Analyzing averages of 2 data sets, we found that the precision of the regularized tensors based on a single scan corresponds to that of unregularized tensors based on 2 scans. The scanning time may therefore be reduced by a factor of 2 if the regularization procedure is used. This is important because a neuroscience laboratory has typically just one scanner but many computers suitable for performing the Bayesian regularization procedure. We also investigated the effect of regularization on the individual elements of the normalized tensors. As an example, Figure 7 shows the value of each of the diagonal elements of the normalized tensors in a full slice before and after regularization together with the value obtained for the average tensor field.

View larger version (38K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 7. The value of the ith diagonal element of the normalized tensors in slice 30 is shown in the ith row, i=1, 2, 3. The value before regularization is shown in the first column, the value after regularization in the second column, and the value from the average tensor field in the third column.

 
Prior to the analysis presented above for synthetic and in vivo data, we made a detailed investigation of the convergence properties of the MCMC procedure. In fact, we performed up to 40 000 sweeps through the full tensor field for one of the individual data set in the in vivo experiment to demonstrate the stability of the convergence which took 5–6 days on a PC equipped with a 2.8-Ghz Intel CPU.

Figure 8 shows for one of the individual data sets the computer time needed before the algorithm has converged, for different values of degrees of freedom in the Wishart distribution. For degrees of freedom in the range 200–300, the algorithm converges in less than 1 h, when the degrees of freedom is around 100 the algorithm converges in around 90 min, and for degrees of freedom less than 50 the algorithm needs more than 2 h for convergence. In the regularization presented in Figures 5–7, the number of degrees of freedom used was 200.

View larger version (10K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 8. Convergence properties of the MCMC algorithm. Minus log-posterior density is shown as a function of computing time (minutes) on a PC with a 2.8-GHz CPU for degrees of freedom in the Wishart distribution equal to n = 25,50,100,200,...,600.

 

	8. AN EXTENSION OF THE PRIOR MODEL

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
It has been stressed recently that the diffusion in voxels with, for example, crossing fibers can be quite complicated and not necessarily well characterized by the standard diffusion tensor model (2.1) (cf., e.g. Frank, 2001).

In Frank (2002) and Alexander and others (2002), a method of extending the tensor model to more complex configurations of fibers is presented. The diffusion function profile is modeled by truncating a spherical harmonic expansion at several orders. If the spherical harmonic expansion is truncated at order 0, an isotropic tensor is obtained, truncating at order 2 gives the tensor model, while higher-order truncations result in more complex diffusion function profiles. In order to describe multi-directional fibers and improve the regularization, such extensions are needed, see also the recent paper by Kreher and others (2005).

The prior model described in the present paper can be generalized to accommodate various complex fiber configurations. One possibility is to use the notion of a marked field of diffusion functions

Here, the mark m_w = 1 refers to a voxel with one predominant fiber direction, while m_w = 2 refers to a fiber crossing. If we let n(₁) be the number of voxels with mark 1, then a generalized prior model is

The parameter ß > 0 determines the fraction of voxels with crossing fibers. The distance function should depend on whether the voxels w and w^' both represent unidirectional fiber bundles or one of them is a fiber crossing location.

	9. PERSPECTIVES

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 
We have developed a framework for utilizing the full 3D information of diffusion tensor data in a regularizing step that allows incorporating prior information about the fiber structure. The actual choice of prior used in Sections 6 and 7 is not a canonical choice. We advise future users of our regularization procedure to experiment with the choice of the function g determining the prior.

For in vivo data, the regularization parameter was estimated by considering an average tensor field which is believed to have a negligible noise. In future applications of the Bayesian regularization procedure, it is of interest to develop a method of estimating from the noisy data. One possibility is to estimate as part of the Bayesian procedure. If p() is the prior distribution of , then the full posterior distribution becomes

where p(|) is given in (3.1) and p(F|) is given in (4.2). Simulation from the full posterior can be done as a Gibbs procedure alternating between updating and .

The full power of the Bayesian setup will be utilized if prior information from atlases of gross fiber architecture is incorporated in the analysis. Without such prior information, the simpler local smoothing techniques may show a similar performance (Glasbey and Horgan, 1995; Banham and Katsaggelos, 1997; Wang and others, 2003, 2004). It should be noted, however, that the Bayesian approach also delivers statements about the uncertainty of the regularization that may be fed forward to subsequent fiber tracking algorithms.

	ACKNOWLEDGMENTS

 
We would like to thank Dr Cyril Poupon, Prof. Carsten Gyldensted, and Dr Fabien Schneider for valuable discussions. The original manuscript was substantially improved by very helpful comments and suggestions from 2 anonymous reviewers. This work is supported by grants from the Danish National Research Foundation and the Danish Natural Science Research Council. Conflict of Interest: None declared.

	FOOTNOTES

 
We say that SR(a,²) if S has the same distribution as , where X and Y are independent random variables, XN(acos,²), and YN(asin,²). In the density of S, a modified Bessel function of the first kind and zero order appears.

	REFERENCES

TOP SUMMARY 1. INTRODUCTION 2. THE FIELD OF... 3. PRIOR MODEL 4. LIKELIHOOD MODEL 5. THE BAYESIAN REGULARIZATION... 6. SYNTHETIC DATA 7. IN VIVO DATA 8. AN EXTENSION OF... 9. PERSPECTIVES REFERENCES

 

Alexander DC, Barker GJ, Arridge SR. Detection and modeling of non-Gaussian apparent diffusion coefficient profiles in human brain data. Magnetic Resonance in Medicine (2002) 48:331–340.[CrossRef][Web of Science][Medline]

Anderson AW. Theoretical analysis of the effects of noise on diffusion tensor imaging. Magnetic Resonance in Medicine (2001) 46:1174–1188.[CrossRef][Web of Science][Medline]

Banham M, Katsaggelos A. Digital image restoration. IEEE Signal Processing Magazine (1997) 14:24–41.

Barndorff-Nielsen O. Information and Exponential Families in Statistical Theory (1978) Chichester: Wiley.

Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophysical Journal (1994) 66:259–267.[Web of Science][Medline]

Basser PJ, Pierpaoli C. Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. Journal of Magnetic Resonance Series B (1996) 111:209–219.[CrossRef][Web of Science][Medline]

Berman JI, Berger MS, Mukherjee P, Henry RG. Diffusion-tensor imaging-guided tracking of fibers of the pyramidal tract combined with intraoperative cortical stimulation mapping in patients with gliomas. Journal of Neurosurgery (2004) 101:66–72.[Web of Science][Medline]

Besag JE. Statistical analysis of dirty pictures. Journal of the Royal Statistical Society Series B (1986) 48:259–302.

Besag JE, York J, Mollie A. Bayesian image restoration with two applications in spatial statistics (with discussion). Annals of the Institute of Statistical Mathematics (1991) 43:1–59.[CrossRef][Web of Science]

Chen B, Hsu EW. Noise removal in magnetic resonance diffusion tensor imaging. Magnetic Resonance in Medicine (2005) 54:393–407.[CrossRef][Web of Science][Medline]

Conturo TE, Lori NF, Cull TS, Akbudak E, Snyder AZ, Shimony JS, McKinstry RC, Burton H, Raichle ME. Tracking neuronal fiber pathways in the living human brain. Proceedings of the National Academy of Sciences USA (1999) 31:10422–10427.

Dubey P, Fatemi A, Huang H, Nagae-Poetscher L, Wakana S, Barker PB, van Zijl PC, Moser HW, Mori S, Raymond GV. Diffusion tensor-based imaging reveals occult abnormalities in adrenomyeloneuropathy. Annals of Neurology (2005) 58:758–766.[CrossRef][Web of Science][Medline]

Frank LR. Anisotropy in high angular resolution diffusion-weighted MRI. Magnetic Resonance in Medicine (2001) 45:935–939.[CrossRef][Web of Science][Medline]

Frank LR. Characterization of anisotropy in high angular resolution diffusion-weighted MRI. Magnetic Resonance in Medicine (2002) 47:1083–1099.[CrossRef][Web of Science][Medline]

Geman S, Geman D. Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Transactions on Pattern analysis and Machine Intelligence (1984) 6:721–741.[Web of Science]

Gilks WR, Richardson S, Spiegelhalter DJ. Introducing Markov chain Monte Carlo. In: Markov Chain Monte Carlo in Practice—Gilks WR, Richardson S, Spiegelhalter DJ, eds. (1996) London: Chapman and Hall. 1–19.

Glasbey CA, Horgan GW. Image Analysis for the Biological Sciences (1995) Chichester: Wiley.

Guttorp P. Stochastic Modeling of Scientific Data (1995) London: Chapman and Hall.

Hurn M, Husby O, Rue H. Advances in Bayesian image analysis. In: Highly Structured Stochastic Systems—Green P, Hjorth NL, Richardson S, eds. (2003) Oxford: Oxford University Press. 302–322.

Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magnetic Resonance in Medicine (1999) 42:515–525.[CrossRef][Web of Science][Medline]

Kreher BW, Schneider JF, Mader I, Martin E, Hennig J, Il'yasov KA. Multitensor approach for analysis and tracking of complex fiber configurations. Magnetic Resonance in Medicine (2005) 54:1216–1225.[CrossRef][Web of Science][Medline]

Lehericy S, Ducros M, Van de Moortele PF, Francois C, Thivard L, Poupon C, Swindale N, Ugurbil K, Kim DS. Diffusion tensor fiber tracking shows distinct corticostriatal circuits in humans. Annals of Neurology (2004) 55:522–529.[CrossRef][Web of Science][Medline]

Liu JS. Monte Carlo Strategies in Scientific Computing (2001) New York: Springer.

McGrew T, Vemuri BC, Chen Y, Rao M, Mareci T. DT-MRI denoising and neuronal fiber tracking. Medical Image Analysis (2004) 8:95–111.[CrossRef][Web of Science][Medline]

Mori S, Crain BJ, Chacko VP, van Zijl PC. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Annals of Neurology (1999) 45:265–269.[CrossRef][Web of Science][Medline]

Pajevic S, Basser PJ. Parametric and non-parametric statistical analysis of DT-MRI data. Journal of Magnetic Resonance (2003) 161:1–14.[CrossRef][Web of Science][Medline]

Poupon C, Clark CA, Frouin V, Régis J, Bloch I, Le Bihan D, Mangin J-F. Regularization of diffusion-based directional maps for the tracking of brain white matter fascicles. NeuroImage (2000) 12:184–195.[CrossRef][Web of Science][Medline]

Ripley B. The use of spatial models as image priors. In: Spatial Statistics and Imaging—Possolo A, ed. (1991) Volume 20. Hayward, CA: Institute of Mathematical Statistics. 309–340. IMS Lecture Notes.

Roberts GO. Markov chain concepts related to sampling algorithms. In: Markov Chain Monte Carlo in Practice—Gilks WR, Richardson S, Spiegelhalter DJ, eds. (1996) London: Chapman and Hall. 45–54.

Sijbers J, den Dekker AJ, Verhoye M, Van Audekerke J, Van Dyck D. Estimation of noise from magnitude MR images. Magnetic Resonance in Imaging (1998) 16:87–90.[CrossRef]

Titterington DM. Some aspects of Bayesian image analysis. In: The Art and Science of Bayesian Image Analysis—Mardia KV, Gill CA, Aykroyd RG, eds. (1997) Leeds: Leeds University Press. 153–160.

Tournier JD, Calamante F, King MD, Gadian DG, Connelly A. Limitations and requirements of diffusion tensor fiber tracking: an assessment using simulations. Magnetic Resonance in Medicine (2002) 47:701–708.[CrossRef][Web of Science][Medline]

Wakana S, Jiang H, Nagae-Poetscher LM, van Zijl PC, Mori S. Fiber tract-based atlas of human white matter anatomy. Radiology (2004) 230:77–87.[Abstract/Free Full Text]

Wang Z, Vemuri BC, Chen Y, Mareci T. A constrained variational principle for direct estimation and smoothing of the diffusion tensor field from DWI. Information Processing in Medical Imaging (2003) 18:660–671.[Medline]

Wang Z, Vemuri BC, Chen Y, Mareci T. A constrained variational principle for direct estimation and smoothing of the diffusion tensor field from complex DWI. IEEE Transactions on Medical Imaging (2004) 23:930–939.[CrossRef][Web of Science][Medline]

Received March 31, 2006; revised September 22, 2006; revised November 24, 2006; revised January 12, 2007; accepted for publication February 7, 2007.


CiteULike    Connotea    Del.icio.us    What's this?





This Article Abstract FREE Full Text (PDF) Supplementary Material All Versions of this Article: 8/4/784    most recent kxm005v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Request Permissions Disclaimer Google Scholar Articles by Frandsen, J. Articles by Vedel Jensen, E. B. Search for Related Content PubMed PubMed Citation Articles by Frandsen, J. Articles by Vedel Jensen, E. B. Social Bookmarking          What's this?

Online ISSN 1468-4357 - Print ISSN 1465-4644

Copyright © 2009 Oxford University Press

Oxford Journals Oxford University Press

• Site Map
• Privacy Policy
• Frequently Asked Questions

Other Oxford University Press sites: Oxford University Press Oxford Journals China Oxford Journals Japan American National Biography Booksellers' Information Service Children's Fiction and Poetry Children's Reference Corporate & Special Sales Dictionaries Dictionary of National Biography Digital Reference English Language Teaching Higher Education Textbooks Humanities International Education Unit Law Medicine Music Online Products Oxford English Dictionary Reference Rights and Permissions Science School Books Social Sciences Very Short Introductions World's Classics