Asymmetric Orientation Distribution Functions (AODFs) revealing intravoxel geometry in diffusion MRI

Magnetic Resonance Imaging 49 (2018) 145–158

Contents lists available at ScienceDirect

Magnetic Resonance Imaging journal homepage: www.elsevier.com/locate/mri

Asymmetric Orientation Distribution Functions (AODFs) revealing intravoxel geometry in diffusion MRI

T

Suheyla Cetin Karayumaka,b, Evren Özarslanc,d, Gozde Unale,* a b c d e

Psychiatry Neuroimaging Laboratory, Brigham and Women's Hospital and Harvard Medical School, USA Faculty of Engineering and Natural Sciences, Sabanci University, Turkey Department of Biomedical Engineering, Linköping University, Linköping, Sweden Center for Medical Image Science and Visualization, Linköping University, Linköping, Sweden Department of Computer Engineering, Istanbul Technical University, Turkey

A R T I C LE I N FO

A B S T R A C T

Keywords: ODF regularization Asymmetric fiber orientations Fiber asymmetry measure Asymmetric ODF (AODF) Directional spatial filtering Steerable filtering Diffusion MRI HARDI

Characterization of anisotropy via diffusion MRI reveals fiber crossings in a substantial portion of voxels within the white-matter (WM) regions of the human brain. A considerable number of such voxels could exhibit asymmetric features such as bends and junctions. However, widely employed reconstruction methods yield symmetric Orientation Distribution Functions (ODFs) even when the underlying geometry is asymmetric. In this paper, we employ inter-voxel directional filtering approaches through a cone model to reveal more information regarding the cytoarchitectural organization within the voxel. The cone model facilitates a sharpening of the ODFs in some directions while suppressing peaks in other directions, thus yielding an Asymmetric ODF (AODF) field. We also show that a scalar measure of AODF asymmetry can be employed to obtain new contrast within the human brain. The feasibility of the technique is demonstrated on in vivo data obtained from the MGH-USC Human Connectome Project (HCP) and Parkinson's Progression Markers Initiative (PPMI) Project database. Characterizing asymmetry in neural tissue cytoarchitecture could be important for localizing and quantitatively assessing specific neuronal pathways.

1. Introduction The diffusional attenuation of MR signal has been exploited to observe tissue anisotropy through a characterization of the diffusion behavior of water molecules in human nervous system. In particular, several diffusion-weighted images with varying gradient directions enable determination of the orientation of coherently organized neural fiber tracts. When there is a single coherent tract per voxel, Diffusion Tensor Imaging (DTI) [1] has been shown to provide the dominant fiber direction [2, 3]. However, in heterogeneous fiber populations with crossing, kissing, fanning or bending fiber configurations, the unimodal Orientation Distribution Function (ODF) assumed by DTI is no longer capable of resolving the orientation of distinct fiber populations within a voxel. Numerous methods have been introduced [4-9] to compute the underlying ODF with less stringent assumptions than those employed in DTI. However, the assumption that the ODF profiles are antipodally symmetric has prevailed through much of the literature. The main purpose of this work is to overcome this assumption; doing so could provide hints as to how the cells are organized within a given voxel. The MR signal is predicted to be complex valued due to the Fourier

*

Corresponding author. E-mail address: [email protected] (G. Unal).

https://doi.org/10.1016/j.mri.2018.03.006 Received 8 August 2017; Accepted 8 March 2018 0730-725X/ © 2018 Elsevier Inc. All rights reserved.

transform relationship between the displacement distribution and the signal. Consequently, some techniques [10-12] have been formulated with the capability of handling complex-valued data to yield asymmetric profiles. However, obtaining accurate phase information in diffusion-weighted acquisitions is a formidable task particularly in in vivo acquisitions. The most important obstacle is the subject motion, which substantially distorts the phase of the signal [13]. Additionally, imperfections in the acquisition (e.g., imperfect B0 shift compensation), flow, and susceptibility variations within the tissue could all influence the detected signal phase substantially. Thus, attributing the observed phase shifts to diffusion alone would be very problematic. Therefore, current diffusion MRI processing pipelines employ magnitude-valued data. When magnitude-valued data are used, the regions containing Yshaped crossings or bending fibers might lead to symmetric profiles (e.g., unimodal or asterisk-shaped ODFs), even though the expected distribution of displacements for water molecules for Y-shaped crossings [10], or curving fibers [14] is asymmetric. To illustrate this point, we show in Fig. 1, the ODF field for a human brain obtained via Q-Ball Imaging (QBI) [4]. This image exhibits symmetric ODF components over a WM region, whose bending is grossly apparent, as well as regions

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

In order to model the sprouting or anti-symmetric crossing structures, Barmpoutis et al. [22] suggested an ODF field construction method by considering inter-voxel information, which creates a field of asymmetric spherical functions that differentiate between crossing and splaying fibers. Similarly, intervoxel information was utilized in Delputte et al. [23] and Duits and Franken [24] to create asymmetric voxel level representations. Reisert et al. [25] attempted to formalize the notion of an asymmetric fiber orientation distribution wherein asymmetric contributions to the FOD result from curvature of the fibers. They introduce a geometric interpretation, which leads to a continuity condition imposed in constrained spherical deconvolution to estimate asymmetric FODs. Similarly, Bastiani et al. [26] showed an asymmetric ODF model at the voxel level by using a spatial regularization in a spherical deconvolution framework [27]. In this paper, we present an efficient technique for voxel-by-voxel reconstruction of AODFs. We continue our preliminary study in [28] on the estimation of the AODFs based on a straightforward directional filtering. Specifically, we employ steerable filters on the ODF field, going beyond the basic Gaussian function-based filtering and bilateral filtering. Furthermore, we address the parameter selection issues, and quantify and demonstrate the results on various synthetic data sets as well as public MRI datasets (HCP [29] and PPMI [30]) including visualization and comparison of constructed maps of the presented asymmetry measure with those of the traditional fractional anisotropy measure. Rather than employing complex-valued data, we exploit the fact that for a voxel with asymmetric intravoxel geometry, the relevant features are expected to prevail only within the neighborhood along the direction of those features. As described in Section 2.1 below, such recovery is possible via a relatively direct filtering of the symmetric ODF field. The level of asymmetry is quantified via a scalar measure, which is presented in Section 2.2. Experimental results using both synthetic and real diffusion MRI data are presented to demonstrate the performance of the proposed algorithm in Section 3, followed by discussions and conclusions in Section 4 andSection 5, respectively.

