Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in Diffusion MRI

Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in Diffusion MRI

Jian Cheng jian.cheng@nih.gov Peter J. Basser pjbasser@helix.nih.gov SQITS, NIBIB, NICHD, National Institutes of Health, United States
Abstract

In Diffusion Tensor Imaging (DTI) or High Angular Resolution Diffusion Imaging (HARDI), a tensor field or a spherical function field (e.g., an orientation distribution function field), can be estimated from measured diffusion weighted images. In this paper, inspired by the microscopic theoretical treatment of phases in liquid crystals, we introduce a novel mathematical framework, called Director Field Analysis (DFA), to study local geometric structural information of white matter based on the reconstructed tensor field or spherical function field: 1) We propose a set of mathematical tools to process general director data, which consists of dyadic tensors that have orientations but no direction. 2) We propose Orientational Order (OO) and Orientational Dispersion (OD) indices to describe the degree of alignment and dispersion of a spherical function in a single voxel or in a region, respectively; 3) We also show how to construct a local orthogonal coordinate frame in each voxel exhibiting anisotropic diffusion; 4) Finally, we define three indices to describe three types of orientational distortion (splay, bend, and twist) in a local spatial neighborhood, and a total distortion index to describe distortions of all three types. To our knowledge, this is the first work to quantitatively describe orientational distortion (splay, bend, and twist) in general spherical function fields from DTI or HARDI data. The proposed DFA and its related mathematical tools can be used to process not only diffusion MRI data but also general director field data, and the proposed scalar indices are useful for detecting local geometric changes of white matter for voxel-based or tract-based analysis in both DTI and HARDI acquisitions. The related codes and a tutorial for DFA will be released in DMRITool.

keywords:
Diffusion MRI, Diffusion Tensor, Orientation Distribution Function, Distortion, Dispersion, Director Field Analysis, Local Orthogonal Frame, Liquid Crystals, dyadic
journal: Medical Image Analysis

1 Introduction

Diffusion MRI is a powerful non-invasive imaging technique widely used to explore white matter in the human brain. Diffusion Tensor Imaging (DTI) (Basser et al., 1994) is used to reconstruct a tensor field from diffusion weighted images (DWIs). High Angular Resolution Diffusion Imaging (Tuch et al., 2002; Frank, 2002; Descoteaux et al., 2007; Tournier et al., 2007; Cheng et al., 2010, 2015, 2014; Özarslan et al., 2013), which makes no assumption of a 3D Gaussian distribution of the diffusion propagator, is used to reconstruct a general function field from DWIs, (e.g., an Orientation Distribution Function (ODF) or Ensemble Average Propagator (EAP) field). Both the ODF and the EAP fields with a given radius are spherical function fields. Exploring microstructural information from the reconstructed tensor field or spherical function field is of interest in many biological and clinical application areas, which makes diffusion MRI a powerful means to study white matter. For example, in an voxel exhibiting anisotropic diffusion, local peaks of the reconstructed spherical function or the principal eigenvector of the reconstructed 2nd-order diffusion tensor normally prescribe the fiber directions in that voxel.

Some scalar indices have been proposed to be estimated voxel-wise from tensors/ODFs/EAPs. For DTI, well-established tensor scalar indices, including the mean diffusivity and Fractional Anisotropy (FA), are widely used as biologically meaningful descriptors (Pierpaoli and Basser, 1996). Kindlmann et al. (2007) proposed two sets of scalar indices (three scalar indices per set), which are orthogonal in terms of tensor changes and the Euclidean inner product. For High Angular Resolution Diffusion Imaging (HARDI), the generalized FA (Tuch, 2004), Orientation Dispersion index (OD) (Zhang et al., 2012), return-to-origin probability (Helmer et al., 2003; Wu and Alexander, 2007), and mean-squared displacement (Basser, 2002; Wu and Alexander, 2007) were all proposed for ODFs and EAPs. These indices indicate some information inside a voxel, but cannot describe local geometric or topological information, including fiber crossing, fanning, bending, and twisting, in a local spatial neighborhood.

Some previous works have extracted local geometric information by considering the local spatial change of tensor fields or ODF fields. Pajevic et al. (2002) demonstrated that the norm of the spatial gradient of the tensor’s isotropic and anisotropic parts can detect boundaries between white matter, CSF, and gray matter. Kindlmann et al. (2007) proposed tangents of scalar invariants and rotation tangents, which are 2nd-order tensors, and also proposed projecting the 3rd-order spatial gradient tensor onto these 2nd-order tangents to obtain the spatial direction with the largest change of scalar indices or rotation of tensors. Based on the rotation tangents of tensors, Savadjiev et al. (2010) proposed fiber curving and fiber dispersion indices. Tax et al. (2016) proposed a sheet probability index to quantify the local sheet structure by using spatial changes of ODF peaks. Duits and Franken (2009, 2011); Portegies et al. (2015) proposed spatial and spherical smoothing to enhance an ODF field in a PDE framework, preserving crossing structures. Reisert and Kiselev (2011); Cheng et al. (2013); Michailovich et al. (2011) considered spatial coherence in ODF estimation.

