White Matter Fiber Tractography Using Nonuniform Rational B-Splines Curve Fitting

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.


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 [2][3][4]. 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 [5][6][7]. At present, there are numerous tracking methods based on HARDI, which could be classified into deterministic and probabilistic algorithms [8].
ey 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. e 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.
ere 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. e 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].
e process starts with the definition of a point in WM region in which a virtual particle is allocated. e particle is iteratively moved along the local diffusion orientations until a stopping criterion is met. e 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].
is 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]. is 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. is data-driven approach was reported that it could improve valid neural connection rate over streamline methods. e other classes are probabilistic approaches. is 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. e 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]. e 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. is paper has comprehensively explored the tracking method based on NURBS curve fitting and has detailed how to determine the related parameters. e 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.
is method was evaluated on tractometer phantom and real brain datasets.

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 DOTmethod [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/mm 2 ) image. SNR of HARDI at b-values � 2000 s/mm 2 were evaluated. e 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. e acquisition parameters were as follows: two images with b � 0 s/mm 2 , 64 DW images with unique, and isotropically distributed orientations (b � 2000 s/mm 2 ). TR � 6700 ms, TE � 85 ms, and voxel dimensions equal to 2 × 2 × 2 mm. e SNR is, approximately, equal to 36.

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.
ere 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]. e 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. is 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.
(1) e term v i,x v i,y v i,z denotes a diffusion direction, and d i 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.

Diffusion Directions along a Fiber Pathway.
Before we conduct NURBS tracking, the consecutive directions along the same pathway have to be extracted. e 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. e 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)). e algorithm is summarized as Algorithm 1. e input parameters, including fiber length threshold L th , angle threshold θ th , and fractional anisotropy (FA) threshold FA th should be determined according to actual situation.

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]. e 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 [34][35][36]. e 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. e consecutive directions were used to compute control points. e weights were computed according to d i . 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). e whole procedure of NURBS tracking is shown in Figure 2.

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. e algorithm is outlined in Algorithm 2.

NURBS-G.
In NURBS-G tracking, we do not consider the tangent relationship between fiber pathway and diffusion direction.
e control points consist of only intersection points between the diffusion directions and the facets of the voxel. e 2D tracking process is demonstrated in Figure 4. e algorithm is outlined in Algorithm 3. 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. rough 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, θ th is set to 60°and L th is 70 mm. FA th was not set for this test, as WM mask was provided in tractometer dataset. Figure 6(a)   Input: ROI, L th , θ th , FA th . Output: Consecutive directions along a pathway.

Results
(1) Select one seed voxel in ROI.
(2) Choose one direction v � v i,x v i,y v i,z in the seed voxel, and its linear equation could be written as follows. cos θ � v · v′/(|v||v′|) (6) If |cos(θ)| ≥ |cos(θ th )|, v′ 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).  Journal of Healthcare Engineering Input: consecutive diffusion directions Output: fiber pathways C(u) (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. e intersection points could be obtained according to the equation given below: where a, b, and c are the length, width, and height of the voxel, respectively. (x, y, z) is the coordinate of the intersection point. Finally, we would get a series of points.
P � [p 1 , p 2 , p 3 . . . , p i , p i+1 , p i+2 . . . , p n−2 , p n−1 , p n ] where p i and p i+2 are the intersection points between the diffusion direction and the facets of the path-through voxel. p i+1 is the center point of the voxel.
(2) Calculate the weights. e 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. e greater the ODF peak along fiber path, the greater the weight.
where w is the weight, m is the number of the consecutive directions, d j is the diffusion probability along the jth consecutive directions, and d i is the diffusion probability along the ith consecutive direction. where n 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. 3 (u) are the third-degree B-spline basis functions defined on the knot vector u. ey could be computed using the Cox-deBoor algorithm [38]. To obtain C(u), we have to compute the basis function N i,3 (u) first. ere are at most four nonzero threedegree B-spline functions at each knot vector interval. So, we could directly get C(u) according to N i−3,3 (u), N i−2,3 (u), N i−1, 3 (u), and N i,3 (u).
ALGORITHM 2: Summary of NURBS-T fiber tracking.

Seed voxel Diffusion directions
Voxel center point Intersection point 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. e former includes spatial metric (SM), tangent metric (TM), and curve metric (CM) [25]. ese metrics focus on the point-to-point performance from a local perspective. e 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. e set of global metrics takes into account the resulting connectivity. In this experiment, we evaluated the results with both local and global metrics. Figures 7-9 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. e method of multidirectional streamline required 27 minutes or so to complete the task at the step of 0.02 mm. ese 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. e values in Table 1 show that the method with the best performance is NURBS-T. Figures 10-12 show the estimated fibers of the in vivo human brain data. In this in vivo experiment, θ th is 60°and L th is 70 mm. FA th is 0.15. We selected three ROIs to trace fiber pathways. e ROI in Figure 10(a) is located in the region of corpus callosum. e ROI in Figure 11(a) lies in the region of parietal lobe. e 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.       Journal of Healthcare Engineering From Figure 10(b), we can easily pick out two fake fiber bundles that are marked by brown arrows. e thin bundle pointed by the left arrow is obviously nonexistent in the region of corpus callosum. e 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). e fiber bundles were stored as .mat file in Matlab 2014b. ese methods were evaluated on the computer possessing 8G RAM and Intel Core i5-7200U CPU.

Discussion
In the presented study, we developed a novel tracking method based on NURBS curve fitting. e method consists of two steps. e first is to obtain the consecutive diffusion directions along a fiber pathway. e 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. e 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. e 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. e 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]. e 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. is 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].

Conclusion
Anatomical connectivity network is important to the investigation of human brain functions. e quality of anatomical connectivity relies on proper tract estimation [6]. In this work, we presented a novel algorithm based on NURBS e proposed methods exhibit promising potential in exploring the structural connectivity of human brain. ey 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 e tractometer and real datasets used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.