Department of Signals and Systems Division of Signal Processing and Biomedical Engineering CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2014 Non-Gaussian Diffusion Imaging: A comparison of the bi-Gaussian and the fourth-order tensor models to diffusion tensor imaging Master’s Thesis in Biomedical Engineering SOFIA GUNNARSSON Abstract Diffusion weighted imaging (DWI) is a magnetic resonance imaging (MRI) technique that can be used to measure, in vivo, the self-diffusion of water molecules in body tissues. This reveals information about the microstructure of the underlying tissue. The simplest approach to modelling the diffusion at each voxel, known as diffusion tensor imaging (DTI), is as a Gaussian probability density function (PDF). The covariance matrix of this PDF defines a second-order tensor. Whilst this model is adequate for modelling fascicles with a single dominant orientation it is not suitable for more complex tissue architectures; e.g. crossing fascicles. This has driven the development of more sophisticated models such as high-order tensors (HOTs). The aim in fitting such models is to determine: (i) the number fascicles and their orientations and (ii) to extract scalar measures, such as fractional anisotropy (FA), to characterise the tissue microstructure. This knowledge is used for two general types of applications: tractography (visualisation/identification of neural fibre tracts) and the detection/characterisation of diseased tissue. This thesis explores two non-Gaussian diffusion-modelling approaches. In particu- lar the bi-Gaussian (multi-tensor model with two tensors) and the fourth-order tensor are compared to conventional DTI. In the first part of the thesis synthetic data were generated for a range of signal-to-noise (SNR) values to simulate, for a single voxel, one fibre (for a variety of FA-values), and two crossing fibres (for six FA combinations and a range of crossing angles). Each model was fitted to this data and the quality of the fit evaluated in terms of the three measures: signal deviation, angle deviation and tensor element deviation (in the latter case only when the model matched the struc- ture). The results show that: (i) the fourth-order tensor model consistently has lower signal deviation than the bi-Gaussian and second-order tensor models; (ii) for one fibre both the second-order tensor and the fourth-order tensor models yield low angular devi- ation, whilst the bi-Gaussian model does not (model-structure mismatch); (iii) for two fibres the angular deviation decreases for both the bi-Gaussian model and the fourth order model with increasing SNR and crossing angle; (iv) the tensor element deviation decreases with increasing SNR. In the second part of the thesis the concept of FA is explored. In the case of the conventional second-order tensor a standard definition exists for computing a scalar mea- sure of FA from the tensor. FA-definitions for the bi-Gaussian and fourth-order tensor models are explored and compared, empirically (using both the synthetic data and real brain DWI data), to the second-order definition. FA estimation for high-order tensors involves the concept of Z-eigenvalue decomposition. The thesis presents an implemen- tation that is more robust and accurate than the existing open-source implementation. The experimental results for the synthetic data show that: (i) the FA-definition for the second-order tensor in general yields FA-values close to true FA for one fibre and close to the mean value of true FA for two fibres; (ii) the bi-Gaussian FA-definitions all tend to overestimate the true FA for both one and two fibres; (iii) the fourth-order definitions yield FA-values that correlate well with true FA for one fibre and for two fibres with a small crossing angle, but decreases with increasing crossing angle. The results for the real brain data show that FA-maps estimated with the bi-Gaussian model (fitted with the Levenberg-Marquardt algorithm) yields many undefined values, but that FA-maps estimated with the fourth-order tensor definitions correlate well with the maps estimated with the second order tensor model. Keywords: Fractional Anisotropy, bi-Gaussian diffusion model, fourth-order diffusion tensor, Diffusion Tensor Imaging, Diffusion Weighted Imaging, Diffusion Weighted MRI Acknowledgements I would like to give my appreciation to Mohammad Alipoor for his infinite patience when answering my many questions and his valuable inputs and to Andrew Mehnert for his engaged feedback. For providing answers to my medical queries and for valuable discussions I would like to thank Ylva Lilja. I would like to thank MedTech West for providing an inspirational environment and the means to carry trough this master thesis. Finally I would like to thank my parents, Rune and Helene Gunnarsson for their endless support and love, and Daniel Gleeson for always making time for me, discussing with me and making me believe he is as interested in my work as I am. The Author, Gothenburg, September 25, 2014 Contents 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Diffusion Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 Scalars derived from DWI data . . . . . . . . . . . . . . . . . . . . 2 1.2 Aims and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Structure of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 DWI 5 2.1 What is Diffusion, and what is DWI? . . . . . . . . . . . . . . . . . . . . . 5 2.2 How to Measure Diffusion with MRI . . . . . . . . . . . . . . . . . . . . . 5 2.2.1 The b−value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Gaussian Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3.1 DTI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 Non-Gaussian Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.1 Multi-Gaussian Diffusion . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.2 High Order Diffusion Tensor . . . . . . . . . . . . . . . . . . . . . 9 3 Empirical Evaluation of Non-Gaussian Diffusion Models 10 3.1 Evaluation Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.1.1 Fascicle structures/architecture defined in this study . . . . . . . . 11 3.1.2 Gradient Directions . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.3 Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.4 Procedure for creating rotations of the crossing fibre signals . . . . 13 3.2 Fitting the Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.1 DTI (Second order tensor model) . . . . . . . . . . . . . . . . . . . 14 3.2.2 Bi-Gaussian Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.3 Fourth Order Tensor Model . . . . . . . . . . . . . . . . . . . . . . 16 3.3 Measures of Goodness-of-Fit . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.3.1 Signal Deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 i CONTENTS 3.3.2 Angle Deviation and Estimation of the Main Direction . . . . . . . 17 3.3.3 Tensor Element Deviation . . . . . . . . . . . . . . . . . . . . . . . 18 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.4.1 Results for One Fibre-Structure . . . . . . . . . . . . . . . . . . . . 19 3.4.2 Results for Two Fibre-Structure . . . . . . . . . . . . . . . . . . . 23 3.5 Discussion of Quality of Fit Results . . . . . . . . . . . . . . . . . . . . . . 33 4 FA - A Measure of Structure 34 4.1 What is FA? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.1.1 Properties of FA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.1.2 What to expect from FA . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2 Generating Signals to Estimate FA from . . . . . . . . . . . . . . . . . . . 36 4.3 Calculating FA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.1 FA for DTI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.2 FA for the Bi-Gaussian Model . . . . . . . . . . . . . . . . . . . . 36 4.3.3 FA for the Fourth Order Model . . . . . . . . . . . . . . . . . . . . 37 4.4 Z-Eigenvalue Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.4.1 Theory of Z-eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . 38 4.4.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.4.3 The Improved Implementation of Z-eigenvalues Decomposition . . 41 4.5 Results of FA-estimations . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.5.1 Results for one Fibre FA-Estimations . . . . . . . . . . . . . . . . 42 4.5.2 Results for two Fibre FA-Estimations . . . . . . . . . . . . . . . . 44 4.5.3 Results of Real Data FA-Estimations . . . . . . . . . . . . . . . . . 47 4.6 Conclusions, Discussion and Improvement of FA-estimations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5 Summary and Conclusions 52 5.1 Thesis Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.2 Key Contributions and Findings . . . . . . . . . . . . . . . . . . . . . . . 53 5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.5 Further Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Bibliography 58 A Plots of Quality Measures 59 A.1 One Fibre Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 A.1.1 Angle Deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 A.1.2 Signal Deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 A.1.3 Tensor Element Deviation . . . . . . . . . . . . . . . . . . . . . . . 63 A.2 Two Fibre Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 A.2.1 Angle Deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 A.2.2 Signal Deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 ii CONTENTS A.2.3 Tensor Element Deviation . . . . . . . . . . . . . . . . . . . . . . . 88 B Plots of FA-values 94 B.1 One Fibre Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 B.2 Two Fibre Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 B.2.1 High-High FA Fibre Structure . . . . . . . . . . . . . . . . . . . . 97 B.2.2 High-Medium FA Fibre Structure . . . . . . . . . . . . . . . . . . . 100 B.2.3 High-Low FA Fibre Structure . . . . . . . . . . . . . . . . . . . . . 103 B.2.4 Medium-Medium FA Fibre Structure . . . . . . . . . . . . . . . . . 106 B.2.5 Medium-Low FA Fibre Structure . . . . . . . . . . . . . . . . . . . 109 B.2.6 Low-Low FA Fibre Structure . . . . . . . . . . . . . . . . . . . . . 112 C FA Maps from Real DWI-Data 115 iii Abbreviations 2DT Second order Diffusion Tensor DTI Diffusion Tensor Imaging DW Diffusion Weighted DWI Diffusion Weighted Imaging FA Fractional Anisotropy Fascicle A bundle of fibres FoO Fraction of Occupancy HOT High Order Tensor In vivo Latin for within the living MRI Magnetic Resonance Imaging PDF Probability Density Function ROI Region of Interest SNR Signal to Noise Ratio Voxel Volume element of a 3D digital image; analogous to a picture element (pixel) of a 2D digital image iv 1 Introduction T he research described in this thesis was undertaken in the MedTech West centre located at Sahlgrenska University Hospital in Gothenburg, Sweden in 2013. It constitutes part of a larger research project concerned with the devel- opment of diffusion imaging tools for the assessment of the visual pathways in the human brain. This chapter starts with a brief introduction to DWI, and then presents the aim and objective of this thesis, the scope, and finally a description of the structure of this thesis. 1.1 Background An understanding of how the human nervous system, particularly in the brain, is or- ganized has been a long standing area of research. Historically only post mortem and invasive analysis was possible, but now with diffusion weighted MRI (Magnetic Reso- nance Imaging) (hereinafter denoted DWI) it is possible to visualise the white matter tracts in vivo [1]. With region of interest (ROI) methods it is possible to separate and identify indi- vidual tracts. The construction of a 3D-model of the white matter tracts was one of the first applications of DWI. Comparisons of the tractography maps derived from DWI, with structures identified using more invasive methods have shown that DWI is a reliable method for relatively large and homogeneous structures such as the corpus callosum [1]. 1.1.1 Diffusion Models The easiest way to model the structure is with DTI (Diffusion Tensor Imaging) which assumes that the diffusion (see section 2.1) PDF can be described with a second order tensor. When a voxel contains crossing, kissing or other multiple fibre directions this model is inappropriate because the second-order tensor describes only an ellipsoid of diffusion with one principal direction. In DWI the resolution is low; one voxel in clinical 1 1.1. BACKGROUND DWI typically has a side of one or two millimetres. Since the likelihood that all fibres within one voxel have the same alignment is low, a more complete model is desirable. DTI utilises a second order tensor to describe the diffusion, and models the sig- nal attenuation as a mono-exponential. However the observed attenuation in tissue is not mono-exponential [2]. Nevertheless, modelling the signal with DTI is still a useful technique and can provide information about the structure of the tissue, with the large advantage of being relatively fast to compute. Numerous other models have been proposed, including the bi-Gaussian model [2], where the sum of two mono-exponentials model the signal; the multi-Gaussian, similar to the bi-Gaussian but instead a sum of any number of mono-exponentials; DSI, diffu- sion spectrum imaging; and the high order tensor model, where the conventional second order tensor is replaced with a higher order tensor. These models seek to more accu- rately describe the hindered diffusion of the water molecules, and thus to yield better reconstructions of the underlying structure. 1.1.2 Scalars derived from DWI data From the data acquired from a DWI-scan and then reconstructed with DTI, it is possible to derive a scalar called FA (Fractional Anisotropy). This is often referred to as a measure of structure or integrity (as discussed in [3]). The possibility to compute a scalar that describes the fibre integrity or fibre structure of each voxel could help determine the progress of diseases like MS (Multiple Sclerosis) or Alzheimer’s Disease that cause demyelination of white matter fibres [4]. It is also proposed that FA-maps can provide a method to demonstrate improvement of white matter integrity after tumour removal. In those cases where the micro-structural integrity is of interest (demyelination, struc- tural damage, etc.), the fibre directions are unimportant [3]. FA is used to help determine a suitable starting point for tractography, because a high FA also is indicative of the structure in the voxel. A high degree of structure is in this context a voxel with highly aligned white matter fibres. However it can also be argued that a voxel with two clearly distinguishable fascicles is highly structured. Similarly a voxel containing three, four or more fascicles are highly structured when compared with, for example, grey matter. It is clear that structure and integrity are difficult terms to define. Therefore, clinical experts engaged in research in DWI have been consulted throughout this project, to gain an understanding of what it is that they perceive FA to mean. A standard definition of FA exists for DTI, but not for other reconstruction models. In this thesis we evaluate the various definitions that have been proposed in the literature for a select set of models as outlined in the Aims and Objectives. For models where it is possible to generate multiple FA values, for example the multi-tensor models, the DWI community is not in agreement about whether FA should be expressed as one or several scalars. In this thesis we have adopted the convention that FA is a single scalar when generating FA maps. 2 1.2. AIMS AND OBJECTIVES 1.2 Aims and Objectives The first aim of this research was to compare the conventional Gaussian diffusion model with two non-Gaussian diffusion models, the bi-Gaussian and the fourth-order tensor, for single and crossing fibre architectures. The second aim of this research was to explore the concept of FA and to evaluate existing FA definitions for the three models chosen in the first aim. To this end the research had the following objectives: i) To generate synthetic data to simulate a single voxel, a single fibre (for three different FA values), and two crossing fibres (for six FA combinations, and a range of crossing angles). ii) To fit each model to this data and evaluate the quality of the fit in terms of the three measures: signal deviation, angle deviation and tensor element deviation. iii) To review the physical meaning and significance of FA. iv) To develop a robust and accurate implementation of the Z-eigenvalue decomposition algorithm needed for the FA estimation for high order tensors. v) To investigate and compare FA definitions for the bi-Gaussian and fourth-order tensor models, empirically (using both the earlier synthetic data and real brain DWI data), to the second-order definition. 1.3 Scope The non-Gaussian diffusion models explored in this research were limited to the bi- Gaussian and fourth-order tensor. The general problem of estimating the number of fascicles and their relative fractions within a voxel was not considered in this thesis. 1.4 Structure of the Thesis The remainder of the thesis is organised as follows: Chapter 2 Briefly reviews the concept of diffusion, describes what DWI (Diffusion Weighted Imaging) is and how it is used to measure diffusion, and presents the three diffusion models explored in this thesis. Chapter 3 Addresses objectives i) and ii) associated with the first aim of the thesis. In particular it evaluates how well the three diffusion models perform in terms of signal deviation, angle deviation and tensor element deviation. The framework used for synthe- sising the signal is presented and how the quality measures are calculated is described. The resulting measures are presented and the models compared. 3 1.4. STRUCTURE OF THE THESIS Chapter 4 Addresses objectives iii), iv) and v) associated with the second aim of the thesis. In particular the concept of FA and what it means is discussed. The protocol for synthesising signals with known FA for use in the experiments is described, as well as how to estimate the FA values for the three different models. Z-eigenvalue decomposition of high order tensors is discussed, and a new implementation is presented. Finally the results of the FA value estimation experiments for the different models are presented and compared. Chapter 5 Summarises the thesis and the conclusions. Future work and improvements are proposed. Appendix The complete set of results are presented here, in three appendices. The first with the results from model evaluation and the second from synthetic FA-estimations and the third with real data FA-maps. 4 2 Diffusion Weighted Imaging T his chapter begins by describing diffusion, and the difference between iso- tropic and anisotropic diffusion. Section 2.2 then describes diffusion weighted imaging (DWI) which is an magnetic resonance imaging (MRI) technique for measuring bulk diffusion in vivo. Section 2.3 describes conventional Gaussian diffusion (the simplest model of diffusion). Finally Section 2.4 describes the two non- Gaussian diffusion models explored in this thesis. 2.1 What is Diffusion, and what is DWI? The term diffusion refers to the random, or Brownian, motion of particles or molecules in a liquid or gas. DWI is an MRI technique that measures the self diffusion of water molecules in body tissue, i.e. the random displacement of the molecules due to thermal agitation. In an unrestricted volume, the diffusion of water is isotropic, the Brownian motion is the same in all directions. However within structured tissue, like the brain white matter, the diffusion is not isotropic, but depends on the direction. Anisotropic diffusion is not restricted to brain white matter, but can be found in other biological tissues also [5]. Herein though the focus is on white matter. DWI can be used to perform tractography, that is to identify and visualise the white matter fibre tracts via the main direction of diffusion in each voxel. DWI can also provide information about the progress of a disease or condition, for example a brain tumour, via both tractography and scalar measures (features) describing tissue structure. 2.2 How to Measure Diffusion with MRI In DWI the signal derives from the signal attenuation due to diffusion or more precisely the proton spins dephasing in a spatially varying magnetic gradient [5]. A scan is made 5 2.2. HOW TO MEASURE DIFFUSION WITH MRI without diffusion weighting (b = 0 1), providing a reference signal, S0. Then DWI mea- surements are executed with a diffusion sensitising gradient applied in several directions (b 6= 0). A model of diffusion is then fitted voxel-wise to this data. The model makes it possible to quantitatively characterise the diffusion. The complexity of the model together with sampling considerations determines the minimum number of gradient di- rections. The b-value is determined by several factors including the desired resolution, the apparatus and the scanning time. DWI is an MRI method where the applied magnetic field is not homogeneous as in T1- and T2-weighted images. Instead the pulsed magnetic field is linearly varied causing the protons in the water to precess at different rates. A refocusing gradient pulse of the same magnitude but opposite direction is applied, and those protons that have moved will not be perfectly refocused, causing a detectable signal loss. This diffusion weighted method was introduced by Stejskal and Tanner, who gave name to the Stejskal-Tanner- equation: S S0 = e−TE/T2e−γ 2G2δ2(∆− δ 3)D (2.1) describing the signal attenuation, where S is the signal with applied gradient G and S0 without, γ is the gyromagnetic ratio, δ the duration of the pulses, ∆ the time from the first pulse to the refocusing pulse, TE the echo time, T2 the spin-spin relaxation time, and D the diffusion-coefficient [5]. To simplify equation (2.1) it is convenient to set b = γ2G2δ2 ( ∆− δ 3 ) , this is called the b-value. Additionally, if the term e−TE/T2 is assumed to be constant over all gradients, the Stejskal-Tanner equation is reduced to: S S0 = e−bD. (2.2) In DWI, the conventional way to estimate the diffusion is to model the ADC (Ap- parent Diffusion Coefficient) with a second order tensor according to ADC = − log ( S S0 ) . (2.3) This has known disadvantages when the underlying structure does not have one main direction of diffusion, but consists of several directed sources of diffusion. Many other approaches have been proposed to address this problem. The bi-Gaussian model and the fourth-order diffusion tensor are two such models. In this thesis they are compared to the second order tensor used in conventional DTI (Diffusion Tensor Imaging) in terms of signal error, angular error and tensor element error. The tensor element error only makes sense when the structure and model matches, for example for DTI when modelling one fibre, and bi-Gaussian when modelling two crossing fibres. 1b = γ2G2δ2 ( ∆ − δ 3 ) usually in s/m2. For further explanation of the b-value see the following paragraphs and Section 2.2.1 6 2.3. GAUSSIAN DIFFUSION 2.2.1 The b−value As described previously, the term b-value was introduced in [6] to simplify the Stejskal- Tanner equation, (2.1). This is somewhat misleadingly suggesting that the b-value is less complex than it is. Commonly in literature the b-value is described as a constant value, e.g. 1500 as used in this report, but b = γ2G2δ2 ( ∆− δ 3 ) . These constituent parameters all have different impacts on the result. The term e−TE/T2 in equation (2.1) is often omitted since it generally is considered to be constant over all gradients. Nevertheless it describes the signal’s exponential dependency on TE . Increasing TE will decrease the signal amplitude, and significantly decrease the SNR (Signal to Noise Ratio) for all measurements regardless of the b-value. Considering only unit-norm gradients, bnominal = γ2δ2 ( ∆− δ 3 ) , the relationship between TE and b-value can be approximated according to TE ≈ 12bnominal γ2 . (2.4) Increasing the b-value, will also result in a larger minimal TE and thus decrease the SNR [7]. This is one of the important trade-off’s in DTI, a higher b-value is known to increase the contrast between the DW (Diffusion Weighted) gradient directions, but forces TE to rise and SNR to drop. If we instead not are considering unit-norm gradients, the b-value can be increased without affecting the bnominal. Since bk = bnominalG 2 k the desired gain can be obtained by increasing the gradients to have norm greater than one. 2.3 Gaussian Diffusion Before DWI became a widely-used imaging technique, it had been observed that in dif- fusion MRI the echo intensity was dependant on the specimen’s orientation with respect to the direction of the applied magnetic field gradient, characterised by the diagonal elements of the effective self-diffusion tensor Deff. Deff is a 3 times 3 symmetric ma- trix with six unique elements. In 1993, Basser, Mattiello and LeBihan in [6] proposed that the off-diagonal elements of Deff must be taken into account to fully characterise anisotropic diffusion. In other words, both the eigenvectors and eigenvalues are of impor- tance; not only the eigenvalues. This provided a more general model than the previous, where the diagonal elements non-equality only provided information about the isotropic or anisotropic nature of the structure. An estimation of the fibre tracts orientations and distributions in vivo, i.e. tractography, was now possible. 2.3.1 DTI In DTI the PDF (Probability Density Function) of the water molecules diffusion is assumed to be Gaussian. The covariance matrix of this PDF is a second order tensor 7 2.4. NON-GAUSSIAN DIFFUSION (2DT) D written D =  Dxx Dxy Dxz Dxy Dyy Dyz Dxz Dyz Dzz  . (2.5) The tensor is symmetric with 3 rows and 3 columns, and thus out of those nine elements only six are unique. This covariance matrix describes the diffusion along and correlation between orthogonal axes [5]. Using the Stejskal-Tanner equation the signal with gradient directions ~g (which is included in b) is S = S0 · e−b:D (2.6) where b : D is defined according to b : D = b · (gxxDxx + gyyDyy + gzzDzz + 2gxyDxy + 2gxzDxz + 2gyzDyz) . (2.7) To determine D at least seven measurements must be performed where one usually is with b = 0. The b = 0 signal is here denoted S0. Since DTI utilises a second order tensor to describe the diffusion, it cannot describe more than one main direction (corresponding to the dominant eigenvalue/-vector pair). For complex fibre structures such as crossing, kissing, branching, etc., this disadvantage motivates the use of more complex models. 2.4 Non-Gaussian Diffusion The first proof that the diffusion of water-molecules in the brain tissue is not Gaussian was presented in the 1990s[1]. Two methods to model the non-Gaussian properties are described here. The first method, the multi-Gaussian model, extends Equation (2.6) to a sum of exponentials. This has been discussed and evaluated in [1], [2], [7]. The second method, the high order diffusion tensor model, can model more complex inter-voxel fibre structures [8]. This model extends the second order tensor to a fourth order tensor, and can thus describe more complex structures. 2.4.1 Multi-Gaussian Diffusion The bi-Gaussian model is an extension of the mono-exponential conventional second order diffusion tensor. A second exponential term is added, thus the signal now is described by S = S0(A1 · e−bg:D1 +A2 · e−bg:D2). (2.8) Here A1 and A2 are the fractions of occupancy, i.e. what proportion of the voxel’s volume that is assigned to that tensor (fascicle). This can be extended to a multi- Gaussian model, with any number of exponentials, each corresponding to one fascicle’s direction. 8 2.4. NON-GAUSSIAN DIFFUSION The simplification utilised in DTI for the relationship between S0 and S, transforming the non-linear problem to a linear problem cannot be used here. Estimating the two tensors D1 and D2 in Equation (2.8) is more difficult than to estimate D in Equation (2.6). Additionally, to estimate the scalars A1 and A2, it was demonstrated by [7] that the signals must be acquired over a range of b−values to avoid an infinite number of solutions. 2.4.2 High Order Diffusion Tensor In 2002, high angular resolution diffusion imaging (HARDI) was proposed in [9], to de- scribe non-Gaussian diffusion. Over a sphere N discrete gradient directions are sampled. Without assuming the nature of the diffusion, the apparent diffusion coefficient is esti- mated in all gradient directions. One of the methods proposed to analyse HARDI data is with high order tensors (HOT’s). An advantage of using the high order diffusion tensor is that the 2DT theory can be generalised. The signal is modelled with Equation (2.6) but D is now a 3-dimensional fourth order tensor, instead of a 3-dimensional second order tensor. 9 3 Empirical Evaluation of Non-Gaussian Diffusion Models T his chapter presents an empirical evaluation of the two selected non-Gaus- sian diffusion models – the bi-Gaussian and the 4th-order tensor – and the conventional 2nd-order tensor model used in DTI (Diffusion Tensor Imaging). The aim of the evaluation was to investigate how well the two non-Gaussian models describe single and crossing fibres compared to DTI. Several experiments are presented that evaluate the performance of the three models in terms of signal devia- tion, angle deviation and tensor element deviation. In the first section the framework used to generate the synthetic data (with prescribed structure, properties, and noise) is described. Section 3.2 describes how the three models were fitted to the synthetic data. In the first section the framework used to generate the synthetic data (with prescribed structure, properties, and noise) is described. Section 3.3 defines the three quality of fit measures and the required adjustments depending on the fibre architecture/structure and model. Section 3.4 presents the fitting results. A discussion of these results is then presented in Section 3.5. 3.1 Evaluation Framework The commonly used evaluation framework of [10] was used in this study. This framework prescribes the steps involved in synthesising signals corresponding to different rotations of a fascicle architecture in a single voxel (for different noise levels), model fitting, and measurement of the goodness of fit. More specifically the framework comprises the following steps: (i) The structure to synthesise is defined, i.e. the number of fibres and their FA-values (FA is short for fractional anisotropy). 10 3.1. EVALUATION FRAMEWORK (ii) The gradient directions and b-value are defined. (iii) The number of rotations NR and Euler angles, noise levels Nσ, number of noise realisations NRN and crossing angles Nφ are defined. (iv) The tensor defined in step (i) is rotated with an angle defined in step (iii). (v) The diffusion signal is simulated, for underlying structure with one fibre according to the Stejskal-Tanner equation (2.1), or for two fibres according to Equation (2.8). (vi) The tensors D is estimated with each model. (vii) The measures of goodness of fit and derived scalar measures are calculated. (viii) Steps (v) to (vii) are repeated for a total of NRN Rician noise realisations. (ix) Steps (iv) to (viii) are repeated for a total of NR different rotations. (x) The mean values and standard deviations for the goodness of fit measures and scalar measures are calculated. In this research the focus lies in how well the FA-maps correlate with the underlying structure. Thus the FoO (Fraction of Occupancy) values have been assumed to be known during the simulations and estimations, and were always set to 0.5 each, such that there are two fascicles of the same fraction in each voxel. 3.1.1 Fascicle structures/architecture defined in this study The fascicle structures/architecture used in this study were all defined in terms of the number of crossing fibres (one or two) and their individual FA-values. Three different FA-values were used: FAhigh = 0.94, FAmedium = 0.51 and FAlow = 0.18. For one fibre this resulted in three FA-structures: high, medium and low, and for two crossing fibres in six: high-high, high-medium, high-low, medium-medium, medium-low and low-low. The three different FA-values defined three diffusion tensors, one high, one medium 11 3.1. EVALUATION FRAMEWORK and one low according to the following Dhigh =  17 0 0 0 1.01 0 0 0 1  · 10−4 (3.1) (3.2) Dmedium =  17 0 0 0 10 0 0 0 5  · 10−4 (3.3) (3.4) Dlow =  1.4 0 0 0 1.1 0 0 0 1  · 10−4. (3.5) Their magnitudes are based on the largest diffusivity 1 experimentally found along a fibre [3]. Each was used as an initialisation tensor Dinit to generate the synthetic signals and were appropriately rotated to create the desired fibre-crossing structure. Due to symmetry, it suffices to define the acute angle of intersection between two fascicles to be from 0◦ to 90◦ in steps of 10◦. 3.1.2 Gradient Directions The gradient directions used in the study are the 81 evenly distributed directions defined in [11]. This was chosen to assure that the tensor-element estimation is robust. For DTI it was shown in [12] that at least 30 unique sampling directions are necessary to estimate the six unique elements. For the fourth order tensor, as described in Section 3.2.3, we have fifteen unique elements to estimate requiring a larger set of sampling directions. The b-value was chosen to be 1500 s/mm2 which is a value that can be found in the clinical settings [7]. 3.1.3 Noise The noise in measured MRI (Magnetic Resonance Imaging) signals is commonly assumed to be of Rician distribution [13]. Therefore we added Rician noise to the synthetically generated signals for values of σ ∈ [0.02, 0.14] in steps of 0.02, where σ = 1/SNR. This resulted in NRN = 7 noise realisations with SNR ∈ [7.1, 50]. In [3] a lowest SNR for practical use and a given b-value was suggested according to SNRmin = 3 e−b·10−3 (3.6) 1Diffusivity has unit m2/s and describes how fast the molecules diffuse 12 3.1. EVALUATION FRAMEWORK where b is in mm2/s. With a b-value of 1500 mm2/s this means that in practical use the SNR for that b-value should be larger than 14. Given that when we generate the synthetic signal we can choose what SNR-range to use, it was chosen to cover this smallest value with margin, thus covering a range from more noise than preferred to very little noise. The Rician noise was generated and added to the noise-free signal S according to Snoise = √ (S + n1)2 + (n2)2 (3.7) where Snoise is the noisy signal, ni = σδi and δ ∈ [0,1] is a normally distributed pseudo- random number. For each noise level the procedure was repeated NRN = 25 times each time for new δ. In this way seven noise level sets of synthetic signals were generated, giving a total of 25× 7 noise-realisations. 3.1.4 Procedure for creating rotations of the crossing fibre signals Since the number of gradient directions is finite, the fibre orientation in relation to the gradient directions has an impact on the signal intensity. For example, in the extreme case with only one gradient direction applied, a fibre perpendicular to the gradient direction generates a stronger signal than a parallel one would, see equation (2.6) and (2.7). To limit the signal-dependence on the fibre-gradient-orientation relation, multiple rotations of each crossing fibre signal were generated. Each diffusion tensor D was rotated with a rotation matrix R to obtain D = RDinitR T, where R was defined by combining the three rotational matrices accord- ing to equation (3.8), one for rotation about each axis, Rx, Ry and Rz. The three rotation matrices are shown in Equations (3.9) – (3.11) where α, β and γ are the angles to rotate around the x-, y- and z-axis respectively. Dinit is an initialisation tensor with elements defined to achieve a desired FA-value, as described in Section 3.1.1. R = RxRyRz (3.8) Rx =  1 0 0 0 cos(α) − sin(α) 0 sin(α) cos(α)  (3.9) Ry =  cos(β) 0 sin(β) 0 1 0 − sin(β) 0 cos(β)  (3.10) Rz =  cos(γ) − sin(γ) 0 sin(γ) cos(γ) 0 0 0 1  (3.11) The rotation procedure was performed in the same way for both one fibre and two fibre synthetic signals. Rotations corresponding to the following intervals α ∈ [0, π/2], 13 3.2. FITTING THE MODELS β ∈ [0, π/2] and γ ∈ [0, 2π] were considered. For α and β the intervals were evenly divided into three steps each and γ for four. Thus the fibres were rotated a total amount of 36 times. 3.2 Fitting the Models The Gaussian diffusion model (DTI) and the two non-Gaussian models (the fourth order tensor model and the bi-Gaussian model) were fitted to each synthetic signal in turn as described below. 3.2.1 DTI (Second order tensor model) The estimate of the tensor D, denoted D̂, was estimated using a linear least squares fit of the diffusivity function d(~gi) with respect to the six unique elements in D̂ according to min D̂ { N∑ i=1 ( d(~gi)− d̂(~gi) )2 } . (3.12) The relation between the diffusivity function and the logarithm of the signal is linear, viz. d(~gi) = −1 b log ( Si S0 ) (3.13) and the relation between the measured signal and estimated tensor is Ŝi = S0e −b~giD̂~gTi . (3.14) 3.2.2 Bi-Gaussian Model The bi-Gaussian model was fitted using non-linear least squares (Levenberg-Marquardt algorithm) to obtain estimates for the two second order tensors D1 and D2 according to min D̂1,D̂2 { N∑ i=1 ( Si − Ŝi )2 } (3.15) where N is the number of gradient directions. The relationship between the measured signal and estimated tensors is given by Ŝi = 0.5e−b~giD̂1~gTi + 0.5e−b~giD̂2~gTi (3.16) which cannot be expressed as a linear least squares problem as was the case for DTI. To initialise the minimisation, three approaches were investigated: 1. Ranged random initialisation [9], where the initial guess is two second order tensors of three dimensions, whose elements are of 10−4 order of magnitude. The order of magnitude is based on the free diffusion of water at 37◦ [7] [6]. 14 3.2. FITTING THE MODELS (a) Initial guess: Ranged ran- dom (b) Initial guess: 2DT (c) Initial guess: 2DT randomly perturbed Figure 3.1: Angle deviation as a function of fibre crossing-angle for the bi-Gaussian tensor estimation of different initialisation methods. (a) Initial guess: Ranged ran- dom (b) Initial guess: 2DT (c) Initial guess: 2DT randomly perturbed Figure 3.2: Signal deviation as a function of fibre crossing-angle for the bi-Gaussian tensor estimation of different initialisation methods. 2. A conventional 2nd-order tensor is estimated, which is used as the initial guess for both bi-Gaussian tensors [14]. 3. To avoid local minima, a conventional 2nd-order tensor is estimated and randomly perturbed (a random number in the range [−1, 1] × 10−4 is added to each tensor element) as initial guess. This is repeated twenty-five times (no improvement was detected with more than 25 repetitions). The angular deviation for the three initialisation methods is shown in Figure 3.1. For angles lower than 50◦ the results are similar, but for larger angles the two random initialisations are superior. With the DTI initialisation the optimisation step found local minima for large crossing angles, which resulted in a larger angular deviation. The signal deviation, see Figure 3.2, shows the same pattern. For low angles the methods’ deviations are similar, but for higher angles the DTI initialisation yields larger deviation. The computation times for method 1 and 2 are significantly lower than for method 3. 15 3.3. MEASURES OF GOODNESS-OF-FIT 3.2.3 Fourth Order Tensor Model One method for indexing the elements of a (fourth order) tensor A is by four indices i1, . . . i4, one for each order, im ∈ [1,2,3] for three dimensions. But the diffusion tensor has an inherent symmetry by definition (similar to the symmetry in DTI) which makes the fourth order, three dimensional tensor possible to represent with fifteen unique elements. Each can be represented by a lower case letter with two indices, here akl where k,l ∈ 0, . . . ,4 and effectively determine the number of 1’s and 2’s (and therefore also 3’s) for Ai1i2i3i4 . One way to interpret the indices for A is to translate 1,2 and 3 to x,y and z respectively. Example The element A1123 of the supersymmetric tensor A can also be represented by the element a21 in the fifteen element long vector a. For the fourth order model the fitting introduced in [15] that ensures that the fourth order tensor is PSD (Positive Semi-Definite) was used. The estimated signal was then fitted, using non-linear least squares, to the Cartesian tensor orientation distribution function (CT-ODF) according to min λ { N∑ i=1 ( Si S0 − Ŝi )2 } (3.17) and the estimated signal is defined according to Ŝi = ∫ S2 M∑ j=1 λjp(~gi;~c) 2B(~v,~g, b)d~v (3.18) where B is the response function and the kernel is a fourth order tensor defined by λ. The author of [15] has developed an open source software package which was used to perform the fourth order tensor fit. 3.3 Measures of Goodness-of-Fit For each synthetic signal, the three models were fitted and three goodness-of-fit measures calculated for each. The calculation of these measurements depends not only on the underlying structure (represented by the synthetic signal), but also on the model used. This section describes how the calculations of the measures were performed. 3.3.1 Signal Deviation Signal deviation ∆S provides a measure of the percentage deviation between the mea- sured signal S and the signal Ŝ estimated by the model; i.e. ∆S = 1 N N∑ i Si − Ŝi Si . (3.19) 16 3.3. MEASURES OF GOODNESS-OF-FIT Ŝ was calculated according to Equation (3.14) for DTI, Equation (3.15) for the bi- Gaussian model and Equation (3.18) for the fourth order tensor model. 3.3.2 Angle Deviation and Estimation of the Main Direction The angle deviation measure ∆φ gives the angular difference between the simulated structure’s direction or directions and the estimated direction. In [16] a mean angle error was used for validation of the results. Since DTI yields one main direction, bi- Gaussian two and the fourth order tensor theoretically up to thirteen directions the angular deviation was defined in three different ways. The number of fibres were assumed to be known when estimating the angle deviation for the fourth order tensor model; i.e. the problem of estimating how many crossing fascicles are contained in a single synthetic signal was not considered in this thesis (as noted in Section 3.1). Angle Deviation for Single Fibre Simulations When simulating one single fibre, the underlying structure can be described with one clearly defined main direction φ. For the simplest case, DTI, the angle deviation was defined as the absolute difference in degrees between the estimated second order tensor’s main direction and the simulated fibre’s direction. The eigenvector ~emax corresponding to the largest eigenvalue λmax is directed in the estimated fibre’s main direction r̂ such that ~emax = r̂. The angle deviation for DTI was defined according to ∆φ = ∣∣∣φ− φ̂∣∣∣ . (3.20) The angle between two vectors was determined according to∣∣∣φ− φ̂∣∣∣ = arccos ( r · r̂√ r2 · r̂2 ) (3.21) where r is the ground truth main direction and equal to the main eigenvalue’s correspond- ing eigenvector. Since arccos(x) ∈ [0, 180] and the acute angle was sought, φ = 180− φ if the resulting angle was larger than 90◦. For the bi-Gaussian model, the two tensors both gave one estimated main direction each, φ̂1 and φ̂2. The angle deviation was calculated according to ∆φ = 1 2 2∑ i=1 ∣∣∣φ− φ̂i∣∣∣ (3.22) where the angular difference and directions were calculated in the same way as for DTI, see Equation (3.21). As previously mentioned, the estimation of the number of fibres was not within the scope of this thesis. Thus so for the fourth order tensor model one main direction was sought. This was, in the same way as for the conventional second order tensors, found 17 3.3. MEASURES OF GOODNESS-OF-FIT by the maximum eigenvalue and it’s corresponding eigenvector. The generalisation of eigenvalues and eigenvectors of two-dimensional matrices into symmetric matrices of even order are called Z-eigenvalues and Z-eigenvectors and is described in Section 4.4. The angle deviation was then calculated in the same way as for DTI, according to Equation (3.20). Angle Deviation for Simulations with two Crossing Fibres When the underlying structure was two simulated crossing fibres with main directions r1 and r2 the angle deviation for DTI was defined according to ∆φ = 1 2 2∑ i=1 ∣∣∣φi − φ̂∣∣∣ (3.23) where |φ− φ̂| was defined in the same way as for one simulated fibre, see Equation (3.21). For the bi-Gaussian and fourth order tensor models the angle deviation for the two estimated directions r̂1 and r̂2 had to be be combined with the right ground truth directions r1 and r2 which was achieved by minimising the sum of the possible differences as follows ∆φ = 1 2 min {∣∣∣φ1 − φ̂1 ∣∣∣+ ∣∣∣φ2 − φ̂2 ∣∣∣ , ∣∣∣φ1 − φ̂2 ∣∣∣+ ∣∣∣φ2 − φ̂1 ∣∣∣} . (3.24) 3.3.3 Tensor Element Deviation The tensor element deviation ∆D was only calculated when the model and underlying structure matched, i.e. when estimating one simulated fibre with DTI or two simulated fibres with the bi-Gaussian model. It is defined according to ∆D = 1 n n∑ i=1 ∣∣∣Di − D̂i ∣∣∣ (3.25) where n is the number of tensor components in the model. 18 3.4. RESULTS 3.4 Results For all setups the three measures of quality were plotted as a function of SNR and in the case of two fibres, the crossing angle as well. The complete set of figures can be found in Appendix A. A representative selection of results is presented below. 3.4.1 Results for One Fibre-Structure 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 Angle Deviation, High FA, FA = 0.94 SNR M e a n A n g le D e v ia ti o n i n d e g re e s DTI Bi−Gaussian Model Fourth Order Model Figure 3.3: Mean angle deviation, one fibre. Ground truth FA= 0.94. The bi-Gaussian model’s large angle deviation is because the model approximates two tensors when the underlying structure is composed of only one. A summary of the number of times that the two non-Gaussian models yielded smaller deviation than DTI, and that the fourth order tensor model yielded smaller deviation than the bi-Gaussian is presented in Table 3.2 (note that the percentage was calculated over all SNR-levels). This shows that the fourth order model yielded smaller angular deviation than the bi-Gaussian for all FA-values, and compared to DTI the deviation is smaller for the high FA-value but not medium and low. The bi-Gaussian model had the largest deviation for all three FA-values. Figure 3.4 presents the angle deviations for each component tensor of the bi-Gaussian model separately. For high and medium FA the best fit-deviation is lower than 20◦ for all SNR-levels and lower than 10◦ for SNR above 14. Compared to the results in Figure 3.3 this is low. So the bi-Gaussian model seems to, in the case of medium or high FA-values, fit one of the two tensors well, but the second one poorly. In Figure 3.5 the mean signal deviation for the low FA-value and all models is shown which indicates that all models give a similar deviation which decreases with increasing SNR. In table 3.2 the models are compared, which shows that DTI gives the smallest signal deviation in a majority of the estimations for one fibre simulations. 19 3.4. RESULTS Table 3.1: The percentage of times the bi-Gaussian and the fourth order tensor models had lower angle deviation than DTI, and that the fourth order model had smaller deviation than the bi-Gaussian model (H = high FA, M = medium FA, and L = low FA). H M L Bi-Gaussian 1% 2% 13% Fourth Order 68% 41% 26% Fourth Order better than bi-Gaussian 99% 95% 81% Table 3.2: The percentage of times the bi-Gaussian and the fourth order tensor models had lower signal deviation than DTI, and that the fourth order model had smaller deviation than the bi-Gaussian model (H = high FA, M = medium FA, and L = low FA). H M L Bi-Gaussian 23% 7% 20% Fourth Order 11% 4% 19% Fourth Order better than bi-Gaussian 26% 18% 34% The tensor element deviation for DTI is shown in Figure 3.6. The magnitude of the tensor element deviation is of the same order of magnitude as the ground truth tensor elements, and the smallest deviation was obtained for low FA. 20 3.4. RESULTS 0 10 20 30 40 50 0 10 20 30 40 50 60 SNR A n g le D e v ia ti o n i n d e g re e s Minimum Angle Deviation, bi−Gaussian model 0 10 20 30 40 50 40 45 50 55 60 65 SNR A n g le D e v ia ti o n i n d e g re e s Maximum Angle Deviation, bi−Gaussian model High FA Low FA Medium FA High FA Low FA Medium FA Figure 3.4: Angle deviation for the two main directions of the two tensors estimated with bi-Gaussian model, to the left the direction with lowest angle deviation, to the right the largest deviation. 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 Signal Deviation, Low FA, FA = 0.18 SNR M e a n S ig n a l D e v ia ti o n i n % DTI Bi−Gaussian Model Fourth Order Model Figure 3.5: Mean signal deviation, one fibre. Ground truth FA= 0.18 21 3.4. RESULTS 5 10 15 20 25 30 35 40 45 50 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x 10 −4 Tensor Element Deviation DTI, One Fibre, All FA SNR M e d ia n T e n s o r E le m e n t D e v ia ti o n Low FA Medium FA High FA Figure 3.6: Median tensor element deviation, one fibre, DTI. Ground truth FAlow = 0.18, FAmedium = 0.51 and FAhigh = 0.94. 22 3.4. RESULTS 3.4.2 Results for Two Fibre-Structure Angle Deviation 0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25 30 35 40 45 50 DTI, Angle Deviation Angle M e a n A n g le D e v ia ti o n SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 3.7: Angle deviation as function of fibre crossing angle. Ground truth FA1 = 0.94 and FA2 = 0.94. In Figure 3.7 a typical result for the angle deviation with DTI is shown. Note that the deviation is computed as the average of the deviations of the direction estimated by the model and the two ground truth directions (see Equation (3.23)). The deviation was shown to increase almost linearly with increasing crossing angle, showing a precise but not accurate behaviour 2. The estimated direction was close to the direction directly in-between the two simulated fibre directions (see Figure 3.8). All simulated FA-values showed the same trends but with decreasing mean-FA the dependence of SNR increased (see Figure 3.9 where both ground truth FA-values are medium and low, FA1 = 0.51 and FA2 = 0.18). The bi-Gaussian model (note that this is the model used to synthesise the signal) performed best when both fibres had high FA and high crossing angle. When the crossing angle is low, the structure approaches the behaviour of one fibre in which case the bi- Gaussian model was shown to give high angular deviation (see Section 3.4.1). For the structure of two high-FA fibres Figure 3.10 shows the angular deviation as a function of the crossing angle. For the fourth order tensor model, as described earlier, the number of fibres and their fraction of occupancy were not estimated but assumed to be known. In the same manner as for the bi-Gaussian model, when the angle decreases and approaches parallel fibres, the structure cannot be distinguished from one fibre. So for parallel fibres the model 2Precision means that the standard deviation is small but says little about how accurate the result is, while accuracy means that the mean value is a good estimate but nothing about the spread of the results (the deviation). 23 3.4. RESULTS Figure 3.8: The dashed line represents the estimated direction with DTI when the under- lying structure, represented by the two filled lines, is two crossing fibres 0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25 30 35 40 45 50 DTI, Angle Deviation Angle M e a n A n g le D e v ia ti o n SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 3.9: DTI, angle deviation as function of fibre crossing angle. Ground truth FA1 = 0.51 and FA2 = 0.18. and the structure does not match. As can be seen in Figure 3.11 the angle deviation was small for crossing angles over 70◦ and large for smaller angles. For all simulations the deviation was larger for smaller crossing angles and vice versa, but with decreasing FA this change is less distinct. A comparison between the angular deviations for all of the models is presented in Table 3.3 where all crossing angles were taken into account. It can be seen that the bi-Gaussian model only gave a majority of lower angular deviations when both fibres had high FA-values, and the fourth order model never had the majority of lowest angu- lar deviations. However, since these models do not estimate the number of fibres and for low crossing angles this is a problem (as shown by the results in Section 3.4.1), the same comparison was done but only for crossing angles larger than 60◦. The results are 24 3.4. RESULTS 0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25 30 Bi−Gaussian, Angle Deviation Angle M e a n A n g le D e v ia ti o n SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 3.10: Bi-Gaussian model, angle deviation as function of fibre crossing angle. Ground truth FA1 = 0.94 and FA2 = 0.94. 0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25 30 35 40 45 50 Fourth Order, Angle Deviation Angle M e a n A n g le D e v ia ti o n SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (a) Fourth order model, angle deviation as function of fibre crossing angle. 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50 Fourth Order, Angle Deviation SNR M e a n A n g le D e v ia ti o n i n d e g re e s Φ = 0 Φ = 10 Φ = 20 Φ = 30 Φ = 40 Φ = 50 Φ = 60 Φ = 70 Φ = 80 Φ = 90 (b) Fourth order model, angle deviation as function of SNR. Figure 3.11: Angle deviation as function of fibre crossing angle for fourth order tensor estimation. Ground truth FA1 = 0.94 and FA2 = 0.94. presented in Table 3.4. For those angles the fourth order model gave lower angular devi- ation for all FA-values, and the bi-Gaussian model as well except when both FA-values are low. Additionally the fourth order model gave, in a majority of the measurements, a lower angular deviation than the bi-Gaussian model. Signal Deviation For DTI the signal deviation decreased with increasing SNR, but as can be seen in Figure 3.12 crossing angle has a large impact. This is expected due to the model-structure mismatch. The crossing angle had most impact when both fibres had high FA-value. 25 3.4. RESULTS Table 3.3: The percentage of times (over all SNR levels and crossing angles) the bi-Gaussian and the fourth order tensor model had a lower angle deviation than DTI, and that the fourth order model had a lower deviation than the bi-Gaussian model. LL, ML, etc. refers the ground truth FA value, Low-Low, Mean-Low and so on. LL ML MM HL HM HH Bi-Gaussian 23% 35% 38% 41% 46% 61% Fourth Order 44% 43% 38% 38 % 42% 32% Fourth Order better than bi-Gaussian 81% 62% 50% 46% 44% 29% Table 3.4: The percentage of times (over all SNR levels and crossing angles larger than or equal to 60◦) the bi-Gaussian and the fourth order tensor model had a lower angle deviation than DTI, and that the fourth order model had a lower deviation than the bi-Gaussian model. LL, ML, etc. refers the ground truth FA value, Low-Low, Mean-Low and so on. LL ML MM HL HM HH Bi-Gaussian 33 % 63% 70% 73% 80% 97% Fourth Order 72 % 88% 80% 85% 97% 81% Fourth Order better than bi-Gaussian 81% 77% 73% 66% 72% 55% When the SNR was high (how high depended on the FA-value, with higher FA the SNR- level was higher) the signal deviation increased with increasing crossing angle, but for low SNR the relation was reversed (this can be seen in Figure 3.13 where the SNR limit is 10). The bi-Gaussian model and fourth order model performed similarly for all com- binations of FA-values. The signal deviation decreased with increasing SNR, and for high-high, high-medium and medium-medium FA-structures, i.e. when no fibre had low FA-value, the signal deviation decreased with increasing crossing angle (see Figure 3.14). When instead one fibre or more had low FA, the crossing angle did not have any effect on the mean signal deviation (see Figure 3.15) and very little on the variance (see Figure 3.16). In Figure 3.15 and Figure 3.14 the bi-Gaussian model and fourth order model are similar but at least one difference can be noticed, i.e. the bi-Gaussian model has for SNR from 17 and above, a lower signal deviation. In the same way as for the angle deviation the total differences were summarised, and the bi-Gaussian and fourth order model were compared to DTI. The results for all angles can be seen in Table 3.5, and as can be seen DTI gives the lowest signal deviation 26 3.4. RESULTS 5 10 15 20 25 30 35 40 45 50 2 4 6 8 10 12 14 16 18 20 22 DTI Signal Deviation SNR M ea n S ig na l D ev ia tio n in % Φ = 0 Φ = 10 Φ = 20 Φ = 30 Φ = 40 Φ = 50 Φ = 60 Φ = 70 Φ = 80 Φ = 90 Figure 3.12: Signal deviation for DTI. Ground truth FA1 = 0.94 and FA2 = 0.94. 0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25 30 DTI Signal Deviation Angle M e a n S ig n a l D e v ia ti o n i n % SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 3.13: Signal deviation as function of fibre crossing angle for DTI. Ground truth FA1 = 0.94 and FA2 = 0.51. in a majority of the cases, with exceptions for high-high and high-low FA-values. In the last row it can be seen that the fourth order model had a lower signal deviation than the bi-Gaussian model in a majority of the estimations. For crossing angles larger than 60◦ see Table 3.6. Here we can see that in agreement with the plots for cases when at least one fibre had low FA, the results for one or more low FA-fibres in Table 3.6 and Table 3.5 are almost the same. For high-high structure, both the bi-Gaussian model and the fourth order model give lower signal deviation, but the high-medium structure still had lowest signal deviation with DTI. In both these tables it can be seen that the bi-Gaussian model in a majority of the estimations had 27 3.4. RESULTS 0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25 30 Bi−Gaussian Signal Deviation Angle M ea n S ig na l D ev ia tio n in % SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (a) Bi-Gaussian model, ranged random ini- tialisation, signal deviation as function of fibre crossing angle. 0 10 20 30 40 50 60 70 80 90 0 5 10 15 20 25 30 Fourth Order Signal Deviation Angle M e a n S ig n a l D e v ia ti o n i n % SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (b) Fourth order model, signal deviation as function of fibre crossing angle. Figure 3.14: Ground truth FA1 = 0.94 and FA2 = 0.94. 0 10 20 30 40 50 60 70 80 90 0 5 10 15 Bi−Gaussian Signal Deviation Angle M e a n S ig n a l D e v ia ti o n i n % SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (a) Bi-Gaussian model, ranged random ini- tialisation, signal deviation as function of fibre crossing angle. 0 10 20 30 40 50 60 70 80 90 2 4 6 8 10 12 14 16 Fourth Order Signal Deviation Angle M e a n S ig n a l D e v ia ti o n i n % SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (b) Fourth order model, signal deviation as function of fibre crossing angle. Figure 3.15: Ground truth FA1 = 0.94 and FA2 = 0.18. lower deviation than the fourth order model. Tensor Element Deviation For the synthetic signal with two fibres, the only tensor element deviation that could be calculated was when the model used was the bi-Gaussian. In all cases the tensor element deviation’s order of magnitude was reasonable, i.e. smaller than the ground truth tensor element. For all simulations except for the medium-medium FA-values, the deviation decreased with increasing SNR, e.g. see Figure 3.17. Similarly to the signal deviation, when one fibre had low FA-value, the deviation was not dependent on the crossing angle (see Figure 3.18). For high-high and high-medium FA, the deviation decreased with 28 3.4. RESULTS 0 20 40 60 80 7 8 10 13 17 25 50 SNR D e v ia ti o n i n p e rc e n t Bi−Gaussian Crossing angle = 0 0 20 40 60 80 7 8 10 13 17 25 50 SNR D e v ia ti o n i n p e rc e n t Fourth order Crossing angle = 0 0 20 40 60 80 7 8 10 13 17 25 50 SNR D e v ia ti o n i n p e rc e n t Bi−Gaussian Crossing angle = 90 0 20 40 60 80 7 8 10 13 17 25 50 SNR D e v ia ti o n i n p e rc e n t Fourth order Crossing angle = 90 Figure 3.16: Boxplots of mean of signal deviation for the bi-Gaussian and Fourth order model. Parallel and perpendicular fibres. Ground truth FA1 = 0.18 and FA2 = 0.18. Table 3.5: The percentage of times (over all SNR levels and crossing angles) the bi-Gaussian and the fourth order tensor model had a lower signal deviation than DTI, and that the fourth order model had a lower deviation than the bi-Gaussian model. LL, ML, etc. refers the ground truth FA value, Low-Low, Mean-Low and so on. LL ML MM HL HM HH Bi-Gaussian 20% 24% 9.6% 55% 29% 55% Fourth Order 20% 19% 5.2% 52% 24% 50% Fourth Order lower than bi-Gaussian 35% 40% 22% 33% 33% 30% increasing crossing angle. The Impact of SNR For the synthetic data with FA1 = FA2 = 0.51, SNR’s impact on the angle deviation depended on the simulated crossing angle (see Figure 3.19). For crossing angles above 29 3.4. RESULTS Table 3.6: The percentage of times (over all SNR levels and crossing angles larger than or equal to 60◦) the bi-Gaussian and the fourth order tensor model had a lower signal deviation than DTI, and that the fourth order model had a lower deviation than the bi-Gaussian model. LL, ML, etc. refers the ground truth FA value, Low-Low, Mean-Low and so on. LL ML MM HL HM HH Bi-Gaussian 20% 24% 10% 56% 42% 84% Fourth Order 20% 19% 6.8% 53% 45% 87% Fourth Order better than bi-Gaussian 35% 41% 24% 34% 41% 33% 70◦ the angle deviation decreases with increasing SNR, but for crossing angle lower than 60◦ the relationship is reversed. A similar behaviour occurred for FA1 = FA2 = 0.18, the low-low fibre configuration, where increasing the SNR increased the angle deviation (see Figure 3.20). 30 3.4. RESULTS 0 10 20 30 40 50 60 70 80 90 0 0.5 1 1.5 2 2.5 x 10 −3 Bi−GaussianTensor Element Deviation Angle M e d ia n T e n s o r E le m e n t D e v ia ti o n SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 3.17: Bi-Gaussian model, tensor element deviation as function of fibre crossing angle. Ranged random initialisation. High-high FA fibre structure. 0 10 20 30 40 50 60 70 80 90 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x 10 −3 Bi−GaussianTensor Element Deviation Angle M ed ia n T en so r E le m en t D ev ia tio n SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 3.18: Tensor element deviation for bi-Gaussian model estimation. Ground truth FA1 = 0.94 and FA2 = 0.18. Bi-Gaussian model, tensor element deviation as function of crossing angle. Ranged random initialisation. 31 3.4. RESULTS 0 10 20 30 40 50 60 70 80 90 5 10 15 20 25 30 35 40 45 Fourth Order, Angle Deviation Angle M e a n A n g le D e v ia ti o n SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 3.19: Angle deviation for fourth order tensor estimation. Ground truth FA1 = 0.51 and FA2 = 0.51. 5 10 15 20 25 30 35 40 45 50 50 52 54 56 58 60 62 64 66 68 Bi−Gaussian, Angle Deviation SNR M e a n A n g le D e v ia ti o n i n d e g re e s Φ = 0 Φ = 10 Φ = 20 Φ = 30 Φ = 40 Φ = 50 Φ = 60 Φ = 70 Φ = 80 Φ = 90 Figure 3.20: Random initialisation, angle deviation for bi-Gaussian model estimation. Ground truth FA2 = 0.18 and FA2 = 0.18. 32 3.5. DISCUSSION OF QUALITY OF FIT RESULTS 3.5 Discussion of Quality of Fit Results The angular deviation plots show that both the bi-Gaussian model and the fourth order model yielded more accurate direction estimates with increasing crossing angle. They also show that for small crossing angles the direction estimate is inaccurate. As previ- ously mentioned, the synthetic signal from two parallel fibres and one fibre cannot be distinguished, and the fourth order model had as low angular deviation as DTI for one fibre. This indicates that the fourth order model used here, needs to be complemented with an estimation of the number of fibres, and preferably their fraction of occupancy as well. For one fibre, the model-structure mismatch for the bi-Gaussian model meant that DTI generated lower deviation in both signal and angle. However the fourth order model had the lowest angular deviation for the high FA-simulations. The tensor element deviation could not be used to compare the different models, but was useful in evaluating the different initialisation methods for the bi-Gaussian model. It also makes sense in this case to compare relative values and not absolute, since some tensor elements could be very small and others big. The value of signal deviation as a quality measure is debatable. All models minimise with respect to the signal, and for SNR of 10 and above, they all result in the same range of deviation. When providing a value of how well your model performs, when the underlying structure is unknown, this measure is often the only available one, since you then do not know the ground truth directions of the fibre bundles, or the tensor element values, etc. The angle deviation is more informative. It provides information about how well the model actually reconstructs/estimates directions, and which is used (important in tractography). To do a more complete analysis of the fourth order model, estimation of the number of fibres should also be taken into consideration, as well as a larger variety of medium FA-values. Not only could the explicit FA values be within a range from e.g. [0.3, 0.7] but several different eigenvalue-relations that all generate the same FA could be used, to evaluate how that relation affects the resulting directions. The effect of noise should also be considered. It was shown in the simulations that crossing fibres with medium and low FA-values can generate a larger angular deviation with increasing SNR. In practice this could lead to the detection of false fibre tracts in tractography. 33 4 FA - A Measure of Structure T his chapter explores the concept of fractional anisotropy (FA). Its meaning and significance from the perspective of diffusion weighted imaging are dis- cussed in Section 4.1 and Section 4.1.2. Section 4.2 describes how signals with known FA were generated for the purpose of empirically evaluating FA esti- mation with the conventional second order tensor model (DTI), the bi-Gaussian model, and the fourth-order tensor model. Section 4.3 describes how FA is estimated for these three models. Section 4.4 presents the Z-eigenvalue decomposition method needed to estimate FA using the fourth-order tensor model, as well as a new implementation of the method developed by the author. Section 4.5 presents the results of the FA estimation experiments for the different models. Finally Section 4.6 discusses these results, draws conclusions, and offers suggestions for improved FA estimation. 4.1 What is FA? Fractional anisotropy, FA, is a scalar that for each voxel provides information about how anisotropic the diffusion within that voxel is. This scalar can then be interpreted in terms of fibre-integrity or structure. A standard definition of FA exists for DTI. Generalisations have been proposed for the bi-Gaussian and fourth order model, but there has not been a general method accepted throughout the DWI (Diffusion Weighted Imaging)-community. FA-maps can be used as a tool to find regions of interest; for example to locate the area in which a tractography should be started (the seed points). In that case the correlation between high FA and white matter-fibres is utilised. It can also be used to determine progression of a disease like MS (Multiple Sclerosis); the condition causes demyelination of the fibres which in turn decreases the FA-values which as shown in [17]. In the same way it could be used to confirm that a treatment has the desired effect, that the white matter integrity is improved which is indicated by increased FA-values. 34 4.1. WHAT IS FA? That FA does not consistently represent fibre integrity is known. In [3] it is pointed out that the decrease of FA in white matter-areas with crossing fibres is natural. In this section five FA-definitions for the non-Gaussian models are evaluated, including how the FA-value is altered with crossing fibres. 4.1.1 Properties of FA There are two properties of FA that must be fulfilled. FA ultimately is a measure of the anisotropy of diffusion within the voxel, which means two things. Firstly, that if the diffusion is isotropic, i.e. the Brownian motion is independent of direction, FA should be very low, preferably zero. Secondly, if the diffusion has one dominant direction, effectively if all white matter fibres are aligned, FA should be high, preferably one, since this is the most anisotropic environment. The first case, with isotropic diffusion, can also be regarded as the case where the white matter fibres are oriented in all (infinite) directions within the voxel. FA can take values in the range of [0, 1], so between the two extreme cases from one to infinite directions when the structure decreases, the FA-value must in some way decrease with increasing number of directions. 4.1.2 What to expect from FA FA is a scalar invariant; i.e. it is rotationally and translationally invariant. In order to evaluate whether a certain method gives a good FA value or not, a ’good FA-value’ must be defined. FA-values have two main applications, in tractography and in progression of a disease. For the tractography application the main issue is to distinguish between structured and unstructured tissue, such that the structured tissue has an FA-value that is sufficiently higher than the FA-value for unstructured tissue. In the disease progression application, the change over time is of interest. The FA-values from one session to another must make it possible to highlight the areas where the degree of myelination has changed since the last scan(s). So all increases or decreases of the FA-value are of importance. In order to specify the desired behaviour of FA, discussions have been held with Dr. Ylva Lilja MD and PhD-student at Sahlgrenska University Hospital, and with Mohammad Alipoor Lic. and PhD-student at Chalmers University of Technology. Y. Lilja’s research is within preoperative imaging of the optic nerves before neurosurgery with the aim to improve the knowledge of the anatomy and function of nerve tracts in the brain as a tool to aid in risk assessment before surgery. M. Alipoor research is within signal processing, he is involved in a project related to DWI analysis, with the aim to develop new and effective methods for tractography and the analysis of fascicles in the brain. Examples of their work include [18] and [10]. It was during these discussions concluded that FA should give information about the fibre integrity; e.g. about the demyelination of the white matter fibres. Even if the voxel of interest contains several crossing fibres, the FA should be able to provide information of the fibre integrity. However, when the number of fibre-direction increases and approaches 35 4.2. GENERATING SIGNALS TO ESTIMATE FA FROM infinity, the FA-values should decrease. In a voxel with an infinite number of direction the FA should approach zero. From one fibre-direction to infinite fibre-directions the FA should decrease from one to zero. Whether this change should be linear, exponential, or something else is an open question. 4.2 Generating Signals to Estimate FA from For the synthetic signal the same synthetic signal and structure as described in Section 3.1 was used. The real data was obtained at Sahlgrenska University Hospital, with one b = 0 acquisition, and a single-shell 1 with b = 800 s/mm2 and 32 gradient directions. An experimentally determined threshold of 20 % of maximum S0 was set to remove the voxels outside the brain. To find an informative and known structure, a region of interest-mask was drawn by a physician to locate a volume containing the optic chiasm 2. The slice that contained the selected voxels was then processed. 4.3 Calculating FA Prior to the FA-estimations the three models were fitted to the signal as described in Section 3.2. An additional method, not previously described, was also used to fit the fourth order model to the signal, described in Section 4.3.3. The FA definitions for the three models are described in Section 4.3.1, 4.3.2 and 4.3.3. 4.3.1 FA for DTI FA values for DTI were computed using the standard definition based on the eigenvalues of the estimated tensor, viz. FA = √ 1 2 (λ1 − λ2)2 + (λ1 − λ2)2 + (λ2 − λ3)2 λ2 1 + λ2 2 + λ2 3 (4.1) where λi is the ith eigenvector’s eigenvalue. By definition FA ∈ [0,1], where 0 reflects an isotropic voxel and 1 a highly anisotropic voxel. 4.3.2 FA for the Bi-Gaussian Model For the bi-Gaussian model two second order tensors D1 and D2 were estimated, from which two FA-values FA1 and FA2, one for each tensor, were calculated according to Equation (4.1). These two values were combined in the following three ways to obtain a single FA value: FAmean = mean (FA1,FA2), FAmax = max (FA1,FA2) and FAmin = min (FA1,FA2). 1The number of shells is determined by how many different b-values are used, often referred to as single-shell acquisition for one b-value and multi-shell acquisitions for more than one b-value. 2The optic chiasm is the area in the brain where the optic nerves, whom stretches from the eyes to the chiasm, partially cross each other and on the other side of the chiasm becomes the optic tracts. 36 4.4. Z-EIGENVALUE DECOMPOSITION 4.3.3 FA for the Fourth Order Model For the fourth order model the concept of Z-eigenvalue decomposition was used to define FA. Two methods were used to find the fourth order tensor, the first method is described in Section 3.2.3. The second method fits the fourth order tensor to the ADC (Apparent Diffusion Coefficient, see Equation (2.3))-profile with the diffusivity function d(g) = log ( S S0 ) = −b ∑3 i,j,k,l=1 (gigjgkglDijkl) . (4.2) This is essentially a generalisation of the method used to fit the tensor in DTI. The fourth order tensor was fitted to the ADC with least squares fit min D { N∑ i=1 ( d(gi)− d̂(gi) )2 } (4.3) where N is the number of gradient directions [13]. For the fourth order model, the conventional eigenvalue-decomposition cannot be used to obtain eigenvalues from which to define FA. Eigenvalues and eigenvectors are de- fined only for square matrices. However a generalisation, called Z-eigenvalue-decomposition was introduced in [13] for symmetric tensors of even order. This decomposition is dis- cussed in more detail in Section 4.4. The resulting Z-eigenvalues λiare the basis of two FA definitions found in the literature for the fourth order model. The first definition, introduced in [13], is FAQi = √ ν ν − 1 √∑ν i=1(λi − λmean)2∑ν i=1 λ 2 i (4.4) where λmean = 1 ν ∑ν i=1 λi. The second definition, introduced in [18] is FAMA = λmax νλmean . (4.5) Both definitions range from 0 to 1. 4.4 Z-Eigenvalue Decomposition In this section the author’s implementation of the Z-eigenvalue decomposition method is described. This implementation is based on previous implementations described in Section 4.4.2. The definitions and theory in this section are drawn from [19] and [13]. 37 4.4. Z-EIGENVALUE DECOMPOSITION 4.4.1 Theory of Z-eigenvalues Anmth order n-dimensional tensor Ai1,...,im is supersymmetric if Ai1,...,im = Aperm{i1,...,im} for all permutations of i1,...,im. A can be described by an mth order polynomial f(x) as follows f(x) ≡ Axm := ∑n i1,...,im=1 Ai1,...,imxi1 . . . xim (4.6) where x can be regarded as an mth order n-dimensional tensor of rank one. If the polynomial f(x) is positive semi-definite (PSD), then A is said to be PSD. The real solutions x and λ of { Ax = λx xTx = 1 (4.7) are called the Z-eigenvectors and Z-eigenvalues respectively. If and only if all Z-eigen- values of A are positive, A is PSD. If a second order tensor is symmetric, then it has only real eigenvalues with real eigenvectors. However this is not true for higher order super- symmetric tensors. Therefore it is important to emphasise that it is the real solutions that are Z-eigenvalues and Z-eigenvectors (the complete set, including complex solutions are denoted H-eigenvalues and H-eigenvectors). The symmetric hyper-determinant det (A) of the super-symmetric tensor A is defined with an irreducible polynomial in Ai1,...,im which vanishes when x ∈ Cn and x 6= 0 such that f(x) = 0 and∇f(x) = 0. If m = 2 this coincides with the conventional determinant. An extension of the Kronecker delta is defined according to δi1,...,im = { 1, if i1 = ... = im 0, otherwise (4.8) An mth order n-dimensional tensor is called the mth unit tensor if its elements are δi1,...,im for i1,...,im = 1 and is denoted I. Additionally, det (I) = 1 to specify sign and size of the symmetric hyper-determinant. If m is even, the Z-eigenvalues of A are the real solutions λ to φ(λ) = det (A− λI) . (4.9) φ is called the characteristic polynomial of A. The number of Z-eigenvalues is strictly less than m2 −m+ 1. Finding the Eigenvalues In this section, the method suggested in [13] for finding the Z-eigenvalues and Z-eigenvectors is summarised. 38 4.4. Z-EIGENVALUE DECOMPOSITION The Z-eigenvalues and eigenvectors are found by finding the real solutions to the equation system  4∑ i=1 4−i∑ j=0 idijx i−1 1 xj2x 4−i−j 3 = 4λx1 4∑ i=0 4−i∑ j=1 jdijx i 1x j−1 2 x4−i−j 3 = 4λx2 4∑ i=0 3−i∑ j=0 (4− i− j)dijxi1x j 2x 3−i−j 3 = 4λx3 x2 1 + x2 2 + x2 3 = 1 (4.10) Here dij denotes the 15 unique elements of a fourth order super-symmetric tensor D. The solutions can be divided into four cases described below, that are solved sepa- rately, and contribute to the set of eigenvalues. Case 1: The elements in the eigenvector x2 = x3 = 0 which correspond to d3,1 = d3,0 = 0 and thus x1 = ±1, λ = d4,0. Case 2: The elements in the eigenvector x1 = x3 = 0 which correspond to d1,3 = d0,3 = 0 and thus x2 = ±1, λ = d0,4. Case 3: If x3 = 0, and x1, x2 6= 0 then by eliminating λ (4.10) is simplified into the equation system  ∑4 i=1 idi,4−it i−1 = ∑3 i=0 (4− i) di,4−iti+1∑3 i=0 di,3−it i = 0 (4.11) where t = x1/x2. If the two equations have common solutions t, then x1 = t√ 1+t2 , x2 = ±1√ 1+t2 , x3 = 0 and λ = d(~x). Case 4: If x3 6= 0 then the resulting equation system is ∑4 i=1 ∑4−i j=0 idiju i−1vj = ∑4 i=0 ∑3−i j=0 (4− i− j) dijui+1vj∑4 i=0 ∑4−i j=1 jdiju ivj−1 = ∑4 i=0 ∑3−i j=0 (4− i− j) dijuivj+1 (4.12) if λ is eliminated and u = x1/x3 and v = x2/x3. After solving this for u and v the eigenvector and eigenvalues are x1 = u/ √ 1 + u2 + v2, x2 = v/ √ 1 + u2 + v2, x3 = ±1/ √ 1 + u2 + v2 and λ = d(~x). These four cases cover all possible combinations. In some implementations there have been suggestions that additional combinations of ~x in Case 3 is necessary (e.g., x1, x3 6= 0 and x2 = 0), but those are actually covered in Case 4 which include all combinations when x3 6= 0. 39 4.4. Z-EIGENVALUE DECOMPOSITION Table 4.1: The relation between the fifteen elements in the vector T and their relation to the fourth order symmetric tensor A used in the implementation of the Z-eigenvalues decomposition from [11] T(1) 1 ·A3333 T(9) 4 ·A1222 T(2) 4 ·A2333 T(10) 6 ·A1133 T(3) 6 ·A2233 T(11) 12 ·A1123 T(4) 4 ·A2223 T(12) 6 ·A1122 T(5) 1 ·A2222 T(13) 4 ·A1113 T(6) 4 ·A1333 T(14) 4 ·A1112 T(7) 12 ·A1233 T(15) 1 ·A1111 T(8) 12 ·A1223 4.4.2 Previous Software Implementations Two implementations were available, the first one can be found in [11]. This implemen- tations finds the set of eigenvalues and eigenvectors in three different cases. First, the vector T that contains the fifteen unique fourth order tensor elements is transformed into a symmetric fourth order tensor A according to Table 4.1. The factor corresponds to the multiplicity of the specific element. Then a fourth order tensor Ā is calculated element-wise according to Āijkl = ∆̄i: ·A:jkl (4.13) In the implementation ∆ = I where I is the identity matrix which means that Ā = A. The eigenvalues are then calculated in three different steps described below. The equations are presented in terms of the vector of the fifteen unique elements to facilitate comparison with the two other implementations. Step 1: If T14 = T13 = 0 then λ = T15 and ~x = (1 0 0). Step 2: For t solve the equation system −1 4T14 · t4 + (T15 − 1 2T12) · t3+ +3 4(T14 −T9) · t2 + (1 2T12 −T5) · t+ 1 4T9 = 0 1 4T6 · t3 + 1 4T11 · t2 + 1 4T8 · t+ 1 4T4 = 0 (4.14) The real solutions ti gives the eigenvectors ~x(i) and eigenvalues λ1 according to ~x(i) = (ti 0 0)/γi where γi = √ t2i + 1 and λ = T15 · (ti/γi)4 40 4.4. Z-EIGENVALUE DECOMPOSITION Step 3: For two variables, u and v equation system −1 4T13 · u4 − 1 4T11 · u3 · v + (T15)− 1 2T10) · u3 − 1 4T8 · u2 · v2+ +(3 4T14 − 1 2T7) · u2 · v + (3 4T13 − ·34T6) · u2 − 1 2T3 · u · v2− 1 4T4 · u · v3 + 1 2T12 · u · v2 + (1 2T11 − 3 4T2) · u · v+ (1 2T10 −T1) · u+ 1 4T9 · v3 + 1 4T8 · v2 + ·14T7 · v + 1 4T6 = 0 −1 4T13 · u3 · v + 1 4T14 · u3 − 1 4T11 · u2 · v2 + (1 2T12 − 1 2T10) · u2 · v+ 1 4T11 · u2 − 1 4T8 · u · v3 + (3 4T9 − 1 2T7) · u · v2 + (1 2T8− 3 4T6) · u · v + 1 4T7 · u+ 3 4T4 · v2 − 1 4T4 · v4 + (T5 − 1 2T3) · v3− 3 4T2 · v2 + (1 2T3 −T1) · v + 1 4T2 = 0 (4.15) is solved with Matlab’s solver ’solve’ and the eigenvectors and eigenvalues are calculated according to x(l) = (ul vl 1)/βl where βl = √ u2 l + v2 l + 1 and λl = ∑3 i=1 ∑3 j=1 ∑3 k=1 ∑3 m=1 Aijkm · x(i)(l)x(j)(l)x(k)(l)x(m)(l). This implementation works in most cases, but is sometimes slow. For T = [1.4751 0.7276 − 3.9916 − 0.0061 12.4934 − 2.7988 − 7.1486 17.7676 −3.6968 5.2683 0.0099 18.6609 − 4.3874 19.5485 5.9023] this takes almost half an hour to solve, but with the later implementation a few sec- onds. Additionally, if T consists of only zeros and ones this returns a set of com- plex solutions with a larger set of eigenvalues than the right amount, e.g. if T = [1 1 0 1 0 0 0 0 0 0 1 0 0 0 0] it returns a set of five complex eigenvectors (where the complex part is small but present) instead of three. The second implementation (by the first author, Alipoor, of [18]) available utilises the case-based method described in Section 4.4.1 using Matlab’s solve function for both Case 3 and Case 4. During simulations it occasionally crashed and it was observed that in Case 4 solve utilises back-substitution which can cause division by zero. 4.4.3 The Improved Implementation of Z-eigenvalues Decomposition This uses the cases described in Section 4.4.1 and is based on the implementation by Alipoor. The method to find the solutions to Case 1-3 was essentially kept, but sub- cases giving repeated solutions were isolated and removed. One additional step was implemented, confirming that the found solutions solved the original problem in Equation (4.10). To avoid the back substitution in Case 4 the computer algebra system MuPAD – which is included in Matlab’s Symbolic Math Toolbox – is called directly instead of via solve. Since all real solutions to (4.12) are requested the MuPAD-function numeric::polysysroots is used. This returns a numerical approximation of both real and complex solutions and to isolate the real solution a tolerance of 10−12 is set. All complex solutions with an 41 4.5. RESULTS OF FA-ESTIMATIONS imaginary part larger than the tolerance are removed, and if the imaginary part is smaller the solution is considered to be real. After removing the repeated solutions and verified that they all are solutions to the original problem (4.10), the eigenvalues and eigenvectors are ordered with largest eigenvalue first. 4.5 Results of FA-estimations Here a representative selection of the FA-estimation-results is presented. The complete set of results can be found in Appendix B. 4.5.1 Results for one Fibre FA-Estimations 5 10 15 20 25 30 35 40 45 50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 SNR F A − v a lu e DTI FA Low FA = 0.18 Medium FA = 0.51 High FA = 0.94 Figure 4.1: FA estimated with DTI, ground truth FA: FAlow = 0.18, FAmed = 0.51 and FAhigh = 0.94. In one fibre simulations the results of FA for DTI are presented in Figure 4.1, where when ground truth FA is high, the estimated FA-value was accurate. For the two lower FA-values, the estimates were more sensitive to noise, and the signal from the low ground truth FA-value was overestimated for SNR below 20. For the lowest SNR-level the low FA-value was estimated to be higher than the medium FA-value. The three FA-values generated from the bi-Gaussian model were for low and medium FA overestimated, but all three were close to ground truth for the high FA-value. The most accurate result was acquired from FAmin (see Figure 4.2). For FAmin the medium FA-value was consistently higher than the ground truth FA-value, and the high-FA-value was slightly underestimated. Both definitions for the fourth order model overestimated the medium and slightly underestimated the low. Qi’s definition was close to ground truth for the high FA-value. 42 4.5. RESULTS OF FA-ESTIMATIONS 5 10 15 20 25 30 35 40 45 50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 SNR F A − v a lu e Minimum of Bi−Gaussian Model Low FA = 0.18 Medium FA = 0.51 High FA = 0.94 Figure 4.2: Minimum of FA’s estimated with bi-Gaussian model, ground truth FA: FAlow = 0.18, FAmed = 0.51 and FAhigh = 0.94. In Figure 4.3 the results from FAQi are plotted. 5 10 15 20 25 30 35 40 45 50 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 SNR F A − v a lu e Fourth order model, Qi−definition FA Low FA = 0.18 Medium FA = 0.51 High FA = 0.94 Figure 4.3: FA-values estimated with Qi’s FA-definition for the fourth order model, ground truth FA: FAlow = 0.18, FAmed = 0.51 and FAhigh = 0.94. 43 4.5. RESULTS OF FA-ESTIMATIONS 4.5.2 Results for two Fibre FA-Estimations Consistently for all models, was that when one fibre had low FA-value, there was a very low dependence on crossing angle, e.g. see Figure B.21, B.23, and many others in Appendix B. DTI For DTI, when no ground truth tensor had low FA, the estimated FA decreased with increasing crossing angle. The higher the ground truth FA-values were, the more the estimated FA decreased with increasing crossing angle. This behaviour can be seen in Figure 4.4 where both synthetic fibres had high FA-value. 0 10 20 30 40 50 60 70 80 90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Crossing Angle F A − v a lu e DTI FA SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 Figure 4.4: Synthetic signal, two fibres with FA1 = FA2 = 0.94, DTI FA-estimation. Bi-Gaussian model The three bi-Gaussian FA-definitions were insensitive to the crossing angle. Excluding the high-high structure, when comparing the estimated FA-value to the mean of the ground truth FA-values both FAmax and FAmean were overestimated. In Figure 4.5 the FAmean of high-low and medium-medium-structure is plotted and the resulting FA values are all higher than the mean of the ground truth FA-values. Fourth Order Model Figure 4.6 show a typical result of the fourth order FA-estimations. When the crossing angle increases from 50◦ to 70◦ the estimated FA-value decreases rapidly. This behaviour was noted in all simulations when no fibre had low ground truth FA-value. In this rapid FA-estimation drop, the relation between SNR and FA changed. For crossing angles 44 4.5. RESULTS OF FA-ESTIMATIONS 0 10 20 30 40 50 60 70 80 90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Crossing Angle F A − v a lu e Bi−Gaussian Model, mean FA SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (a) Synthetic signal, two fibres with FA1 = 0.94 and FA2 = 0.18. Mean of FA’s esti- mated with the bi-Gaussian model. Ranged random initialisation. 0 10 20 30 40 50 60 70 80 90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Crossing Angle F A − v a lu e Bi−Gaussian Model, mean FA SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (b) Synthetic signal, two fibres with FA1 = 0.94 FA2 = 0.51. Mean of FA’s estimated with the bi-Gaussian model. Ranged ran- dom initialisation. Figure 4.5: Bi-Gaussian model, synthetic signal, two fibres with FA1 = 0.94 FA2 = 0.51. lower than 50◦ the FA-value increased with increasing SNR, but for angles higher than 70◦ the FA-value decreased with increasing SNR. If we compare DTI and the fourth order model, both methods’ FA-values decrease with increasing SNR, FAQi lowest value is approximately the same as the lowest DTI, and FAMA is lower than DTI. 0 10 20 30 40 50 60 70 80 90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Angle of separation F A − va lu e MA−definition HOT FA SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (a) Synthetic signal, two fibres with FA1 = FA2 = 0.94. FA’s estimated with MA’s FA-definition for the fourth order model. 0 10 20 30 40 50 60 70 80 90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Angle of separation F A − va lu e Qi−definition HOT FA SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (b) Synthetic signal, two fibres with FA1 = FA2 = 0.94. FA’s estimated with Qi’s FA-definition for the fourth order model. . Figure 4.6: Fourth order tensor model, synthetic signal, two fibres with FA1 = FA2 = 0.94. It was also observed that for perpendicular fibres the fourth order estimated FA- values was lower for the high-high fibre configuration than for high-medium. Compare Figure 4.6 and Figure 4.7. 45 4.5. RESULTS OF FA-ESTIMATIONS 0 10 20 30 40 50 60 70 80 90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Angle of separation F A − va lu e MA−definition HOT FA SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (a) Synthetic signal, two fibres with FA1 = 0.94 and FA2 = 0.51 FA’s estimated with MA’s FA-definition for the fourth order model. 0 10 20 30 40 50 60 70 80 90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Angle of separation F A − va lu e Qi−definition HOT FA SNR = 50 SNR = 25 SNR = 17 SNR = 13 SNR = 10 SNR = 8 SNR = 7 (b) Synthetic signal, two fibres with FA1 = 0.94 and FA2 = 0.51. FA’s estimated with Qi’s FA-definition for the fourth order model. Figure 4.7: Fourth order tensor model, synthetic signal, two fibres with FA1 = 0.94 and FA2 = 0.51. 46 4.5. RESULTS OF FA-ESTIMATIONS 4.5.3 Results of Real Data FA-Estimations In Figure 4.8 the conventional FA-map is shown. The non-Gaussian FA-maps were compared to this, given that no ground truth FA-values is available for real data. Figure 4.8: The conventional FA from DTI. For the bi-Gaussian model (see Figure 4.9), in many voxels an FA-value could not be estimated since neither of the estimated tensors was PSD. For these voxels the FA-value was set to zero, which turns out black in the image. This problem was not as significant for the synthetic signal, since many estimates and noise realisations for each structure was performed, and the non-PSD tensors could simply be ignored. However, if that artefact is put aside, the map is brighter than Figure 4.8 and has less contrast (ignoring black voxels). The fourth order model was fitted to the signal in two different ways. The first method described in Section 3.2.3, used for model evaluation of the fourth order tensor model in Chapter 3, yielded the FA-map shown in Figure 4.10. The contrast is much lower than in Figure 4.8. The second method, described in Section 4.3.3, yielded an FA-map with higher contrast (see Figure 4.11). The areas with high anisotropy (the bright areas) in this map, corresponds well with the DTI-FA-map. However the fourth order FA-map has fewer high anisotropic areas. Remember from the synthetic signal FA-maps for the fourth order model, that when the crossing angle increased to above 70◦ the FA-value rapidly decreased. 47 4.5. RESULTS OF FA-ESTIMATIONS Figure 4.9: The mean FA of the two FA from bi-Gaussian model. Figure 4.10: The Qi definition of FA from HOT, filtered with Gaussian window of 4 voxel size, fitted to the CT-ODF. 48 4.5. RESULTS OF FA-ESTIMATIONS Figure 4.11: The Qi definition of FA from HOT, filtered with Gaussian window of 4 voxel size, fitted to the ADC-profile. 49 4.6. CONCLUSIONS, DISCUSSION AND IMPROVEMENT OF FA-ESTIMATIONS 4.6 Conclusions, Discussion and Improvement of FA-estimations The first and easiest conclusion to draw from these experiments is that the FA-estima- tions for the three models give different results. FA computed from DTI, as previously mentioned, decreases with increasing crossing angle. It is well known that the conven- tional FA measure is lower in complex structures such as fibre crossings [20]. The three FA-definitions for the bi-Gaussian model are all invariant to the crossing angle. If this could be utilised the FA-value would not decrease in a complex area. However, since the bi-Gaussian FA-definitions consistently overestimates low anisotropy areas it is difficult to know if the FA-value reflects the ground truth or if it overestimates it. Consequently the bi-Gaussian FA-estimates are less informative than DTI FA-values. One possible improvement is to include a PSD-constraint in the model-fit used for the bi-Gaussian model. This would eliminate the negative eigenvalues and at the very least, remove the black voxels from the real data FA-maps. Furthermore, since the negative eigenvalues from the synthetic data simply are ignored and the mean values are calculated out of those estimated eigenvalues that are positive, the resulting FA-maps could be rather dif- ferent from those presented above. The effect of a PSD constraint would be interesting to evaluate and could be a continuation of the results presented here. The two FA-estimates derived from the fourth order tensor model were both ap- proximately invariant to the crossing angle below 50◦ and above 70◦ and exhibit rapid decrease in-between those two angles. The real data FA-maps for the fourth order model have fewer areas of bright areas, indicating less information than in the DTI FA-maps. One possible reason for this is that for DTI only six parameters were estimated, com- pared to 13 parameters for the fourth order model. The relation between number of measurements and parameters to estimate is for DTI larger than five, for the fourth order tensor model the relation is smaller than three. For future improvements it would be interesting to see if a larger number of gradient directions would improve the fourth order FA-maps and if they would provide more (but still accurate information) than the DTI FA-map. In discussions with dr. Ylva Lilja MD and PhD-student and Mohammad Alipoor, Lic. and PhD-student (see Section 4.1.2 for a presentation of their research) during the this research it was concluded that, clinically speaking, the desired behaviour is that an increased crossing angle naturally decreases the FA-value. The distinction between white matter-integrity and fibre-structure is key here. If the characteristic of interest is integrity, then the increase of crossing angle, probably should not decrease the integrity- value because whether or not the white matter fibres are crossing or not does not decrease their integrity. If instead the characteristic of interest is fibre structure, then increased crossing angle (and an increase of directions) is consistent with loss of structure and the measure should decrease. This was highlighted in [3]. For future work, the inclusion of a fraction of occupancy estimation for both the bi-Gaussian model and the fourth order model should be explored. For the bi-Gaussian model this requires multi-shell signal acquisition, but in a synthetic environment that 50 4.6. CONCLUSIONS, DISCUSSION AND IMPROVEMENT OF FA-ESTIMATIONS is just a matter of more simulation time. This could then lead to a more sophisticated FA-measure where the fraction between the fibres (Fraction of Occupancy – FoO) could be included as a weight to their respective eigenvalues. For the fourth order model, the current normalising factor ν (see Equation (4.4) and Equation (4.5)) give the same weight to all eigenvalues. If the FoO’s of these eigenvalues were known many small eigenvalues with low FoO would have less impact on the resulting FA value. The concept of a scalar measure that provides information about the underlying structure is good, but what exactly it should reflect is still an open question. Nevertheless this does not preclude the possibility of generalising the definition of FA for DTI to non- Gaussian models. However if such a measure has good accuracy and can differentiate between integrity and structure, white and grey matter, it would be superior to FA from DTI. 51 5 Summary and Conclusions I n this chapter the thesis and the work presented in it is summarised in Section 5.1. The major contributions and findings are presented in Section 5.2 and the conclusions in Section 5.3, and finally opportunities for further research in Section 5.5. 5.1 Thesis Summary Chapter 1 Introduces DWI (Diffu