Fig. 1. ODF field that is reconstructed by Q-ball imaging [4] on a sample HARDI data.

with potentially asymmetric junctions. Several works have been published with varying levels of relevance to our study. In a number of studies, antisymmetric fanning structures and degree of anisotropy in fanning were modeled with no reference to voxel-level asymmetric representation of fibers in [15-17]. Schultz [18] demonstrated that indicating likely fiber directions can be estimated from a small number of diffusion-weighted measurements by taking into account information from local neighborhoods. In Ref. [19], the challenge of partial volume averaging of fiber directions in diffusion MRI data due to multiple fiber orientations within a voxel is noted for tracts with high curvature, crossing, branching, bottlenecks and splaying or sprouting. To overcome the challenge, the data is regularized by considering all pairwise neighboring voxel ODFs in a local neighborhood weighted with an index of co-helicity. Pizzolato et al. [20] propose to recover the ensemble average propagator (EAP) directly from the complex diffusion weighted signal in order to also describe eventual diffusional asymmetry, thus obtaining an asymmetric EAP. Their results are limited to 3D Monte Carlo simulations of diffusion in white matter. In Ref. [21], a regularization on the ODF field is carried out using a cone model like we employ in this paper. However, the asymmetric ODF regularization employed by the authors is based on asymmetric weights assigned to the utilized pair of antipodally-symmetric cones placed at a given voxel during filtering, and no asymmetric ODF at a voxel level is constructed or depicted.

2. Materials and methods Our method relies on well-established symmetric ODF profiles where the reconstruction of the profiles is performed by following any of the classical HARDI techniques. Here, we utilized Diffusion Orientation Transform (DOT) [7] and Spherical Deconvolution (SD) [5, 9] using HARDI tools [31]. After the estimation of a multi-directional representation of the local fiber orientations (ODFs), we develop a method that exploits voxel-based ODFs in a conic spatial neighborhood to capture underlying asymmetry of the fiber populations in a voxel. The proposed method is depicted in a flowchart in Fig. 2. The ODF is defined at each voxel location x = (x,y,z) ∈ℝ3 by ϕ (x , u )̂ for the direction vector u ̂ ∈ 2 . At the given voxel x and the given direction u ,̂ a cone Cx , u ̂ (Fig. 3 (a)) is constructed by a certain radius r and height h. Although the directional asymmetry of the fiber pathway cannot be deduced from the ODF of a single voxel, it becomes apparent by exploiting the neighborhood diffusion profiles holistically,

Fig. 2. Flowchart of the proposed AODF reconstruction method.

146

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 3. (a) The cone Cx , u ̂ is constructed with a certain height and a radius; (b) an asymmetric ODF is reconstructed by rotating the cone model along all sampled directions on the unit sphere at each voxel point, and executing a directional regularization on the ODF field in Cx , u .̂

Fig. 6. (a) Asymmetry measure map; (b) gFA map calculated for circularly bending fibers: (dark red-to-yellow:high-to-low values).

steerable filter. 2.1.1. Gaussian function The idea of a cone-based directional ODF regularization is first realized by building a Gaussian function Gσ (x , u )̂ whose center lies on the axis of the cone Cx , u ̂ with a variance parameter σ [28]. At the voxel x, once the cones are constructed along each sampling direction u ̂ as in Fig. 3 (b), the weighted smoothing of the ODF field is carried out. The weights of Gaussian functions go from the highest on the central axis to the lowest on the side surface of the cone (Fig. 4). The regularization using a Gaussian function can be formulated by:

Fig. 4. lllustrative example of Gaussian weighted cones for directional regularization model.

and considering the ODFs at all voxels in the local volume delimited by the cone Cx , u .̂ 2.1. Asymmetric ODFs revealed by a cone model In this section, the directional regularization of ODF fields using the cone model described above, is presented through three different filtering approaches: (i) a Gaussian function; (ii) a bilateral filter; (iii) a

f (x , u )̂ =

Fig. 5. Illustration of steerable Gaussian basis filters: Top row shows the conic steerable basis filters ( corresponding basis filters.

147

1 Nx

∂ ∂x 2

ϕ (y, u )̂ Gσ (y, u )̂

∑

(1)

y ∈ Cx , u ̂

G (⋅) ,

∂ ∂y 2

G (⋅) ,

∂ ∂ G (⋅) ), ∂x ∂y

bottom row shows the cross-sectional views of

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 7. Circularly bending fibers (zoomed views of the same regions can be found below): (a) AODF field created by the analytical results in Özarslan et al. [14]; (b) ODF field created by the DOT method [7]; (c) AODF field created by the proposed technique.

(a)

(b)

(c)

Fig. 8. Visualizations of parameter sweeping for cone radius r, cone height h and σ. The heights of the graphs are the angle difference in degrees between peaks of the analytical ODFs and computed AODFs. Left: r, h parameter sweeping at σ = 0.5 for Gaussian filtering; Middle: r, h parameter sweeping at σ = 0.5 for steerable filtering based experiments; Right: σ sweeping experiments for Gaussian filtering (black curve), bilateral filtering (blue curve) and steerable filtering (red curve).