The terms “splay”, “bend” and “twist” have been used to qualitatively describe complex local white matter structural configuration in literature of diffusion MRI for about 20 years (Basser, 1997; Pajevic et al., 2002; Johansen-Berg and Behrens, 2009). However, to our knowledge, there is no existing work that quantitatively describes the degree of local orientational change of white matter, including splay, bend, and twist, from general ODF fields in dMRI, although the fiber curving and dispersion indices in Savadjiev et al. (2010) can be seen to quantify “splay” and “bend” for a tensor field in DTI.

Basser (1997) discussed the initial idea to study the torsion and curvature of a fiber tract by using the Frenet frame 111https://en.wikipedia.org/wiki/Frenet-Serret_formulas along the tract. Torsion and curvature from the Frenet frame were later used in diffusion data analysis in Batchelor et al. (2006). Savadjiev et al. (2007) used the Frenet frame as a prior in the relaxation labeling algorithm to regularize the data and estimate ODFs in voxels. These works on the Frenet frame studied geometric information along a single tract. However, tractography is known to be sensitive to a large number of parameters, and any flaws in the reconstructed tracts due to noise or parameter selection will inevitably be reflected in the geometric information that is extracted subsequently. Piuze et al. (2015) proposed moving frames determined by the geometry of cardiac data, and calculated Maurer-Cartan connections. However this method is not applicable to general diffusion MRI data, and does not consider the sign ambiguity in the frame.

There exist some connections between diffusion MRI data analysis and liquid crystals. Orientational order parameter is well-established to describe the degree of alignment in liquid crystals Andrienko (2006). Lasič et al. (2014); Szczepankiewicz et al. (2015) calculated the order parameter map by estimating variance of microscopic diffusion parameters from the contrast between diffusion signals measured by directional and isotropic diffusion encoding. However, it cannot be used for general DTI and HARDI data. Topgaard (2016) used a diffusion tensor method to estimate the director orientations of a lyotropic liquid crystal as a spatially resolved field of Saupe order tensors.

In this paper, inspired by orientation and distortion analyses applied to liquid crystals, we propose a unified framework, called Director Field Analysis (DFA), to study the local geometric information of white matter from the reconstructed spherical function field. DFA works both for tensor fields obtained from DTI and for spherical function fields from HARDI. At the voxel level, 1) the Orientational Order index (OO) and Orientational Dispersion index (OD) are defined for the spherical function in a voxel with a given axis (e.g., the ODF with its principal direction); and 2) the principal direction is extracted from the spherical function in such a voxel exhibiting anisotropic diffusion. At a local neighborhood level, 1) an orthogonal coordinate frame is defined for each voxel with anisotropic diffusion, where the first axis is the extracted principal direction; 2) OO is defined for spherical functions in a local neighborhood with the given principal direction; and 3) three distortion indices (splay, bend, twist) and a total distortion index are defined based on the spatial directional derivatives of the principal direction. An overview of the DFA pipeline for a spherical function field is shown in Fig. 1.

Figure 1: Director Field Analysis (DFA) pipeline for an ODF field obtained from DTI or HARDI. DFA provides total six scalar indices calculated from a spherical function field at the voxel level and at the local neighborhood level.

This paper is organized as follows. Section 2 provides a unified overview of existing works on tensor field analysis for exploring local geometric information (Pajevic et al., 2002; Kindlmann et al., 2007; Savadjiev et al., 2010), which is also a motivation for the proposed DFA. Section 3 proposes the DFA framework that works for both diffusion tensor fields and ODF fields. Section 4 demonstrates some results of synthetic and real data experiments by using DFA. Section 5 discusses some issues on implementing DFA.

2 Tensor Field Analysis

This section provides an overview of existing works exploring the local geometric features of a 2nd-order diffusion tensor field by using the spatial gradient of tensors (Pajevic et al., 2002; Kindlmann et al., 2007; Savadjiev et al., 2010) in a unified framework. It also proposes a new 4th-order structure tensor applied to 2nd-order tensor data that generalizes the conventional 2nd-order structure tensor applied to scalar fields.

