Fiber Orientation Estimation Using Nonlocal and Local Information

Abstract

Diffusion magnetic resonance imaging (dMRI) enables in vivo investigation of white matter tracts, where the estimation of fiber orientations (FOs) is a crucial step. Dictionary-based methods have been developed to compute FOs with a lower number of dMRI acquisitions. To reduce the effect of noise that is inherent in dMRI acquisitions, spatial consistency of FOs between neighbor voxels has been incorporated into dictionary-based methods. Because many fiber tracts are tube- or sheet-shaped, voxels belonging to the same tract could share similar FO configurations even when they are not adjacent to each other. Therefore, it is possible to use nonlocal information to improve the performance of FO estimation. In this work, we propose an FO estimation algorithm, Fiber Orientation Reconstruction using Nonlocal and Local Information (FORNLI), which adds nonlocal information to guide FO computation. The diffusion signals are represented by a set of fixed prolate tensors. For each voxel, we compare its patch-based diffusion profile with those of the voxels in a search range, and its nonlocal reference voxels are determined as the k nearest neighbors in terms of diffusion profiles. Then, FOs are estimated by iteratively solving weighted $\ell _{1}$-norm regularized least squares problems, where the weights are determined using local neighbor voxels and nonlocal reference voxels. These weights encourage FOs that are consistent with the local and nonlocal information. FORNLI was performed on simulated and real brain dMRI, which demonstrates the benefit of incorporating nonlocal information for FO estimation.

1 Introduction

By capturing the anisotropy of water diffusion in tissue, diffusion magnetic resonance imaging (dMRI) enables in vivo investigation of white matter tracts. The fiber orientation (FO) is a crucial feature computed from dMRI, which plays an important role in fiber tracking [5].

Voxelwise FO estimation methods have been proposed and widely applied, such as constrained spherical deconvolution [16], multi-tensor models [9, 13, 17], and ensemble average propagator methods [10]. In particular, to reduce the number of dMRI acquisitions required for resolving crossing fibers, sparsity assumption has been incorporated in the estimation problem. For example, it has been used in the multi-tensor framework [1, 9, 13], leading to dictionary-based FO estimation algorithms that have been shown to reconstruct FOs of good quality yet using a lower number of dMRI acquisitions [1].

Because of image noise that adversely affects FO estimation, the regularization of spatial consistency has been used in FO estimation problems. For example, smoothness of diffusion tensors and FOs has been used as regularization terms in the estimation in [12, 15], respectively, but no sparsity regularization is introduced. Other methods incorporate both sparsity and smoothness assumption. For example, in [11, 14] sparsity regularization is used together with the smoothness of diffusion images in a spherical ridgelets framework, where FO smoothness is enforced indirectly. More recently, [4, 18] manage to directly encode spatial consistency of FOs between neighbor voxels with sparsity regularization in the multi-tensor models by using weighted $\ell _{1}$-norm regularization, where FOs that are consistent with neighbors are encouraged. These methods have focused on the use of local information for robust FO estimation. However, because fiber tracts are usually tube-like or sheet-like [19], voxels that are not adjacent to each other can also share similar FO configurations. Thus, nonlocal information could further contribute to improved FO reconstruction by providing additional information.

In this work, we propose an FO estimation algorithm that improves estimation quality by incorporating both nonlocal and local information, which is named Fiber Orientation Reconstruction using Nonlocal and Local Information (FORNLI). We use a dictionary-based FO estimation framework, where the diffusion signals are represented by a tensor basis so that sparsity regularization can be readily incorporated. We design an objective function that consists of data fidelity terms and weighted $\ell _{1}$-norm regularization. The weights in the weighted $\ell _{1}$-norm encourage spatial consistency of FOs and are here encoded by both local neighbors and nonlocal reference voxels. To determine the nonlocal reference voxels for each voxel, we compare its patch-based diffusion profile with those of the voxels in a search range, and select the k nearest neighbors in terms of diffusion profiles. FOs are estimated by minimizing the objective function, where weighted $\ell _{1}$-norm regularized least squares problems are iteratively solved.

2 Methods

2.1 Background: A Signal Model with Sparsity and Smoothness Regularization

