The study of neural connectivity has grown rapidly in the past decade. Revealing brain anatomical connection improves not only clinical measures but also cognition understanding. In order to achieve this goal, we have to track neural fiber pathways first. Aiming to estimate 3D fiber pathways more accurately from orientation distribution function (ODF) fields, we presented a novel tracking method based on nonuniform rational B-splines (NURBS) curve fitting. First, we constructed ODF fields from high angular resolution diffusion imaging (HARDI) datasets using diffusion orientation transform (DOT) method. Second, under the angular and length constraints, the consecutive diffusion directions were extracted along each fiber pathway starting from a seed voxel. Finally, after the coordinates of the control points and their corresponding weights were determined, NURBS curve fitting was employed to track fiber pathways. The performance of the proposal has been evaluated on the tractometer phantom and real brain datasets. Based on several measure metrics, the resulting fiber pathways show promising anatomic consistency.

1. Introduction

An outstanding characteristic of white matter (WM) is its fibrillar construction. WM consists of tightly packed and coherently aligned axons that are surrounded by glial cells and that often are organized in bundles [1]. Axons are protected by myelin sheaths, which restricts the free diffusion of water molecules. As a result, the micrometric movements of water molecules are hindered to a greater extent in a direction perpendicular to the axonal orientation than parallel to it. It is now generally accepted that microscopic boundaries to diffusion in WM coincide with the local orientations of WM fiber pathways [24]. With this feature, we can trace fiber pathways and then reveal anatomical connection between brain functional areas.

Compared to diffusion tensor imaging (DTI), high angular resolution diffusion imaging (HARDI) could resolve multiple intravoxel fiber orientations contained in a WM voxel. Moreover, HARDI just needs to sample the diffusion signal on a spherical shell as opposed to a complete three-dimensional Cartesian grid of DSI [57]. At present, there are numerous tracking methods based on HARDI, which could be classified into deterministic and probabilistic algorithms [8]. They exploit the diffusion anisotropy to follow fiber tracts from voxel to voxel through the brain [9]. Recently, multishell multitissue (MSMT) models have been proposed to deal with partial volume effects and can remarkably increase the precision of fiber orientations over single-shell models [10].

Streamline tracking is an important deterministic approach. Streamline tracking propagates paths within the vector field of local fiber orientations [9], providing deterministic connectivity information between different brain functional areas. Later, many variants of the streamline method have been presented. The streamline-based tracking technique is the one most commonly used in tractography, and it appears to give excellent results in many instances if the vector field is smooth and the fibers are strongly oriented along a certain direction. However, the major drawback of streamline-based methods is that the estimation error accumulates along the tracking length [11, 12]. However, the partial volume effects such as crossing, kissing, merging, and splitting in imaging voxels increase the complexity in streamline tracking.

There are also some nonstreamline tractography algorithms. In the graph-based method, each voxel is treated as a graph node, and graph arcs connect neighboring voxels. The weights assigned to each arc are the representative of both structural and diffusivity features [13]. When partial volume exists, the algorithm treats the image as a multigraph and distributes the connectivities in a weighted manner. Aranda et al. presented a particle method which was proposed to estimate fiber pathways from multiple intravoxel diffusion orientations (MIVDO) [14]. The process starts with the definition of a point in WM region in which a virtual particle is allocated. The particle is iteratively moved along the local diffusion orientations until a stopping criterion is met. The estimation of fiber pathways is determined by the particle trajectory. Galinsky and Frank proposed a method for estimating local diffusion and fiber tracts based upon the information entropy flow that computes the maximum entropy trajectories [15]. This novel approach to fiber tracking incorporates global information about multiple fiber crossings in each individual voxel. Malcolm et al. used Watson function to analyze ODF construction, which provides a compact representation of the diffusion-anisotropic signal [16]. This algorithm models the diffusion as a discrete mixture of Watson directional functions and performs tractography within a filtering framework. Recently, global tractography was proposed in [17], which aims to find the full track configuration that best explains the measured diffusion weighted imaging (DWI) data. This data-driven approach was reported that it could improve valid neural connection rate over streamline methods.