The 2nd-order diffusion tensor field. In a diffusion tensor field denoted as , there is a diffusion tensor at the voxel x, where , and is the set of symmetric positive definite matrices. The diffusion tensor is symmetric with six unique (i.e., independent) elements.

The 3rd-order spatial gradient of the diffusion tensor. For a tensor field denoted as with elements at the voxel x, its spatial gradient at voxel x, denoted as , is a 3rd-order tensor with elements , where . Since the diffusion tensor is symmetric with six unique elements, the 3rd-order spatial gradient has 18 unique elements.

Mapping the 3rd-order spatial gradient to a vector. Let be a designed 2nd-order weighting tensor in tensor space, then the tensor inner product

(1)

produces a vector in the image space that is the spatial gradient of the scalar field at the voxel x. Note that the inner product in Eq. (1) is performed at voxel x, and the x dependency is omitted in the notation if there is no ambiguity. There are several ways to design a physically meaningful weighting tensor, . could be a constant independent of spatial position, x, or a function of x. 1) If , then is the mean diffusivity field, and Eq. (1) is its spatial gradient. 2) If is a scalar function that maps to a scalar value, then is the gradient of the scalar function, which is also a 2nd-order tensor with elements . If we set , then the vector in Eq. (1) is the spatial gradient of the scalar field , because of by the chain rule. If is the mean diffusivity function, then . We can also use other scalar invariants of tensors, e.g., FA. 3) If we choose as the rotation tangent  (Kindlmann et al., 2007) defined as the change of tensor value due to infinitesimal rotations around the -th eigenvector, then Eq. (1) denotes the direction in which the tensor orientation around the -th eigenvector varies the fastest.

Mapping the 3rd-order spatial gradient to a scalar value. Let be a 2nd-order weighting tensor in tensor space, and let v be a vector, then the tensor inner product

(2)

produces a scalar value that is the directional derivative of the scalar field at voxel x along the vector v, and is also the weighted mean of the directional derivative of along the vector v. 1) If we set , then Eq. (2) is the squared norm of the tensor gradient, which is useful for detecting boundaries of a tensor field (Pajevic et al., 2002). 2) By choosing v as the three eigenvectors of , and as three rotation tangents around three eigenvectors, respectively, we have total 9 scalar values to distinguish 9 configurations of tensor fields (Savadjiev et al., 2010). 3) The above 9 scalar indices can be combined to devise the fiber curving and fiber dispersion indices (Savadjiev et al., 2010).

The 4th-order structure tensor. We propose a new 4th-order structure tensor with elements , which is analogous, but generalizes the structure tensor of a scalar field 222https://en.wikipedia.org/wiki/Structure_tensor. The above 4th-order structure tensor is minor symmetric (Moakher, 2009), i.e., , . Thus, there are 36 unique elements out of a possible total of elements, and there is one-to-one mapping between this 4th-order tensor and a 2nd-order tensor (i.e., a matrix). However, since the 4th-order tensor is minor symmetric, the corresponding matrix is not symmetric in general. Thus, eigenvalues and the 2nd-order left and right eigenvectors can be calculated based on eigen-decomposition of the non-symmetric matrix. We may define some scalar invariants from these six eigenvalues of the 4th-order structure tensor, which can be used as features in this high dimensional space. We can also contract the 4th-order structure tensor to a scalar value by using the tensor inner product , which is the weighed mean of the squared spatial directional derivative along vector v. When setting v as three eigenvectors and corresponding weighting as rotation tangents divided by the spatial gradient, then the tensor inner product produces nine scalar indices that are the squares of the corresponding nine indices in Savadjiev et al. (2010). Thus, the curving and dispersion indices in Savadjiev et al. (2010) can also be obtained by choosing the corresponding .

3 Method: Director Field Analysis

Section 2 provides a unified framework to explore geometric structure information (e.g., boundary, curving, dispersion, etc.) from a tensor field, by considering a different weighting matrix on the spatial gradients. However, it is challenging to generalize this framework to ODFs in HARDI, where ODFs are normally general spherical functions with antipodal symmetry.

In this section, we propose a novel mathematical framework, called Director Field Analysis (DFA). Section 3.1 defines director related concepts to deal with vectors with sign ambiguity. Section 3.2 provides a set of mathematical tools for DFA. Section 3.3 proposes OO and OD for tensors and ODFs in voxels and in a spatial neighborhood, and gives closed-form results in some specific cases. Section 3.4 extracts the principal direction and its related local orthogonal frame in voxels exhibiting anisotropic diffusion. Section 3.5 defines four orientational distortion indices and demonstrates the implementation of the calculation by using the local orthogonal frame in Section 3.4. Section 3.1 and 3.2 are the theory part of DFA. Section 3.3, 3.4, and 3.5 are the application part of DFA in diffusion MRI. Fig. 1 demonstrates DFA to a spherical function field obtained from DTI or HARDI.