Sparsity regularization has been shown to improve FO estimation and reduce the number of gradient directions required for resolving crossing fibers [1]. A commonly used strategy to incorporate sparsity is to model the diffusion signals using a fixed basis. The prolate tensors have been a popular choice because of their explicit relationship with FOs [1, 9, 13]. Specifically, let $\{\mathbf {D}_{i}\}_{i=1}^{N}$ be a set of N fixed prolate tensors. The primary eigenvector (PEV) $\varvec{v}_{i}$ of each $\mathbf {D}_{i}$ represents a possible FO and these PEVs are evenly distributed on the unit sphere. The eigenvalues of the basis tensors can be determined by examining the diffusion tensors in noncrossing tracts [9]. Then, the diffusion weighted signal $S_{m}(\varvec{g}_{k})$ at voxel m associated with the gradient direction $\varvec{g}_{k}$ ($k = 1,2,\ldots ,K$) and b-value $b_{k}$ can be represented as

$$\begin{aligned} S_{m}(\varvec{g}_{k}) = S_{m}(\varvec{0})\sum \limits _{i=1}^{N}f_{m,i}e^{-b_{k}\varvec{g}_{k}^{T}\mathbf {D}_{i}\varvec{g}_{k}} + n_{m}(\varvec{g}_{k}), \end{aligned}$$

(1)

where $S_{m}(\varvec{0})$ is the baseline signal without diffusion weighting, $f_{m,i}$ is $\mathbf {D}_{i}$’s unknown nonnegative mixture fraction ($\sum _{i=1}^{N} f_{m,i} = 1$), and $n_{m}(\varvec{g}_{k})$ is noise.

We define $y_{m}(\varvec{g}_{k}) = S_{m}(\varvec{g}_{k})/S_{m}(\varvec{0})$ and $\eta _{m}(\varvec{g}_{k}) = n_{m}(\varvec{g}_{k})/S_{m}(\varvec{0})$, and let $\varvec{y}_{m}=(y_{m}({\varvec{g}_1}),y_{m}({\varvec{g}_2}),\ldots ,y_{m}({\varvec{g}_K}))^{T}$ and $\varvec{\eta }_{m}=(\eta _{m}({\varvec{g}_1}),\eta _{m}({\varvec{g}_2}),\ldots ,\eta _{m}({\varvec{g}_K}))^{T}$. Then, Eq. (1) can be written as

$$\begin{aligned} \varvec{y}_{m} = \mathbf {G}\varvec{f}_{m} + \varvec{\eta }_{m}, \end{aligned}$$

(2)

where $\mathbf {G}$ is a $K\times N$ dictionary matrix with $G_{ki}=e^{-b_{k}\varvec{q}_{k}^{T}\mathbf {D}_{i}\varvec{q}_{k}}$, and $\varvec{f}_{m}=(f_{m,1},f_{m,2},\ldots ,f_{m,N})^{T}$. Based on the assumption that at each voxel the number of FOs is small with respect to the number of gradient directions, the mixture fractions can be estimated using a voxelwise sparse reconstruction formulation

$$\begin{aligned} \hat{\varvec{f}}_{m} = \mathop {{{\mathrm{arg\,min}}}}\limits _{\varvec{f}_{m}\ge \varvec{0},||\varvec{f}_{m}||_{1}=1}||\mathbf {G}\varvec{f}_{m}-\varvec{y}_{m}||_{2}^{2} + \beta ||\varvec{f}_{m}||_{0}. \end{aligned}$$

(3)

In practice, the constraint of $||\varvec{f}_{m}||_{1}=1$ is usually relaxed, and the sparse reconstruction can be either solved directly [8] or by approximating the $\ell _{0}$-norm with $\ell _{1}$-norm [1, 9, 13]. Basis directions corresponding to nonzero mixture fractions are determined as FOs.

To further incorporate spatial coherence of FOs, weighted $\ell _{1}$-norm regularization has been introduced into dictionary-based FO estimation [4, 18]. For example, in [18] FOs in all voxels are jointly estimated by solving