The other classes are probabilistic approaches. This class of methods utilizes a stochastic process to estimate the connection probability between brain areas. A Bayesian approach was presented in [18], and it handled noise in a theoretically justified way. The persistent angular structure (PAS) of fiber bundles was used to drive probabilistic tracts, and PDF is incorporated into the method to estimate the whole-brain probability maps of anatomical connection [19]. Using automatic relevance determination in a Bayesian estimation scheme, the tracking in a multivector field was performed with significant advantages in sensitivity [20]. The residual bootstrap method made use of spherical harmonic (SH) representation for HARDI data in order to estimate the uncertainty in multimodal q-ball reconstructions [21]. However, these methods cannot directly delineate the fiber paths in 3D brain space. Furthermore, they are very time consuming in resolving the complexity of the diffusion pattern within each HARDI voxel.

In [22, 23], the authors argued that NURBS provides a framework to characterize WM pathways. However, the determination of the parameters including control points and weights has not been discussed. This paper has comprehensively explored the tracking method based on NURBS curve fitting and has detailed how to determine the related parameters. The tracking method consists of three steps: first is the computation of ODF field from HARDI datasets; second is the selection of consecutive diffusion directions along a fiber pathway; and the last is NURBS pathway fitting. This method was evaluated on tractometer phantom and real brain datasets.

2. Materials and Methods

2.1. HARDI Datasets

Two different types of HARDI datasets are used to evaluate our approach: from the physical diffusion phantom of tractometer and from an in vivo human brain. For each dataset, we firstly constructed ODF fields using DOT method [24] and then applied the proposed algorithms to estimate fiber paths.

Phantom study was performed using data acquired from a physical diffusion phantom of tractometer. Imaging parameters for the 3 × 3 × 3 mm acquisition were as follows: field of view FOV = 19.2 cm, matrix 64 × 64, slice thickness TH = 3 mm, read bandwidth RBW = 1775 Hz/pixel, partial Fourier factor 6/8, parallel reduction factor GRAPPA = 2, repetition time TR = 5 s, and echo times TE = 102 ms. A SNR of 15.8 was measured for the baseline (b = 0 s/mm2) image. SNR of HARDI at b-values = 2000 s/mm2 were evaluated. The diffusion sensitization was applied along a set of 64 orientations uniformly distributed over the sphere [25]. For comparative study, the ground truth fibers are available on the website http://www.lnao.fr/spip.php?rubrique79 [25].

A healthy volunteer was scanned on a Siemens Trio 3T scanner with 12 channel coils. The acquisition parameters were as follows: two images with b = 0 s/mm2, 64 DW images with unique, and isotropically distributed orientations (b = 2000 s/mm2). TR = 6700 ms, TE = 85 ms, and voxel dimensions equal to 2 × 2 × 2 mm. The SNR is, approximately, equal to 36.

2.2. ODF Fields

Compared with diffusion tensor, ODFs reflect the diffusion probability along any given angular direction, and higher values indicate higher consistency between the fiber orientation and diffusion direction. ODFs can be seen as a continuous function over the sphere that encodes diffusion anisotropy of water molecules within each voxel. There are two definitions of ODF. One is Tuch’s nonmarginal ODF that is defined as the radial integration of PDF and does not represent a true probability density [26, 27]. The other is marginal ODF that is introduced by Wedeen, and it is a true probability density since its integral over the sphere is one [28]. ODF peaks are assumed to correspond to the underlying fiber orientations. At present, there are several algorithms to compute ODFs from HARDI datasets. Tuch presented a simple linear matrix formulation that was provided to construct ODFs using radial basic function (RBF) [26]. Diffusion orientation transform (DOT) converts water diffusivity profiles into probability profiles under the monoexponential signal decay assumption through computing PDF at a fixed distance from the origin [24, 29, 30]. Spherical deconvolution (SD) estimates fiber orientations by assuming that a single response function can adequately describe HARDI signals measured from any fiber bundle [31]. Compared to other methods, DOT can improve the angular resolution, make the ODF sharper, and keep its accuracy and robustness to noise [27, 30]. In our work, we used DOT to construct ODFs from HARDI datasets.