3.1 Director and Director Field

We define a director as a unit norm vector v that is equivalent to . The director term is borrowed from studies of liquid crystals 333https://en.wikipedia.org/wiki/Liquid_crystal. We also define a director with weight, or weighted directors, as a vector associated with a weight , which is equivalent to , where , . If , then a weighted director can be represented as . See Fig. 2. A director v can be uniquely represented as a dyadic tensor, , which avoids the sign ambiguity, and a director with weight can be uniquely represented as a dyadic tensor, . We define a director field as , where there are some weighted directors in each voxel x.

Directors occur very often in diffusion MRI studies. Eigenvectors of diffusion tensors, local maxima of ODFs, and local fiber directions are all directors. Based on eigen-decomposition, a tensor is the sum of three dyadic tensors that represent three weighted directors. A spherical function which satisfies antipodal symmetry, i.e., , can be seen as infinite weighted directors . Thus, a spherical function field, , is a director field by definition.

An ODF in a voxel exhibiting anisotropic diffusion is anisotropic, and the orientations where the ODF takes its local peak (i.e., local maximal values) are normally considered to be local fiber directions in that voxel. A normal peak detection algorithm for ODFs performs a grid search in a spherical mesh, and then refines the solution by using a gradient ascent on the continuous sphere (Tournier et al., 2004). Note that peak detection is only performed for voxels exhibiting anisotropic diffusion (e.g., where ODFs have Generalized FA (GFA) (Tuch, 2004) values larger than ). Moreover, in order to avoid including small peaks produced by noise, only peaks whose values are larger than a threshold percentage (e.g., ) of the largest ODF value are counted. After peak detection, for each voxel x, we obtain a discrete spherical function from the continuous spherical function , where are local peaks. This discrete spherical function field is also a director field, or called a peak field, which emphasizes local peaks and suppresses weights for other directors. A peak field can also be extracted from a tensor field. In each voxel for a tensor field, there is 0 or 1 peak, and the principal eigenvector of the tensor is considered as a peak, if the tensor has a large FA value (e.g., larger than ).

3.2 Mathematical Tools for Directors

We provide a set of mathematical tools for analyzing directors and director fields. These tools are useful not only for this paper, but also for other applications which deal with continuous or discrete director data.

(a) (b) (c)
Figure 2: The mean, main, and difference of directors. Directors , , are visualized as vectors and , and the length of is the positive weight . (a) the mean (in blue) and main (in red) directors of two directors (in black), where , and . See Proposition 2. (b) the mean (in blue) and main directors of three directors (in black) which are orthogonal to each other, where , . The mean director is not unique, because with an arbitrary sign assignment can be the mean director. The main director is . (c) the difference of two directors. The two blue vectors denote the director representation of the difference, i.e., . The two red arcs denote the rotation matrix representation of the difference, i.e., , where is a scaled rotation matrix such that . The rotation matrix representation has no sign ambiguity, while the director representation has the sign ambiguity.

3.2.1 Mean Director of a Set of Directors

For a given weighted directors , if we convert a director to a vector by assigning a sign, we have total possible sign assignments. Thus, we have possible Euclidean mean vectors for the vectors. We define a mean weighted director of a set of weighted directors as the Euclidean mean vector with the sign assignment that takes the maximal norm among the mean vectors.

Definition 1.

Mean director of weighted directors. A mean weighted director of a set of weighted directors is defined as , where the signs , and is the mean director operator.

It is obvious that and have the same mean director. Thus, without loss of generality, we assume non-negative weights for calculating the mean director. If the angle between any two vectors and is no more than , then the sign assignment for the mean director can be proved to be , . See Proposition 1 whose proof is based on the proof of the mean director of two directors, which is trivial. The mean director may be not unique when some directors are orthogonal.

Proposition 1.

Mean director of weighted directors in a cone. For a set of weighted directors with non-negative weights, if all directors are in a cone, i.e., , , then the mean weighted director is .

3.2.2 Main Director of a Set of Directors

We define the main director of a set of weighted directors as the main axis in Principal Component Analysis (PCA) by using eigen-decomposition. This concept is from the average director of molecule orientations in liquid crystals.

Definition 2.

Main director of weighted directors. A main weighted director of a set of weighted directors is defined as , where is the eigenvalue of the tensor which has the largest absolute value among all eigenvalues, and is its corresponding eigenvector, and is the main director operator.