$$\begin{aligned} \{\hat{\varvec{f}}_{m}\}_{m=1}^{M}=\mathop {{{\mathrm{arg\,min}}}}\limits _{\varvec{f}_{1},\varvec{f}_{2},\ldots ,\varvec{f}_{M}\ge \varvec{0}}\sum \limits _{m=1}^{M}||\mathbf {G}\varvec{f}_{m}-\varvec{y}_{m}||_{2}^{2} + \beta ||\mathbf {C}_{m}\varvec{f}_{m}||_{1}, \end{aligned}$$

(4)

where M is the number of voxels and $\mathbf {C}_{m}$ is a diagonal matrix that encodes neighbor interaction. It places smaller penalties on mixtures fractions associated with basis directions that are more consistent with neighbor FOs so that these mixture fractions are more likely to be positive and their associated basis directions are thus encouraged.

2.2 FO Estimation Incorporating Nonlocal Information

In image denoising or segmentation problems, nonlocal information has been used to improve the performance [3, 6]. In FO estimation, because fiber tracts are usually tube-shaped (e.g., the cingulum bundle) or sheet-shaped (e.g., the corpus callosum) [19], voxels that are not adjacent to each other can still have similar FO patterns, and it is possible to use nonlocal information to improve the estimation. We choose to use a weighted $\ell _{1}$-norm regularized FO estimation framework similar to Eq. (4), and encode the weighting matrix $\mathbf {C}_{m}$ using both nonlocal and local information.

Finding Nonlocal Reference Voxels. For each voxel m, the nonlocal information is extracted from a set $\mathcal {R}_{m}$ of voxels, which are called nonlocal reference voxels and should have diffusion profiles similar to that of m. To identify the nonlocal reference voxels for m, we compute patch-based dissimilarities between the voxel m and the voxels in a search range $\mathcal {S}_{m}$, like the common practice in nonlocal image processing [3, 6]. Specifically, we choose a search range of a $11\times 11\times 11$ cube [3] whose center is m. The patch at each voxel $n\in \mathcal {S}_{m}$ is formed by the diffusion tensors of its 6-connected neighbors and the diffusion tensor at n, which is represented as $\mathbf {\Delta }_{n}=(\mathbf {\Delta }_{n,1},\ldots ,\mathbf {\Delta }_{n,7})$.

We define the following patch-based diffusion dissimilarity between two voxels m and n

$$\begin{aligned} d_{\mathbf {\Delta }}(\mathbf {\Delta }_{m},\mathbf {\Delta }_{n}) = \frac{1}{7}\sum \limits _{j=1}^{7} d(\mathbf {\Delta }_{m,j},\mathbf {\Delta }_{n,j}), \end{aligned}$$

(5)

where $d(\cdot ,\cdot )$ is the log-Euclidean tensor distance [2]

$$\begin{aligned} d(\mathbf {\Delta }_{m,j},\mathbf {\Delta }_{n,j}) = \sqrt{\mathrm {Trace}(\{\log (\mathbf {\Delta }_{m,j})-\log (\mathbf {\Delta }_{n,j})\}^{2})}. \end{aligned}$$

(6)

For each m we find its k nearest neighbors in terms of the diffusion dissimilarity in Eq. (5), and define them as the nonlocal reference voxels. k is a parameter to be specified by users. Note that although we call these reference voxels nonlocal, it is possible that $\mathcal {R}_{m}$ contains the neighbors of m as well, if they have very similar diffusion profiles to that of m. We used the implementation of k nearest neighbors in the scikit-learn toolkit¹ based on a ball tree search algorithm.

Guided FO Estimation. We seek to guide FO estimation using the local neighbors and nonlocal reference voxels. Like [18], we use a 26-connected neighborhood $\mathcal {N}_{m}$ of m. Then, the set of voxels guiding FO estimation at m is $\mathcal {G}_{m}=\mathcal {N}_{m}\cup \mathcal {R}_{m}$.

Using $\mathcal {G}_{m}$, we extract a set of likely FOs for m to determine the weighting of basis directions and guide FO estimation. First, a voxel similarity between m and each voxel $n\in \mathcal {G}_{m}$ is defined