Table 1 Parameter setting experiments with peak angle errors (in degrees) between analytical ODFs and computed AODFs over the synthetic bending fiber data (mean ± standard deviation). The optimal r and h are estimated as 2mm, and σ is estimated as 0.50 for all filtering and ODF reconstruction methods. DOT Gaussian filtering

Bilateral filtering

Steerable filtering

Number Number Number Number Number Number Number Number Number

Spherical deconvolution (SD) of of of of of of of of of

gradients = 064,error(degrees) = 10.7 ± 2.3 gradients = 128,error(degrees) = 9.8 ± 3.2 gradients = 256,error(degrees) = 8.5 ± 2.9 gradients = 064,error(degrees) = 10.9 ± 3.9 gradients = 128,error(degrees) = 11 ± 3.2 gradients = 256,error(degrees) = 8.9 ± 2.9 gradients = 064,error(degrees) = 6.3 ± 2.5 gradients = 128,error(degrees) = 5.8 ± 3.2 gradients = 256,error(degrees) = 5.1 ± 2.9

∑ y ∈ Cx , u ̂

of of of of of of of of of

gradients = 064,error(degrees) = 9.6 ± 2.6 gradients = 128,error(degrees) = 8.7 ± 2.5 gradients = 256,error(degrees) = 8.0 ± 2.3 gradients = 064,error(degrees) = 10.3 ± 2.1 gradients = 128,error(degrees) = 9.3 ± 3.2 gradients = 256,error(degrees) = 8.2 ± 2.8 gradients = 064,error(degrees) = 5.1 ± 2.2 gradients = 128,error(degrees) = 4.9 ± 3.1 gradients = 256,error(degrees) = 4.0 ± 2.7

that is capable of seeing further along a given direction in a larger neighborhood. Moreover, the integration over only the corresponding orientation's ODF values leads to sharpening of the gross orientations while the non-essential (e.g., noise induced) or less dominant directions are suppressed. More importantly, this results in extraction of existing local asymmetries in a fiber distribution at the voxel level.

with

Nx =

Number Number Number Number Number Number Number Number Number

Gσ (y, u ), ̂ (2)

where f (x , u )̂ is the regularized AODF field value at voxel x along u .̂ Note that filtering is performed only along the orientation of the cone u ,̂ i.e., the ODF values only in the direction corresponding to that of the cone are smoothed out over the voxels in the cone. This creates a directional regularization effect and constructs an enhanced field of view

2.1.2. Bilateral function based on a Gaussian distribution Secondly, we exemplify the directional ODF regularization via a 148

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 9. Circularly bending fibers ODF field reconstructed by DOT [7]): (Top row) AODF field created by steerable filtering technique, (Middle row) AODF field created by bilateral filtering technique, (Bottom row) AODF field created by the Gaussian filtering for (a-d-g) r=2 mm, h=2 mm, (b-e-h) r=3 mm h=4mm, (c-f-i) r= 5mm, h= 5 mm parameter settings.

bilateral filter, which has been developed as a non-linear, edge-preserving low pass filter for images [32]. It involves an averaging of intensity values from nearby pixels to set the intensity value at each pixel in an image. We adapt the bilateral filtering to ODF fields as follows: for all voxels in a given cone, again the ODF values at the orientation corresponding to the axial direction of the cone are used as the intensity values in an image. The weighted average of those ODF values along the corresponding direction is calculated by:

f b (x , u )̂ =

1 Wx

∑

orientation-selective convolution kernel that can be expressed via a linear combination of a small set of rotated versions of itself. Steerable filters are popularly used for image enhancement and feature extraction [35-37]. We adapt the steerable filtering to ODF regularization by taking the ODF values of the voxels weighted by the kernel coefficients in the given cone. The linear combination of steerable Gaussian cone along the unit z direction can be defined as:

Sx , z =

ϕ (y, u )[ ̂ Gσs (‖x − y‖)

Gσr (‖ϕ (x , u )̂ − ϕ (y , u )‖)], ̂

∑ y ∈ Cx , u ̂

[Gσs (‖x − y‖) Gσr (‖ϕ (x , u )̂ − ϕ (y, u )‖)] ̂

b y (z )⋅Hy, z (5) ∂ ⎡ 2 G (y, ⎣ ∂x

(3)

∂ G (y, ∂y 2

∂ ∂ G (y, ∂x ∂y

z ), z ), z )⎤ is the function of where Hy, z = ⎦ Gaussian second order partial derivatives for y = [x,y,z] and represented as basis filters, and ⋅ is the dot product operator. The functions by(z) are trigonometric functions that interpolate the rotation of kernel along the z direction. In general, steerable cone filtering can be defined as:

where Gσr is the Gaussian range kernel function for smoothing differences in the ODF values and Gσs is the Gaussian spatial kernel function for smoothing differences in coordinates. Wx is the normalization term:

Wx =

∑ y ∈ Cx , z

y ∈ Cx , u ̂

. (4)

f s (x , u )̂ =