Note 1) The dyadic tensor of the largest eigenvalue and eigenvector is the best rank-1 approximation of in terms of the L2 norm. 2) The main director may not be unique, considering there may be more than one eigenvalues which are equal, and are all the largest eigenvalue. 3) Unlike the mean director which is independent of the signs of the weights, the main director is dependent on the weight signs. 4) Although, in general, the mean and the main directors are not the same, in some cases, they may give the same direction. See Proposition 2 and Fig. 2 (a) for the two director case. See Fig. 2 (a) and (b) for an illustration of the mean and main directors.

Proposition 2.

Two weighted directors with the same weight. For two weighted directors with the same weight, denoted as and , the main director is if , and is if . The mean director is , if , and is , if .

The mean director and main director describe different meaningful information about directors. Take a diffusion tensor , , as an example. There are three weighted directors , or represented as because of non-negative weights. The mean director is with any sign assignment of , while the main director is . See Fig. 2 (b). Additionally, small changes in , , , and may change the mean director, but not the main director, if and are still the largest eigenvalue and eigenvector. This example clearly shows that the mean director concept is a generalization of the mean vector concept in vector space, while the main director emphasizes the main axis in PCA. Please note that in a general case, the change of any director may cause the change of the main director (i.e., the largest eigenvector and eigenvalue of the tensor ) and also the mean director.

3.2.3 Two Representations of the Difference between Two Directors

We aim to generalize the tensor field analysis in Section 2 to director fields (i.e., ODF fields), and explore geometric structure information by using spatial derivatives which rely on the concept of difference between two directors.