$$\begin{aligned} w(m,n)=\left\{ {\begin{array}{*{20}l} {\exp \{-\mu d^{2}(\mathbf {D}_{m},\mathbf {D}_{n})\}, \,\,\mathrm {if}\,\,n\in \mathcal {N}_{m}}\\ {\exp \{-\mu d^{2}_{\mathbf {\Delta }}(\mathbf {\Delta }_{m},\mathbf {\Delta }_{n})\},\,\,\mathrm {otherwise}}\\ \end{array} } \right. , \end{aligned}$$

(7)

where $\mu =3.0$ is a constant [18], and $\mathbf {D}_{m}$ and $\mathbf {D}_{n}$ are the diffusion tensors at m and n, respectively. When n is a neighbor of m, the voxel similarity is exactly the one defined in [18]; when n is not adjacent to m, the voxel similarity is defined using the patches $\mathbf {\Delta }_{m}$ and $\mathbf {\Delta }_{n}$. Second, suppose the FOs at a voxel n are $\{\varvec{w}_{n,j}\}_{j=1}^{W_{n}}$, where $W_{n}$ is the number of FOs at n. For each m we can compute the similarity between the basis direction $\varvec{v}_{i}$ and the FO configurations of the voxels in the guiding set $\mathcal {G}_{m}$

$$\begin{aligned} R_{m}(i) = \sum \limits _{n\in \mathcal {G}_{m}} w(m,n)\max \limits _{j=1,2,\ldots ,W_n} |\varvec{v}_{i}\cdot \varvec{w}_{n,j}|, \quad i=1,2,\ldots ,N. \end{aligned}$$

(8)

When $\varvec{v}_{i}$ is aligned with the FOs in many voxels in the guiding set $\mathcal {G}_{m}$ and these voxels are similar to m, large $R_{m}(i)$ is observed, indicating that $\varvec{v}_{i}$ is likely to be an FO. Note that $R_{m}(i)$ is similar to the aggregate basis-neighbor similarity defined in [18]. Here we have replaced the neighborhood $\mathcal {N}_{m}$ in [18] with the guiding set $\mathcal {G}_{m}$ containing both local and nonlocal information. These $R_{m}(i)$ can then be plotted on the unit sphere according to their associated basis directions, and the basis directions with local maximal $R_{m}(i)$ are determined as likely FOs $\mathcal {U}_{m}=\{\varvec{u}_{m,p}\}_{p=1}^{U_{m}}$ ($U_{m}$ is the cardinality of $\mathcal {U}_{m}$) at m [18].

With the likely FOs $\mathcal {U}_{m}$, the diagonal entries of $\mathbf {C}_{m}$ are specified as [18]

(9)

where $\alpha$ is a constant controlling the influence of guiding voxels. Smaller weights are associated with basis directions closer to likely FOs, and these directions are encouraged. In this work, we set $\alpha =0.8$ as suggested by [18].

We estimate FOs in all voxels by minimizing the following objective function with weighted $\ell _{1}$-norm regularization,

$$\begin{aligned} E(\varvec{f}_{1},\varvec{f}_{2},\ldots ,\varvec{f}_{M}) = \sum \limits _{m=1}^{M}||\mathbf G \varvec{f}_{m}-\varvec{y}_{m}||_{2}^{2} + \frac{\beta }{W_{m}} ||\mathbf C _{m}\varvec{f}_{m}||_{1}, \end{aligned}$$

(10)

where $\varvec{f}_{m}\ge \varvec{0}$ and $\beta$ is a constant. Note that we assign smaller weights to the weighted $\ell _{1}$-norm when the number of FOs is larger, which in practice increases accuracy. In this work, we set $\beta =0.3$, which is smaller than the one used in [18] because the number of gradient directions in the dMRI data is smaller than that in [18]. Because $\mathbf C _{m}$ is a function of the unknown FOs, to solve Eq. (10) we iteratively solve $\varvec{f}_{m}$ sequentially. At iteration t, for each $\varvec{f}_{m}$ we have

$$\begin{aligned} \hat{\varvec{f}}^{t}_{m}= & {} \mathop {{{\mathrm{arg\,min}}}}\limits _{\varvec{f}_{m} \ge \varvec{0}} E(\hat{\varvec{f}}^{t}_{1},\ldots ,\hat{\varvec{f}}^{t}_{m-1},\varvec{f}_{m}, \hat{\varvec{f}}^{t-1}_{m+1},\ldots ,\hat{\varvec{f}}^{t-1}_{M}) \nonumber \\= & {} \mathop {{{\mathrm{arg\,min}}}}\limits _{\varvec{f}_{m}\ge \varvec{0}}||\mathbf G \varvec{f}_{m}-\varvec{y}_{m}||_{2}^{2} + \frac{\beta }{W^{t-1}_{m}} ||\mathbf C _{m}^{t}\varvec{f}_{m}||_{1}, \end{aligned}$$

(11)

which is a weighted Lasso problem that can be solved using the strategy in [17].

3 Results

3.1 3D Digital Crossing Phantom

A 3D digital phantom (see Fig. 1) with the same tract geometries and diffusion properties used in [18] was created to simulate five tracts. Thirty gradient directions ($b=1000\,\text {s}/\text {mm}^{2}$) were used to simulate the diffusion weighted images (DWIs). Rician noise was added to the DWIs. The signal-to-noise ratio (SNR) is 20 on the b0 image.

Fig. 1.

3D rendering of the digital phantom.

Fig. 2.

FO estimation errors. (a) Means and standard deviations of the FO errors of CSD, CFARI, FORNI, and FORNLI; (b) mean FORNLI FO errors with different numbers of nonlocal reference voxels in regions with noncrossing or crossing tracts.

FORNLI with $k=4$ was applied on the phantom and compared with CSD [16], CFARI [9], and FORNI [18] using the FO error proposed in [18]. CSD and CFARI are voxelwise FO estimation methods, and FORNI incorporates neighbor information for FO estimation. We used the CSD implementation in the Dipy software², and implemented CFARI and FORNI using the parameters reported in [9, 18], respectively. The errors over the entire phantom and in the regions with noncrossing or crossing tracts are plotted in Fig. 2(a), where FORNLI achieves the most accurate result. In addition, we compared the two best algorithms here, FORNI and FORNLI, using a paired Student’s t-test. In all four cases, errors of FORNLI are significantly smaller than those of FORNI ($p<0.05$), and the effect sizes (Cohen’s d) are between 0.5 and 0.6.

Next, we studied the impact of the number of nonlocal reference voxels. Using different k, the errors in regions with noncrossing or crossing tracts are shown in Fig. 2(b). Note that $k=0$ represent cases where only the local information from neighbors is used. Incorporation of nonlocal information improves the estimation quality, especially in the more complex regions with three crossing tracts. When k reaches four, the estimation accuracy becomes stable, so we will use $k=4$ for the brain dMRI dataset.

3.2 Brain dMRI

We selected a random subject in the publicly available dataset of COBRE [7]. The DWIs and b0 images were acquired on a 3T Siemens Trio scanner, where 30 gradient directions ($b=800~\mathrm {s}/\mathrm {mm}^{2}$) were used. The resolution is 2 mm isotropic. The SNR is about 20 on the b0 image.

To evaluate FORNLI (with $k=4$) and compare it with CSD, CFARI, and FORNI, we demonstrate the results in a region containing the crossing of the corpus callosum (CC) and the superior longitudinal fasciculus (SLF) in Fig. 3. We have also shown the results of FORNLI with $k=0$, where no nonlocal information is used. By enforcing spatial consistency of FOs, FORNI and FORNLI improve the estimation of crossing FOs. In addition, in the orange box FORNLI ($k=4$) achieves more consistent FO configurations than FORNI; and in the blue box, compared with FORNI and FORNLI ($k=0$), FORNLI ($k=4$) avoids the FO configurations in the upper-right voxels that seem to contradict with the adjacent voxels by having sharp turning angles.

Fig. 3.

FO estimation in the crossing regions of SLF and CC overlaid on the fractional anisotropy map. Note the highlighted region for comparison.

4 Conclusion

We have presented an FO estimation algorithm FORNLI which is guided by both local and nonlocal information. Results on simulated and real brain dMRI data demonstrate the benefit of the incorporation of nonlocal information for FO estimation.