After ODF fields were constructed, we detected ODF local maxima by thresholding over the sampling shell. Only those above ODF mean value would be retained. This operation can avoid the noise interference effectively [28]. Finally, ODF fields are transformed into vector fields, and we can describe a voxel using a matrix containing diffusion vectors and its corresponding diffusion probability.

The term denotes a diffusion direction, and is the diffusion probability along this orientation. In the next section, we would use this matrix to compute the control points and weights for NURBS pathway fitting.

2.3. Diffusion Directions along a Fiber Pathway

Before we conduct NURBS tracking, the consecutive directions along the same pathway have to be extracted. The orientations of fiber populations within a voxel coincide with the local maxima of ODFs [28]. ODF value along a direction is the reflection of diffusion probability of all the water molecules in a voxel, so it is reasonable to assume that the diffusion directions always pass through the voxel center. The aim of this step is to find the consecutive directions among the neighbors of a seed voxel. Here, we presented a new algorithm to achieve the goal. For the sake of simplicity, we used a two-dimensional diagram as an example to illustrate the process, shown as Figure 1(a). Compared to FACT algorithm [32], it can improve the extraction accuracy of discrete consecutive directions along a pathway. As we can see from Figure 1(b), in FACT, an unreasonable path was found (marked by red dashed lines). But if the distance between V1 (blue line in the seed voxel) and the center points of its neighbor voxel is considered here, we could get a more reasonable pathway (marked by blue dashed lines in Figure 1(b)). The algorithm is summarized as Algorithm 1. The input parameters, including fiber length threshold , angle threshold , and fractional anisotropy (FA) threshold should be determined according to actual situation.

Input: ROI, , , .
Output: Consecutive directions along a pathway.
(1) Select one seed voxel in ROI.
(2) Choose one direction in the seed voxel, and its linear equation could be written as follows.
(3) Calculate the distance between and the center of 26-connected neighbor voxels according to the equation given below.
(4) If , the neighbor voxel is considered as a candidate voxel.
(5) The angles between and the direction in the candidate voxel is calculated.
(6) If , is preserved as a consecutive direction. If no consecutive direction is obtained, save current pathway and stop tracking, then turn to (1). Otherwise, go to (7).
(7) , (2)∼(6) are repeated until .
2.4. NURBS Fitting

NURBS is a powerful tool to describe complex curves using a small number of parameters. It is a wonderful modeling method of curves and can control the object more conveniently and efficiently than traditional modeling method [33]. The order of a NURBS curve defines the number of nearby control points that could influence any given point on the curve. In practice, cubic curves are the ones most commonly used. Higher order curves are seldom used because they may lead to internal numerical problems and require disproportionately large computation time [3436]. The number of control points must be greater than or equal to the order of the curve. In this work, we traced nerve fiber pathways based on NURBS curve fitting. In the fitting, the parameters including control points and weights are needed. The consecutive directions were used to compute control points. The weights were computed according to . In NURBS tracking, we could use both control points and weights to hold local shape control of fiber pathways. We present two tracking methods based on NURBS according the fitting rule, including general NURBS fitting (NURBS-G) and tangent NURBS fitting (NURBS-T). The whole procedure of NURBS tracking is shown in Figure 2.

2.5. NURBS-T

A fiber pathway can be considered as a 3D curve, and its local tangent vector is consistent with the diffusion orientation [37]. According to this premise, we presented NURBS-T algorithm to trace fiber paths. To make it easier to explain, the 2D tracking process is illustrated in Figure 3. The algorithm is outlined in Algorithm 2.