We propose two ways, i.e., the director representation denoted as and the rotation matrix representation denoted as , to represent the difference between two weighted directors with non-negative weights, and . These two directors can be converted to the vectors and by assigning the sign (or ), and such that . Thus, there are two different cases because of the sign ambiguity. We can represent the difference as a director, i.e., . We can also represent the difference as a scaled rotation matrix which rotates to , i.e., , . The rotation matrix can be calculated from the rotation axis (i.e., the cross product of and ) and the rotation angle (i.e., the angle between and 444https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle, and the scale can be calculated from the weights and . Note that this rotation matrix is the same for the above two cases, without sign ambiguity. See Fig. 2 (c) as an illustration. The director representation of the difference has the sign ambiguity, but it gives a vector without a sign which can be projected onto axes. The rotation matrix representation is unique without sign ambiguity, but cannot be projected onto axes.

3.2.4 Spatial Gradient and Directional Derivative of a Director Field

Considering a director field where there is only one director with non-negative weight (simplified notation for ) at each position , a directional derivative along at x is defined as

(3)

Thus, there are also director and rotation matrix representations of the directional derivative because of the two representations of Diff, i.e., and .

For the director field where directors are only obtained in a integer lattice, the central difference can be used to approximate the spatial gradient , where

(4)

, , are the unit norm vectors along spatial axes, and Mean is the mean director operator in Definition 1.

We normally use the rotation matrix representation for the spatial gradient , considering this representation is unique. Let be the rotation matrices for the spatial gradient at x along axes , i.e., is used in Eq. (4). Then, analogously to the spatial gradient of a vector field, the director at position with small can be approximated as the sum of rotated directors, i.e.,

(5)

Note that every director in the above sum is a rotated in a small local rotation, thus we assume all directors are in a cone to obtain a simple sum of vector representation. See Proposition 1. If Eq. (4) is used to approximate , then the angle between and is no more than , thus all three rotated directors are indeed in a cone.

3.3 Orientational Order and Dispersion

Before working on a field of ODFs, an ODF in a voxel can provide some geometric information at the voxel level, including GFA (Tuch, 2004), and orientation dispersion (Zhang et al., 2012).

3.3.1 Orientational Order Transform and Orientational Tensor

The NODDI model is increasingly used to study neurite orientation dispersion (Zhang et al., 2012). NODDI uses the Watson distribution in Eq. (6) to model the ODF with a single orientation, where is the confluent hypergeometric function, is a given axis

(6)

Note that the original formula of the Watson distribution in Zhang et al. (2012) has no unit integral in , because it missed . An orientation dispersion index (OD) was defined as Eq. (7), where we denote it as because it only applies to the Watson distribution.

(7)

Note that in order to obtain good contrast in the dispersion index map, in the NODDI toolbox provided by the authors, a scaled ( in the codes) is used in Eq. (7) to calculate , instead of the estimated from the NODDI model. can not be used for ODFs that have general shapes, have more than one peak, or are not antipodally symmetric. Some other works also proposed dispersion indices based on different models of ODFs, e.g., Bingham distributions Tariq et al. (2016) and mv- distributions in DIAMOND Scherrer et al. (2015). These dispersion indices cannot work for general ODFs. Inspired by liquid crystals, we would like to define the degree of dispersion for general ODFs, independent of microscopic diffusion signal models.

For a general spherical function , , we define the orientational tensor as

(8)

which is related to the Q-tensor in liquid crystal modeling (Andrienko, 2006) 3. is a symmetric matrix dependent on . If is non-negative, then is positive semidefinite. We propose the orientational order index (OO) from the theory of liquid crystals (Andrienko, 2006) to describe the orientation or dispersion of a general spherical function along a given axis n:

(9)

where is the 2nd-order Legendre polynomial. By definition, Eq. (9) is an integral transform in which converts the spherical function to another spherical function , and the kernel is , similar to the Funk-Radon transform used in Q-Ball imaging (Tuch, 2004), where the kernel is . We call Eq. (9) the Orientational Order Transform (OOT), i.e., . Note that we have

(10)

By definition, is antipodally symmetric and has a global maximum and a global minimum on the unit sphere, which correspond to the largest and smallest eigenvectors of , respectively. Based on Definition 2, the main director of infinite weighted directors is the maximum point of .

Figure 3: A cross-section view of a spherical function along axis n. The projection of onto n is .

Although is a spherical function, it is a scalar index when n is chosen as a physically meaningful axis, e.g., takes its maximal value at n. Let be the angle between and axis n, then . Thus, if is a Probability Density Function (PDF) on the unit sphere, then is , where signifies the expectation operation. As shown in Fig. 3, is the expectation of the squared projected length of onto the axis n. If is more concentrated along n, then is larger, so is . Based on the definition, when is a PDF, then we have . If , i.e., the delta function along a given axis, then . If , such that , i.e., is non-zero only in the plane orthogonal to , then . If is the isotropic PDF, i.e., , then . In practice, if we choose the axis n such that takes its local or global maximal value, then is normally non-negative.

We define the orientational dispersion along axis n as

(11)

Then .

Note that the proposed OO is different from the order parameter in Lasič et al. (2014); Szczepankiewicz et al. (2015) which was also inspired by liquid crystals (Andrienko, 2006). In Lasič et al. (2014); Szczepankiewicz et al. (2015), the order parameter is calculated by estimating the variance of microscopic diffusion parameters from the contrast between signals measured by directional and isotropic diffusion encoding. However, it cannot be used for general DTI and HARDI data. The proposed OO in this paper is defined for general spherical functions (i.e., ODFs) along a given axis, independent of microscopic diffusion models and reconstruction of the ODFs.

3.3.2 Axisymmetric Spherical Functions

When is axisymmetric, and its axis is given by , i.e., , where is the corresponding scalar function defined in , then OOT has a closed form:

(12)

where is the angle between n and the axis , and is the 2nd-order Legendre coefficient of . Note that if , when , , then is the global maximum of . When , , then is the global minimum of . In the following development, without any ambiguity, we will use OO to denote , and OD to denote , for axisymmetric spherical functions.

3.3.3 Watson Distributions

The Watson distribution defined in Eq. (6) is axisymmetric with the axis . Thus, based on the above analysis of axisymmetric spherical functions, we have , and

(13)

where is the imaginary error function. Then . The left part of Fig. 4 shows the above OD and as functions of , where the axis n is set as the Watson distribution’s axis . Both dispersion indices decrease as increases. Based on the derivation of , is more sensitive to changes of when is small (), while it is less sensitive when is large (). Compared with , the change of OD is smoother for the change of over the entire range of .

3.3.4 Tensors

For the tensor model in DTI, denoted as , OOT is defined for its ODF, i.e.,

(14)

which is a PDF on the unit sphere. When the three eigenvalues of satisfy , is an axisymmetric function with the axis that is the principal eigenvector of . OOT has a closed-form expression in Eq. (12), and

(15)

The right panel of Fig. 4 shows OO and FA as functions of , where we set . Both OO and FA increases as increases. Thus, OO can be seen as a type of anisotropy index for tensors. For general tensors with , no such closed form solution like Eq. (15) and Eq. (12) exists, but we can calculate OO using the spherical harmonic representation of the ODF.

Figure 4: Left: dispersion indices (OD and ) of a Watson distribution as functions of . Right: OO and FA of prolate tensors () as functions of .

3.3.5 Spherical Harmonic Representation

For a general spherical function , OO and OD can be analytically calculated from the spherical harmonic coefficients of the rotated function. Considering is a real function on the unit sphere, it can always be linearly represented by the real Spherical Harmonic (SH) basis , i.e.,

(16)
(17)

where is the complex SH basis, is the associated Legendre polynomial. For any rotation matrix, the SH coefficients of the rotated function can be calculated with very high accuracy based on the Wigner D-matrix 555https://en.wikipedia.org/wiki/Spherical_harmonics, or based on fitting rotated function samples (Lessig et al., 2012). Let be the rotation matrix which rotates the axis n to z-axis, and be the real SH coefficients of the rotated function , considering the orthogonality of the real SH basis and , we have

(18)

Note that is only determined by the rotated SH coefficient that is only related to and the axis n, based on the rotation property of the SH basis. Thus, is only related to the SH coefficients of with , and also the axis n.

3.3.6 Relationship Between OO, OD, and GFA

For an ODF in an SH representation in Eq. (16), its GFA (Tuch, 2004) is

(19)

Note that the rotation of a spherical function does not change the shape of the function and the norm of SH coefficients, thus we have . Based on Eq. (18), we have

(20)

Combining Eq. (19) and Eq. (20), we have

(21)

The above inequality gives an upper bound of as a function of GFA which is independent of n. Note that the above upper bound is tight, and the equality holds when , for , and for after rotation. If the ODF has unit integral, i.e.,  666Note that ODFs estimated by some methods (e.g., constrained spherical deconvolution (Tournier et al., 2007), Q-Ball Imaging (Tuch, 2004; Descoteaux et al., 2007)), do not have the unit integral, if there is no normalization after estimation., then , and we have

(22)

Thus, for an ODF with low GFA, OO is also low, and OD is high, no matter how we choose the axis n. Note that Eq. (22) does not imply that an ODF with high GFA tends to have high OO, because it is an upper bound of , not a lower bound.

3.3.7 Mixture Model

OOT in Eq. (9) is a linear transform. Thus, if is the PDF of a mixture of models, where is the PDF for the -th model, and is the weight, then is also a mixture of OO functions. Fig. 5 illustrates OO for a two-tensor model with a crossing angle , where two tensors share the same eigenvalues , the weights are and , and one tensor component is along the -axis and the other one rotates from the -axis to the -axis. Based on Eq. (12) and Eq. (15), OO for the mixture model can be analytically calculated.

3.3.8 OO and OD for a General ODF Along the Principal Peak

In the above context, we focus on and as spherical functions. A physically meaningful axis is needed to obtain scalar indices of OO and OD from and . For an axisymmetric function , its axis can be used as described above. For a general function (e.g., an ODF), we can set the axis as the local maxima of (e.g., detected peaks of the ODF), because the peaks of ODFs are considered as local fiber directions in dMRI. A general ODF may have more than one peak. The principal peak of the anisotropic ODF , where the ODF takes its global maximum , i.e., , , is used to calculate OO and OD for the ODF. Note that peaks are detected from ODFs with all orders of SH coefficients, not only SH coefficients with . Thus, the scalar indices of OO and OD are actually dependent on SH coefficients of ODFs with all orders. See Algorithm 1 for the pipeline to calculate OO and OD maps from a given ODF field with SH representation, where peaks are detected for voxels whose GFA values are larger than a given threshold (e.g., ). It is also possible to calculate OO and OD for all voxels by setting the GFA threshold as . As shown in Section 3.3.6, for the voxels with , we have and .

Input: ODF field in SH representation.
Output: OO map, OD map.
Peak detection using gradient ascent for ODFs in voxels with the anisotropy higher than a given threshold (e.g., ). See Section 3.1;
for each voxel x with detected peaks do
       1) Find the principal peak with the largest ODF value, i.e., , ;
       2) Calculate rotation matrix , which rotates to the z-axis ;
       3) Calculate the rotated SH coefficient from under the rotation ;
       4) as shown in Eq. (18), and ;
      