u )̂ is the AODF obtained via the directional bilateral filter of Thus, ODFs in the neighborhood of x. f b (x ,

∑

{Sx , u ̂ϕ (y, u )} ̂

, (6)

y ∈ Cx , u ̂

where {⋅} represents the selection of the steerable filter with the maximum response. Steerable Gaussian cone filters along the unit z direction are illustrated in Fig. 5. Top row illustrates the cone filters along ∂ the second derivative with respect to x ( 2 G (y, z ) ), second derivative

2.1.3. Steerable Gaussian function The orientation selective excitations in the human visual cortex have been exploited by early image processing techniques such as [33], which presented filters with isotropic and anisotropic components for restoration of noisy images. Anisotropic part of the filters involved a set of detectors at different orientations. Such early works constituted the origins of the idea of steerable filters in [34], which include an

∂x

with respect to y (

(

∂ ∂ G (y, ∂x ∂y

)

∂ G (y, ∂y 2

z ) ) and the derivative with respect to x and y

z ) . Bottom row shows the cross-sectional views of the cor-

responding filters. Linear combination of these filters represents the

149

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 10. Circularly bending fibers (ODF field reconstructed by spherical deconvolution Dell’Acqua et al. [9]): (Top row) AODF field created by steerable filtering technique, (Middle row) AODF field created by bilateral filtering technique, (Bottom row) AODF field created by the Gaussian filtering for (a-d-g) r=2 mm, h=2 mm, (b-e-h) r=3 mm h=4mm, (c-f-i) r= 5mm, h= 5 mm parameter settings.

other steerable filters.

a gray-scale map, which provides contrast dependent on the inherent intra-voxel asymmetry of fiber populations. A sample asymmetry map calculated for a circularly bending fiber configuration (Section 3.1.1) is shown in Fig. 6. In contrast to the most widely used scalar measure for WM fibers, the generalized Fractional Anisotropy (gFA) [4], which is not very informative here other than showing a high level of constant anisotropy, the asymmetry map signifies the changing level of asymmetry at different radii in concentric circular fibers.

2.2. Measuring asymmetry After we construct an AODF field using one of the methods in Section 2.1, an asymmetry measure can be computed at each voxel in order to create a scalar map that quantifies asymmetry over WM regions. For this purpose, we propose a measure [28] based on the angular similarity (covariance) metric in Özarslan et al. [11]. Envisioning functions to be vectors on an ∞-dimensional Hilbert space, the cosine of the angle between the two functions is a measure of their similarity. When symmetry is concerned, the functions can be taken to be the AODF denoted by f (u )̂ , and its reflection about the origin, f (−u )̂ . To this end, it is natural to consider the inner product to define a symmetry-seeking measure as in the following:

S

∫S

f (u )̂ 2du ̂

In this section, we will show the results of the cone-based directional ODF regularization and maps of the asymmetry measure on both synthetic and real data. 3.1. Synthetic fibers

̂ û ∫ f (u )̂ f (−u )d cos γ =

3. Experiments and results

=

l l ∑lmax ∑m =−l (−1)l |alm |2 =0 l max l ∑l = 0 ∑m =−l |alm |2

3.1.1. Synthetic bending fibers The first synthetic data is simulated for a full circular fiber geometry. Diffusion taking place in a curving fiber within a full circle geometry was considered in Özarslan et al. [14], wherein the physical signal attenuation for each portion of the curving circular fiber is derived. Using this result, we simulate the q-space (maximum q = 128mm−1, Δq = 1mm−1, number of points in the q-space = 128 × 128 × 128) complex signal for a set of concentric circularly bending fibers, which are reconstructed over a three-dimensional (3D) volumetric grid. A Fourier transform is taken to produce the true EAP.

. (7)

Here, the last expression provides the measure in terms of alm, which are the coefficients obtained when f (u )̂ is represented in a series of spherical harmonics. Then an asymmetry measure α can be defined at each voxel:

α = sin γ =

1 − cos2 γ .

(8)

This scalar measure can be computed on a voxel-by-voxel basis to create 150

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 11. Simulations for strong curvature fibers. (a): ODFs obtained via the method in Dell’Acqua et al. [9]; (b) AODF maps obtained by Gaussian filtering. (c) AODF maps obtained by bilateral filtering; (d) AODF maps obtained by steerable filtering.

underlying ODF field. First, via a sweeping of values for σ parameter, σ is set to 0.5 and kept fixed during all experiments. Fig. 8-(c) shows the σ-sweeping graph only for r = h = 2, and number of gradient directions = 256 case. The parameter sweeping for h and r at the selected σ value is demonstrated in Fig. 8-(a) and (b) for both Gaussian and steerable filtering. Inspecting this graph, the parameters cone height h and radius r are both selected as 2 mm (2 voxels) and fixed in all synthetic experiments. In Table 1, the optimal parameters (r, h, σ) for 64, 128, 256 number of gradient directions, and for each ODF reconstruction and filtering technique are summarized. In Fig. 9 and Fig. 10, the effects of various parameter settings on AODFs for DOT and SD reconstruction techniques are visualized.

Then, we use the DSI approach [6], where the EAP is integrated radially along given gradient directions to obtain the “true” ODFs, which display the inherent asymmetry in the fiber. One such sample simulated experiment is shown in Fig. 7. In~(a), the reconstructed ODFs for a set of circularly bending fibers using the just described DSI method is shown. In (b) the symmetric ODFs using one of the classical methods, more specifically, the DOT approach [7], are depicted. In (c) the AODFs at each voxel constructed by our method in Section 2.1.3 using r = h = 3mm are shown. Clearly, the proposed spatial directional regularization of the symmetric ODFs resulted in an asymmetric ODF field, which is naturally in line with the underlying curving geometry, as also verified by the “true” ODF field obtained from the DSI method. The bending structures, hence the asymmetry is clearly visible through the antipodally asymmetric peaks of the AODF model depicted at the voxel level in Fig. 7 (c). The asymmetry is decreasing for larger circles as expected, which is observed in Fig. 6 (a) that corresponds to the asymmetry map of the same example. Using the synthetic data for which the analytical ODFs are available, we carried out a set of experiments for selecting the parameters h, r, σ (defined in Section 2.1) for each filtering (Gaussian, bilateral, steerable) technique and ODF reconstruction method (SD, DOT) for different gradient directions. Particularly, we swept the parameters in the range of: 1 mm to 6 mm with 1 mm steps for h and r; 0.1 to 1.0 with 0.1 step for σ; with all possible combinations for gradient directions of 64, 128 and 256 [31]. In each experiment, we first estimate the peaks of analytical ODFs and AODFs. Then, at each voxel, we compute the angle difference between the analytical ODF and the AODF in degrees as the agreement criterion between the estimated AODF and the ‘ true’

3.1.2. Synthetic crossing fibers Three more synthetic examples, where the data are from the simulations of Descoteaux et al. [38], demonstrate the capabilities of the proposed technique. HARDI data was generated with b-value of 3000s/ mm2, N = 81 gradient directions, SNR = 35. The volumes are typically of size 30 × 30 × 18mm3 and they are representative of various fiber crossing scenarios including asymmetric junctions. In Fig. 11, a circular pattern of fibers meet a bending line of fibers in non-symmetric junctions. The ODF field shown in Fig. 11 (a) clearly displays the symmetry assumption imposed in conventional HARDI processing as seen from the visualized symmetric ODFs. The AODF field results are visualized for each technique: Gaussian (Fig. 11 (b)), bilateral (Fig. 11 (c)), and steerable (Fig. 11 (d)) filtering. Ellipses with matching colors show the comparison among each technique on the same fiber region of interest. In Fig. 12, a bottleneck or a kissing fiber geometry is depicted. 151

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 12. Simulations for kissing fibers. (a): ODFs obtained via the method in Dell’Acqua et al. [9]; (b) AODF maps obtained by Gaussian filtering; (c) AODF maps obtained by bilateral filtering; (d)AODF maps obtained by steerable filtering.

Fig. 13. Simulations for branching fibers. (a) ODFs obtained via the method in Özarslan et al. [7]; (b) AODF maps obtained by Gaussian filtering; (c) AODF maps obtained by bilateral filtering; (d) AODF maps obtained by steerable filtering.

152

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 14. AODFs constructed via steerable filtering on (a) PPMI and (b) HCP HARDI data.

Fig. 15. gFA (a) and asymmetry (b, c, d) maps (bright-to-dark:between 1 and 0) calculated for strong curvature fibers; (b) Filtering with a Gaussian function; (c) Bilateral filter based on a Gaussian distribution; (d) Steerable Gaussian filter.

Fig. 12 (a) shows the ODF representation of the kissing fibers whereas Fig. 12 (b) depicts the corresponding reconstructed asymmetric AODFs by Gaussian filtering. Fig. 12 (c,d) shows the reconstructed asymmetric AODFs by bilateral and steerable filtering, respectively. Noisy AODFs are depicted by yellow ellipses. On the right side of the figures, zoomed views of the regions are depicted. In order to demonstrate the synthetic data performance using the

DOT ODF reconstruction method, the last synthetic example represents the branching fibers with ODFs based on DOT method. Fig. 13 (a) depicts the ODF representation of the branching fibers while the corresponding reconstructed asymmetric AODFs by Gaussian filtering are shown in Fig. 13 (b). Fig. 13 (c,d) demonstrates the results of the reconstructed asymmetric AODFs by bilateral and steerable filtering, respectively. In areas with smaller crossing angles, the distinct peaks are 153

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 16. gFA (a) and asymmetry (b, c, d) maps (bright-to-dark:between 1 and 0) calculated for kissing fibers; (b) Filtering with a Gaussian function; (c) Bilateral filter based on a Gaussian distribution; (d) Steerable Gaussian filter.

on fiber pathways that exhibit high curvature, we selected the corpus callosum (CC) and uncinate fasciculus (UF). It is well known that noninvasive extraction of the visual pathways including CC and UF through Diffusion MRI is challenging due to the strong bending, crossing and kissing geometric patterns in the relevant anatomy. For the MRI data experiments, AODFs reconstructed by the steerable filter-based ODF regularization are utilized. In the proposed method, the radius parameter of the cone intuitively corresponds to the radius or thickness of the fiber bundles of interest. The parameters for r and h are selected as 4 mm depending on the average anatomical thickness of the bundles CC and UF [39, 40]. Fig. 14 shows the reconstructed AODFs for both the HCP and PPMI HARDI data for the CC fiber bundle.

not resolved in the ODF field. Reconstructing the AODFs alleviated this problem. 3.2. Real MRI data fibers HARDI data are obtained from two datasets: MGH-USC Human Connectome Project (HCP) (https://ida.loni.usc.edu/login.jsp)1 and Parkinson's Progression Markers Initiative (PPMI) Project database (http://www.ppmi-info.org). HCP data are acquired from 288 gradient directions with 1.25 × 1.25 × 1.25 mm3 voxel size whereas PPMI data are acquired from 64 gradient directions with 1.98 × 1.98 × 2 mm3 voxel size. The ODF field for both data are created by the DOT method [7]. As a preprocessing operation, we inserted a sigmoid function on the ODF field to enhance very high components and suppress very low components, which are likely to be due to noise. To focus

3.3. Asymmetry measure results Figs. 15, 16, and 17 demonstrate the comparison of gFA and asymmetry maps that are calculated for each filtering technique over a slice of the following: circular fibers crossing a slightly bending line fibers, kissing fibers and branching fibers. Figs. 18, 19 and 20 show results for white matter fibers CC and UF.

1

Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

154

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 17. gFA (a) and asymmetry (b, c, d) maps (bright-to-dark:between 1 and 0) calculated for kissing fibers; (b) Filtering with a Gaussian function; (c) Bilateral filter based on a Gaussian distribution; (d) Steerable Gaussian filter.

available in analytical form. This facilitated an inverse Fourier transform to obtain the EAPs, and then the reconstruction of the “true” ODF fields based on the DSI method. However, as well known, this approach based on the complex q-space signal could only be utilized in this special case, in contrast to real scenarios where the complex signal is unavailable or unreliable. The bending intravoxel fiber geometry, which cannot be captured by the symmetric ODFs, is captured for the true ODF as well as the AODF field obtained from the regularization. Moreover, as expected, the asymmetry is observed to visibly decrease from smaller to larger circles for both the true ODFs and the AODFs (Fig. 6). In the other synthetic experiments, the resulting representation of the local fiber orientations, particularly in splaying sections of the kissing fibers is clearly observed to exhibit the Y-pattern as desired. Steerable filtering performs the best among the filtering techniques as can be observed from Table 1 with its lowest peak angle error in different settings. This is due to its directional filtering through combinations of spatial derivatives of the Gaussian filter that yields the most sharply reconstructed AODF field, which also exhibits a visually smooth spatial variation. Furthermore, as the number of cone probing directions is increased, the peak angle errors decrease. This is expected due to an increased obtained resolution of fiber orientations with larger

In Fig. 19, the asymmetry and gFA maps for HCP data are given over an axial slice where CC fibers connecting the two brain hemispheres are distinctly visible. Fig. 20 depicts 3D visualizations of both the gFA and the asymmetry maps of Left UF and Right UF over a sample volume from the PPMI dataset.

4. Discussions The proposed algorithm is implemented in Matlab™ and C (Mex), and the experiments are run on a computer with Intel Processor Xeon™ X560 @ 2.67 GHz CPU computer of 64 GB memory. Our algorithm takes less than 12 min to create the AODF field based on Gaussian filtering for the whole brain HCP diffusion data with 288 gradient directions and takes less than 5 min for PPMI data with 64 gradient directions. Time doubles for bilateral filtering and it takes 2.5 times more for steerable filtering. As these filters are suitable for parallel implementations, the AODF directional regularization method can be potentially sped up for practical clinical use. In all experiments above, clearly the bending in the crossing was captured at the voxel-level representations of the junction regions after the directional regularization of the ODFs. We provided a validation using the circular bending fibers, whose complex q-space signal is 155

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 18. An axial slice depicting the asymmetry and gFA maps on a sample PPMI data in (a) and (b) respectively.

Fig. 19. An axial slice depicting the computed asymmetry and gFA maps on a HCP data sample in (a) and (b) respectively (bright-to-dark:between 1 and 0).

geometry of fibers is apparent. Besides asymmetry, our ODF regularization method offers another advantage: as SNR of the diffusion signal decreases, the appearance of spurious peaks as well as deviations from original fiber peak orientations increases as reported by [41]. The proposed technique naturally is expected to overcome noise and spurious peaks that exist in the ODFs through the introduced spatial regularization approach. The construction of AODFs is a step forward in delineating the subvoxel fiber architecture. The most significant improvement over traditional ODFs is seen in regions with bending, kissing, and splaying fibers. Though not illustrated here, employing AODFs in tractography could lead to an improvement in the results. In deterministic tractography based on symmetric ODFs, as the traced fiber enters a new voxel, it typically encounters a larger number of alternative directions

number of cone orientations sweeping diffusion data finely in a local neighborhood around the current voxel. The proposed method does not depend on any given fiber orientation distribution model, but only assumes that ODFs are estimated at the voxel level by any classical method such as SD or DOT. The appearance of the AODFs naturally depends on the method that is used to generate the initial ODFs, though our synthetic quantitative experiments, as in Table 1, demonstrate small variations in the peak angle errors between the two representative ODF methods. For the real data experiments, similarly, the reconstructed AODFs clearly display the inherent bending and asymmetry at the voxel level for the corresponding white matter tracts. Comparing the results (Fig. 1 vs Fig. 14), the unique representation capability of the AODF field against the ODFs for capturing an improved underlying intravoxel 156

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

Fig. 20. Comparison of gFA (Top-row) and asymmetry map (Bottom-row) for Left and Right Uncinate Fasciculus (UF) from PPMI Project database (dark red to yellow:high to low vaLues between 0 and 1)

Comparing the asymmetry map to the gFA map, note that the latter is less sensitive to bendings on the fiber bundles as seen in Fig. 18 and Fig. 19, particularly visible in the central parts of the CC with increased brightness in the asymmetry map while the gFA map is more homogeneous over the CC. For the UF bundle, the higher expressive quality in terms of bending in the asymmetry map compared to the gFA map is distinctively visible in 3D in Fig. 20, which is demonstrated for a subject from the PPMI database. Here the contrast in the gFA map is visibly different, which is most likely due to its insensitivity to curvature of the fibers.

along which the tracking could continue (see Fig. 12 (a) and (b) as an example). This is because every likely direction has its antipodally symmetric counterpart as another option. If AODFs are employed instead, fewer options are available; in the case of a T-junction or bending fibers, even the direction of entrance into the voxel is not available. Consequently, when tractography is performed, the directions of fibers as they enter and exit the voxel have a higher likelihood of being different when it is based on AODFs. In fact, the results could be markedly different when tractography methods that employ the same directional information over the entire voxel are employed. The difference is expected to be less if some spatial interpolation is performed during tractography though the exact nature of the differences would depend on the particular interpolation scheme employed. For the circularly bending fibers, the asymmetry measure is as expected rotationally invariant and increases from high to low radii. This feature is very clearly observed in Fig. 6 (a), whereas the gFA measure is nearly totally uniform in (b). Similarly, the proposed measure becomes larger in the asymmetric regions of interest: the circle itself as well as the crossing regions (Fig. 15), and in the splaying parts of the kissing fibers (Fig. 16), respectively. For the fibers bending slightly, the asymmetry map values are decreased compared to those more substantially bending, as desired (Fig. 15). In Fig. 16, the asymmetry values approach zero for the parts of the straight line of fibers that are away from the junction, as expected. The asymmetry values increase over the splaying fibers and the junctions. For the junction region, observe the expected fall in gFA measure. In the remaining regions, the gFA varies more homogeneously and less informatively in contrast to the asymmetry values. The results also show that the steerable filter performs better on all cases as it is observed to yield smoother results than the other approaches. For real data volumes, the asymmetry measure values are substantial for regions with fibers featuring relatively high curvature.

5. Conclusions In regions with bending, crossing, or sprouting fibers, cytoarchitecture within the voxel may exhibit asymmetry. For representing such inherent asymmetry, we employ spatial regularization through a cone model, which facilitates an inter-voxel filtering effect with a sharpening of the ODFs in some directions while suppressing peaks in other directions, thus yielding an AODF field. The feasibility of the techniques is demonstrated and compared quantitatively on a complex synthetic bending fiber simulation, and qualitatively on INRIA synthetic crossing fibers examples and in vivo data obtained from the HCP and PPMI Project databases. Experiments on synthetic geometries of circular, crossing, and kissing fibers show that the estimated AODFs successfully recover the asymmetry of the underlying geometry. AODFs reveal more information regarding the microstructure within the voxel compared to symmetric ODFs. The level of white matter voxel-level asymmetry is quantified via a scalar measure that is based on the AODFs. The asymmetry measure, which is observed to provide a variation with respect to intra-voxel curvature and underlying asymmetric crossings of fibers, has the potential to be utilized in analyzing white matter structures in population 157

Magnetic Resonance Imaging 49 (2018) 145–158

S. Cetin Karayumak et al.

studies complementing traditional diffusion anisotropy and diffusivity measures. This can lead to an improvement of the sensitivity and specificity of diffusion MRI to neurodegenerative diseases. Indeed, various neuroimaging studies describe pathology-induced white matter alterations and disruption of white matter integrity in presence of neurodegeneration (e.g. [42, 43]). Furthermore, a study by [44] stated the proportion of WM voxels containing crossing fibers to be between 63% and 90% of the whole WM. A significant proportion of those fiber populations are expected to make asymmetric crossings at an intra-voxel level due to current resolution levels of the diffusion MRI. Hence, as it is natural to expect an asymmetric micro-organization of fiber populations at the intra-voxel level, capturing the intra-voxel asymmetry and its disruptions due to biological or pathophysiological processes has a potential to correlate with neurodegenerative diseases. This is certainly a stimulating conjecture, which hopefully leads to future studies in this direction.

[21] [22] [23]

[24]

[25] [26]

[27]

[28]

References [1] Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophys J 1994;66:259–67. [2] Basser PJ, Mattiello J, LeBihan D. Estimation of the effective self-diffusion tensor from the NMR spin echo. J Magn Reson B 1994;103:247–54. [3] Basser PJ, Pajevic S, Pierpaoli C, Duda J, Aldroubi A. In vivo fiber tractography using DT-MRI data. Magn Reson Med 2000;44:625–32. [4] Tuch DS. Q-ball imaging. Magn Reson Med 2004;52:1358–72. [5] Tournier JD, Calamante F, Gadian DG, Connelly A. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage 2004;23:1176–85. [6] Wedeen VJ, Hagmann P, Tseng WYI, Reese TimothyG, Weisskoff RobertM. Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magn Reson Med 2005;54:1377–86. [7] Özarslan E, Shepherd TM, Vemuri BC, Blackband SJ, Mareci TH. Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT). NeuroImage 2006;31:1086–103. [8] Jian B, Vemuri BC, Özarslan E, Carney PR, Mareci TH. A novel tensor distribution model for the diffusion-weighted MR signal. NeuroImage 2007;37:164–76. [9] Dell’Acqua F, Scifo P, Rizzo G, Catani M, Simmons A, Scotti G, et al. A modified damped Richardson-Lucy algorithm to reduce isotropic background effects in spherical deconvolution. NeuroImage 2010;49:1446–58. [10] Liu C, Bammer R, Acar B, Moseley ME. Characterizing non-Gaussian diffusion by using generalized diffusion tensors. Magn Reson Med 2004;51:924–37. [11] Özarslan E, Koay CG, Shepherd TM, Komlosh ME, İrfanoğlu MO, Pierpaoli C. Mean apparent propagator (MAP) MRI: a novel diffusion imaging method for mapping tissue microstructure. NeuroImage 2013;78:16–32. [12] Ozcan A. Complete Fourier direct magnetic resonance imaging (CFD-MRI) for diffusion MRI. Front Integr Neurosci 2013;7:18. [13] Anderson AW, Gore JC. Analysis and correction of motion artifacts in diffusion weighted imaging. Magnetic Resonance in Medicine 1994;32:379–87. [14] Özarslan E, Koay CG, Basser PJ. Remarks on q-space MR propagator in partially restricted, axially-symmetric, and isotropic environments. Magn Reson Imaging 2009;27:834–44. [15] Nedjati-Gilani S, Alexander DC. Models for fanning and bending sub-voxel structures in diffusion MRI. Proceedings of the DMFC MICCAI 2009 Workshop. 2009. [16] Sotiropoulos S, Behrens T, Jbabdi S. Ball and rackets: inferring fiber fanning from diffusion-weighted {MRI}. NeuroImage 2012;60:1412–25. [17] Savadjiev P, Campbell JS, Descoteaux M, Deriche R, Pike GB, Siddiqi K. Labeling of ambiguous subvoxel fibre bundle configurations in high angular resolution diffusion {MRI}. NeuroImage 2008;41:58–68. [18] Schultz T. Towards resolving fiber crossings with higher order tensor inpainting. In: Laidlaw DH, Vilanova A, editors. New Developments in the Visualization and Processing of Tensor Fields. 2012. p. 253–65. [19] Campbell J, Savadjiev P, Siddiqi K, Pike G. Validation and regularization in diffusion MRI tractography. Int Symp Biomedical Imaging ISBI. 2006. p. 351–4. [20] Pizzolato M, Wassermann D, Boutelier T, Deriche R. Exploiting the phase in

[29] [30] [31]

[32] [33]

[34] [35]

[36] [37]

[38]

[39]

[40]

[41]

[42]

[43] [44]

158

diffusion MRI for microstructure recovery: towards axonal tortuosity via asymmetric diffusion processes. Springer International Publishing; 2015; pp. 109–116. Ehricke HH, Otto KM, Klose U. Regularization of bending and crossing white matter fibers in MRI q-ball fields. Magn Reson Imaging 2011;29:916–26. Barmpoutis A, Vemuri B, Howland D, Forder JR. Extracting tractosemas from a displacement probability field for tractography in DW-MRI. MICCAI 2008:9–16. Delputte S, Dierckx H, Fieremans E, DAsseler Y, Achten E, Lemahieu I. Postprocessing of brain white matter fiber orientation distribution functions. Int Symp Biomedical Imaging ISBI. 2007. p. 784–7. Duits R, Franken E. Left-invariant diffusions on the space of positions and orientations and their application to crossing-preserving smoothing of HARDI images. Int J Comput Vis 2011;92:231–64. Reisert M, Kellner E, Kiselev VG. About the geometry of asymmetric fiber orientation distributions. IEEE Trans Med Imag 2012;31:1240–9. Bastiani M, Cottaar M, Dikranian K, Aurobrata Ghosh, Zhang H, Alexander DC, et al. Improved tractography by modelling sub-voxel fibre patterns using asymmetric fibre orientation distributions. ISMRM- Int Soc Magnetic Resonance in Medicine. 2016. Bastiani M, Cottaar M, Dikranian K, Ghosh A, Zhang H, Alexander DC, Behrens TE, Jbabdi S, Sotiropoulos SN. Improved tractography using asymmetric fibre orientation distributions. NeuroImage 2017;158:205–18. Suheyla Cetin, Evren Ozarslan, Gozde Unal. Elucidating intravoxel geometry in diffusion-MRI: asymmetric orientation distribution functions (AODFs) revealed by a cone model. Medical Image Computing and Computer-Assisted Intervention MICCAI 2015: 18th International Conference, Munich, Germany, October 5–9, 2015, Proceedings, Part I. 2015. p. 231–8. Essen DCV, Smith SM, Barch DM, Behrens TE, Yacoub E, Ugurbil K, et al. The WU-Minn Human Connectome Project: an overview. Neuroimage 2013;80:262–79. PPMI. The Parkinson progression marker initiative (PPMI). Prog Neurobiol 2011;95:629–35. Canales-Rodríguez EJ, Melie-Garca L, Yasser Iturria-Medina, Yasser Alemn-Gmez. High angular resolution diffusion imaging (HARDI) tools. 2013 Available at http:// neuroimagen.es/webs/hardi_tools/. Tomasi C. Bilateral filtering for gray and color images. 1998. p. 839–46. Knutsson H, Wilson R, Granlund G. Anisotropic nonstationary image estimation and its applications: part I-restoration of noisy images. IEEE Trans. Commun. 1983;31:388–97. Freeman WT, Adelson EH. The design and use of steerable filters. IEEE Trans. Pattern Anal. Mach. Intell. 1991;13:891–906. Pietro Perona. Computer Vision – ECCV’92: Second European Conference on Computer Vision Santa Margherita Ligure, Italy, May 19–22, 1992 Proceedings. Berlin, Heidelberg: Springer Berlin Heidelberg; 1992. p. 3–18. Simoncelli EP, Farid H. Steerable wedge filters for local orientation analysis. IEEE Trans. Image Process. 1996;5:1377–82. Simoncelli EP, Freeman WT. The steerable pyramid: a flexible architecture for multi-scale derivative computation. ICIP - IEEE Int Conf Image Processing. 1995. p. 444–7. Descoteaux M, Deriche R, Knsche TR, Anwander A. Deterministic and probabilistic tractography based on complex fibre orientation distributions. IEEE Trans Med Imaging 2009;28:269–86. Ardekani BA, Bachman AH, Figarsky K, Sidtis JJ. Corpus callosum shape changes in early Alzheimer's disease: an MRI study using the OASIS brain database. Brain Struct Funct 2014;219:343–52. VonDerHeide R, Skipper L, Klobusicky E, Olson IR. Dissecting the uncinate fasciculus: disorders, controversies and a hypothesis. Brain: A Journal of Neurology 2013;136:1692–707. Parker G, Marshall D, Rosin P, Drage N, Richmond S, Jones D. A pitfall in the reconstruction of fibre ODFs using spherical deconvolution of diffusion MRI data. NeuroImage 2013;65:433–48. Kuceyeski A, Zhang Y, Raj A. Linking white matter integrity loss to associated cortical regions using structural connectivity information in Alzheimer's disease and fronto-temporal dementia: the loss in connectivity (loco) score. NeuroImage 2012;61:1311–23. Campanella M, Ius T, Skrap M, Fadiga L. Alterations in fiber pathways reveal brain tumor typology: a diffusion tractography study. PeerJ 2014;e497. Ben Jeurissen, Alexander Leemans, JacquesDonald Tournier, Jones DerekK, Jan Sijbers. Investigating the prevalence of complex fiber configurations in white matter tissue with diffusion magnetic resonance imaging. Hum Brain Mapp 2013;34:2747–66.