Input: consecutive diffusion directions
Output: fiber pathways
(1) Determine the control points. Because that the diffusion direction (denoted by red thin line in Figure 3) is tangent to the fiber path, there should be three control points that situate on the same direction [33], including the center point (blue dot in Figure 3) of the voxel and the two intersection points (yellow dot in Figure 3) between the diffusion direction and the facets of the voxel. The intersection points could be obtained according to the equation given below:
where , , and are the length, width, and height of the voxel, respectively. is the coordinate of the intersection point. Finally, we would get a series of points.
where and are the intersection points between the diffusion direction and the facets of the path-through voxel. is the center point of the voxel.
(2) Calculate the weights. The weight indicates the attraction of the control points to a path, and we can locally modify the path by adjusting it. In this work, the weight was set according to ODF peaks. The greater the ODF peak along fiber path, the greater the weight.
where is the weight, is the number of the consecutive directions, is the diffusion probability along the jth consecutive directions, and is the diffusion probability along the ith consecutive direction.
(3) The knot vector was normalized, and the nodes are distributed evenly. According to the number of control points of the path, the number of nodes of knot vector changes dynamically. The knot vector could be written as
where is the number of control points.
(4) NURBS fitting. In this procedure, we trace pathways which do not necessarily satisfy the control points precisely, but only approximately.
where are the third-degree B-spline basis functions defined on the knot vector . They could be computed using the Cox-deBoor algorithm [38]. To obtain , we have to compute the basis function first. There are at most four nonzero three-degree B-spline functions at each knot vector interval. So, we could directly get according to , , , and .
2.6. NURBS-G

In NURBS-G tracking, we do not consider the tangent relationship between fiber pathway and diffusion direction. The control points consist of only intersection points between the diffusion directions and the facets of the voxel. The 2D tracking process is demonstrated in Figure 4. The algorithm is outlined in Algorithm 3.

Input: consecutive diffusion directions
Output: fiber pathways
(1) Determine the control points. The set of control points is only composed of intersection points between diffusion direction and voxel border. The intersection points are acquired according to the equation given in the Algorithm 2.
(2) Compute the weights according to the equation given in the Algorithm 2.
(3) The knot vector is the same as the equation given in the Algorithm 2.
(4) NURBS fitting according to the equation given in the Algorithm 2.

3. Results

Figure 5 shows the ODF and vector fields estimated from HARDI images of tractometer. Panel (a) is the mask of fiber pathways. We extracted the diffusion directions corresponding to ODF local maxima that are above the mean value of ODFs. Through this filtration, spurious peaks could be effectively reduced [28].

After the vector fields were obtained, the control points and weights were computed. Next, the fiber pathways were traced with multidirectional streamline, NURBS-T, and NURBS-G. In this phantom experiment, is set to 60° and is 70 mm. was not set for this test, as WM mask was provided in tractometer dataset. Figure 6(a) shows 16 seed points selected according to [25], and 6(b) shows the ground truth fiber pathways. Figures 6(c), 6(d), and 6(e) show the tracking results.

In order to evaluate the proposed algorithms, two kinds of measure methods were taken. One is the point-to-point performance measures; the other is the connection measures. The former includes spatial metric (SM), tangent metric (TM), and curve metric (CM) [25]. These metrics focus on the point-to-point performance from a local perspective. The latter contains valid connections (VC), invalid connections (IC), no connections (NC), valid bundles (VB), and invalid bundles (IB) [39]. From a global point of view, the connections generated by the estimated trajectories are relevant. The set of global metrics takes into account the resulting connectivity. In this experiment, we evaluated the results with both local and global metrics. Figures 79 show the summation of the points per metric for each method. Table 1 shows the evaluation by using the global metrics: VC, IC, NC, VB, and IB.

We can come to that for the spatial metric NURBS-T obtains the best score except Fiber 3 and 10. For the tangent metric, NURBS-T also gets the best position except Fiber 10. For the curve metric, NURBS-T obtains the best place except for Fiber 9 and 15. Summarizing the overall performance over the three metrics, we can conclude that NURBS-T is best on the fiber pathway estimation of the phantom. For the computation time, NRBS-T recovered the previous results in about 23 minutes, and NURBS-G took about 20 minutes. The method of multidirectional streamline required 27 minutes or so to complete the task at the step of 0.02 mm. These methods were all implemented in Matlab R2014b running on the computer possessing 8G RAM and Intel Core i5-7200U.

From the above analysis, NURBS-T presents competitive results for both kinds of measure metrics. Furthermore, we used the mask (Figure 5(a)) to evaluate the resulting connectivity. The values in Table 1 show that the method with the best performance is NURBS-T.

Figures 1012 show the estimated fibers of the in vivo human brain data. In this in vivo experiment, is 60° and is 70 mm. is 0.15. We selected three ROIs to trace fiber pathways. The ROI in Figure 10(a) is located in the region of corpus callosum. The ROI in Figure 11(a) lies in the region of parietal lobe. The ROI in Figure 12(a) is in the region of bilateral mesial temporal lobes. As there is no golden standard of fiber distribution map with high resolution, we can only qualitatively analyze the results.