end
Algorithm 1 Calculation of OO and OD for ODFs with SH representation along principal peaks:
Figure 5: OO for the mixture tensor model. Left: OO as a function of the angle between two tensor components. Right: ODF glyphs of the two-tensor model for different crossing angles, where the yellow tube shows the -axis which is used to calculate OO.

3.3.9 OO, OD and the Orientational Tensor in a Spatial Region

The above OO, OD, and the orientational tensor are defined for a single voxel. They can also be defined for voxels in a spatial region of voxels. A linear weighting generalization of OO can be defined as

(23)

The orientational tensor in a spatial region is

(24)

Because of the linearity of the integration, Eq. (23) is actually OOT in Eq. (9) performed on the region smoothed spherical function , and Eq. (24) is the orientational tensor for the region smoothed function. The largest eigenvector of in Eq. (24) indicates the main orientation of all ODFs in the region .

3.4 Local Orthogonal Frame

(a) (b) (c)
Figure 6: Sketch to determine local orthogonal frames from an ODF field, where an ODF may have 0, 1, or more than 1 peaks. (a) an ODF field with peaks, where yellow tubes denote peaks. (b) the orthogonal plane for the principal peak, where red tubes denote principal peaks. (c) local orthogonal frames, where three tubes in red, green, and blue colors denote three directors in local orthogonal frames.

As described in Section 3.1, after peak detection on a spherical function field or a tensor field, the obtained peak field is also a director field. We propose extracting a local orthogonal frame in each voxel exhibiting anisotropic diffusion from the detected peak field. The orthogonal frame has three orthogonal orientations. Denote the peaks at voxel x as . The first orientation is the principal peak where the ODF takes its global maximum , i.e., , . We call it the principal director of the voxel x. The other two orientations are in the orthogonal plane of the principal direction. Considering is normally antipodally symmetric in diffusion MRI, all these orientations are equivalent to their antipodal ones. Thus, we project all peaks in a spatial local neighborhood onto the orthogonal plane, and define a weighted sum of dyadic tensors in voxel x:

(25)
(26)

where is a local neighborhood of voxel x, is the spatial weight which is normally set to be proportional to , which is normally set as 1 voxel controls spatial weight concentration, is the projected orientation onto the orthogonal plane of . The above matrix is actually the orientational tensor of all projected peaks in region based on Eq. (24), where the continuous integral is replaced by a discrete summation over all projected peaks in region . Note that although we can define using continuous ODF with projected directors in a continuous integration like Eq. (24), we choose a discrete summation over peaks, which actually focuses only on peaks and sets zero weights for orientations that are not peaks in the continuous integration. in Eq. (25) has at most two non-zero eigenvalues, because it is defined by using in the orthogonal plane. The eigenvector for the largest absolute eigenvalue of is set as the second orientation of the orthogonal frame, which is the main director of directors , and indicates the main orientation of the local spatial change of in the orientational plane. Note that we define using the isotropic spatial weight to capture the general spatial change of the principal director in the orthogonal plane. If one has a good motivation and specific spatial prior knowledge (e.g., to capture local change only in a specific region like hippocampus), an anisotropic spatial weight with consideration of spatial prior knowledge may be useful. The third orientation in the orthogonal frame is set as the cross product of the first two orientations. These three orientations are three orthogonal directors due to sign ambiguity. Please see the sketch map in Fig. 6 to determine local orthogonal frames from a given ODF field. If these two eigenvalues of are equal or their difference is very small, then we set the second and third orientations in the orthogonal frame to zero, which means any two orthogonal vectors in the orthogonal plane can be the second and third axes in the orthogonal frame.

3.5 Local Distortion Indices: Splay, Bend, and Twist

splay bend twist
Figure 7: Demonstration of three types of distortions, i.e., splay, bend, and twist.

Three types of orientational distributions in liquid crystals. Based on the liquid crystal analogy, there are three fundamental types of distortions 3 for the director field as demonstrated in Fig. 7. 1) splay: bending occurs perpendicular to the director; 2) bend: bending is parallel to the director and molecular axis; 3) twist: neighboring directors are rotated with respect to one another, rather than aligned. These three fundamental distortions can be used to describe a myriad of complex geometric patterns that liquid crystals can assume. We would like to quantify these fundamental distortion patterns in dMRI by exploring the local spatial changes of principal directors.

Spatial derivatives of the local orthogonal frame. With the local orthogonal frame at each voxel x obtained above, we can define the spatial directional derivatives of along a direction v as

(27)

is the director representation of the difference of two directors as described in Section 3.2.3. Note that we use the director representation for the spatial derivative, instead of a rotation matrix representation, because we would like to project the director onto different axes. See Section 3.2.3 and 3.2.4.

Spatial derivatives of vectors, and Maurer-Cartan connection forms in the moving frame method. If we assume are all well-aligned unit vectors (i.e., no sign ambiguity), then we have

(28)

with elements , where is the -th element of , and is the spatial gradient matrix of . Similarly with Eq. (2), we can extract some features by devising v and a weighting vector w:

(29)

Eq. (29) can be seen as a generalization of Eq. (2) in tensor field analysis. When we set and , Eq. (29) is the projection of the directional derivatives onto , denoted as :

(30)

is the Maurer-Cartan connection form in the moving frame method 777https://en.wikipedia.org/wiki/Maurer-Cartan_form.