From Figure 10(b), we can easily pick out two fake fiber bundles that are marked by brown arrows. The thin bundle pointed by the left arrow is obviously nonexistent in the region of corpus callosum. The pathway pointed by the right arrow is unreasonable since it should not spread along the vertical direction. In Figure 10(d), from the morphological perspective, the fiber bundles are excessively messy and fluffy in the regions pointed by the two arrows because there are fewer constraints on the NURBS-G fitting. In Figures 11(b) and 11(d), there are too many crossing bundles, which disorderly emerge into the edge of WM in the region marked by arrows. In Figure 12(b), some unreasonable bundles could be found as their pathways spread out WM region. From Figure 12(d), we could see there are some minor bundles winds around the main bundles in the region pointed by the up-down arrow. In addition, the existence of the bundles in the regions pointed by the other three arrows is unreasonable.

From these in vivo tracking results, we can qualitatively validate our method. At last, to quantitatively analyze the proposed methods, we compared the results in the aspects of number of bundles, computation time, and storage (Table 2). The fiber bundles were stored as .mat file in Matlab 2014b. These methods were evaluated on the computer possessing 8G RAM and Intel Core i5-7200U CPU.

4. Discussion

In the presented study, we developed a novel tracking method based on NURBS curve fitting. The method consists of two steps. The first is to obtain the consecutive diffusion directions along a fiber pathway. The second is to carry out NURBS curve fitting. For the first step, we proposed a more effective way to find the consecutive vectors for a seed voxel among its 26-connected voxels. The comparison to FACT is shown in Figure 1. In the second step, the control points were obtained according to the equation given in the Algorithm 2. The corresponding weights are computed according to the equation given in the Algorithm 2. From the experimental results, we can conclude that the proposed method is well suited for exploring WM pathways.

The proposed method aims to reveal the connectivity among brain function areas. It is important to realize that our method does depend heavily on the parameters of control points and weights. Although we presented here both the theoretical foundation and a number of practical examples that characterize performance and accuracy of our approach, the main limitation of our work is the lack of a system wide analysis of the two parameters that can influence the fitting results. In NURBS fitting, we would continue to study the mathematical relationship between the weights and ODF peaks.

In general, there are two main factors influencing the tracking results: the noise in HARDI images and partial volume effects [40]. The noise could cause the inconsistency, and the incomplete information about partial volume effect could confuse the tacking process. In consequence, some fiber paths are incorrectly estimated [6]. Before the construction of ODF fields, we used NLPCA to denoise HARDI dataset. In the regions of fiber crossing, branching, and merging, the multiple compartments within a voxel make it hard to find out the fiber orientation from ODF fields for such entangled structures. In fact, the sensitivity to detect multiple fiber populations depends not only on the datasets but also on specifics of the construction technique of ODF. If the resolution capability of the construction method is low, the deviation between ODF maxima and the ground truth directions would become large. This error can limit the fiber tracking technique to fully delineate a fiber tract.

Another important factor that can influence the tracking results is stop criteria. FA could not be considered as one of the tracking stop criteria because FA is generally less than 0.2 in a voxel with crossing fibers [40]. Except for that, we considered the fiber length and the angle as stop criteria. However, validation of fiber tractography remains an open question [25].

5. Conclusion

Anatomical connectivity network is important to the investigation of human brain functions. The quality of anatomical connectivity relies on proper tract estimation [6]. In this work, we presented a novel algorithm based on NURBS curve fitting. The proposed methods exhibit promising potential in exploring the structural connectivity of human brain. They are easily implemented and proved efficient through phantom and real experiments. However, it is still difficult to identify the fiber bundles that are diverging, converging, and kissing. In future, our study will be mainly focused on how to solve this problem with NURBS fitting. More anatomical constraints should be used to guide tracking processes.

Data Availability

The tractometer and real datasets used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.


This study was supported by the Natural Science Foundation of Zhejiang Province (project no. LY17E070007) and National Natural Science Foundation of China (project no. 51207038).