Keywords

1 Introduction

There are three approaches to analyze the fiber bundles, namely, ROI-based analysis (RBA) [1], voxel-based morphometry (VBM) [2, 3], and tract-based analysis (TBA) [4,5,6]. The ROI-based analysis identifies the differences about the region of interest based on the brain map. Voxel-based morphometry searches local changes in the grey matter density according to the T1-weighted structural MRI brain images. Although the two aforementioned approaches have been widely used, there is insufficient information about fiber bundles contained, thus making it difficult to correlate fibers directly with the specific brain function [7]. To detect more fiber features, quantitative analysis of DTI data along the white matter tracts is presented and is known as tract-based analysis. Furthermore, tract-based spatial statistics (TBSS) [8] is created combining VBM and TBA for estimating a group mean FA skeleton which represents the centers of all fiber bundles. It is more sensitive to the FA value than VBM, but the FA skeleton is easily break off at the crossing fiber bundles which generates an inexact result. Spatial registration and analysis at the fiber level is a promising research diretion to overcome the disadvantage of TBSS.

Registration on fiber bundles is essential to compare and analyze fibers from different subjects directly. The reliability of the quantitative analysis relies on the accuracy of registration. During the past three decades, a large amount of studies have been conducted on the registration of fibers. Consequently, three registration methods are evolved, namely, scalar registration, tensor-based registration, and fiber-based registration. Based on the previous studies [9], the transformation matrix of registration is calculated based on one or more channels of the scalar DTI images, such as FA and B0. There are certain limitations of the scalar or vector registration for lacking the information of orientation, only based on the intensity message. Tensor-based registration methods take the directional signal into account in contrast to the scalar registration as mentioned [10,11,12]. However, tensor reorientation may introduce extra errors for tractography, thus restricting the application of tensor-based methods. Direct fiber registration is also proposed in [13, 14], where each fiber is projected into a high dimensional feature space leading to model and target feature points sets. The paper utilized the Coherent point drift (CPD) [15] algorithm to perform a non-rigid registration in order to make a point-to-point correspondence along the fiber tract. According to the CPD algorithm, every point on the fiber bundles is treated as a normal distribution while the optimal matching is completed by the optimization procedure [16].

2 Methods

2.1 Overview and Preprocessing

The main steps of the proposed method for registration of fibers are presented in Fig. 1. Preprocessing includes eddy current correction, tensor calculation, channel image generation and tractography by Diffusion toolkit (DTK) [17]. In tractography, the FA threshold is set as 0.15, and the turn angle threshold is set as 35°. After preprocessing, automated fiber tract clustering is applied to obtain a set of reliable and compact fiber tracts. Meanwhile, the central fiber for each fiber tract is extracted to represent the whole bundle. Then, the center fibers are registered to establish the correspondence among bundles for tract-based analysis. Each step is described in detail in the following subsections.

Fig. 1.
figure 1

Main steps of the proposed method.

Full size image

2.2 Clustering and Central Fiber Extracting

Here, using the white matter query language method described in [18], we extract the corpus callosum (CC) connecting left/right frontal lobe, the uncinate fasciculus (UNC), and the superior longitudinal fasciculus (SLF) for analysis. The model fibers and reference fibers are described as point sets in 3D coordinates.

After clustering, the center fiber for each fiber bundle is extracted. The center fiber fc of the fiber bundle F is defined as the fiber with the smallest average distance from all other fibers in F, where the distance is measured by Hausdorff distance. The Hausdorff distance refers to the maximum distance between one point set and the nearest point set of another. dm,n is the Hausdorff distance of fiber m and fiber n.

$$ d_{m,n} = max_{{p_{k} \subset F_{m} }} min_{{p_{l} \subset F_{n} }} \left\| {p_{k} - p_{l} } \right\| $$
(1)

where pk and pl are two point sets belong to fiber m and fiber n, respectively.

$$ f_{c} = Arg\,{ \hbox{min} }\,d\left( f \right) $$
(2)

Figure 2 shows the fiber tract clustering and central fiber extracting results for the SLF, UNC, and CC connecting left/right frontal lobe are showed, and the red line represents the central fiber for each bundle.

Fig. 2.
figure 2

Results of the central fiber extracting, red lines represent the central fibers. (Color figure online)

Full size image

2.3 Non-rigid Registration

CPD is a non-rigid matching algorithm based on a uniform velocity field, which uses the calculus of variations to register the maximum likelihood estimation after normalization [15]. The basic idea of the proposed registration method is to measure the correlation between two fiber sets by considering their continuous approximations. To solve the flexible deformation and relative relation of the fiber point set data, we use Gaussian mixture model (GMM) to fit the fiber point sets to get a probability distribution. Then the EM (Expectation Maximization) algorithm is used to optimize the parameters with the theory of the maximum likelihood estimation. The transformation is derived after several terms of optimal iteration. Specifically, CPD not only can achieve the registration of non-linear and non-rigid point sets, but also has an accurate result in case of large noise and frequent changes.

We assume the point set of the model central fiber Mξ*i = (M1, …, Mi)T as the kernel of GMM and the point set of the target central fiber data Tξ*j = (T1, …, Tj)T as the data set of GMM. The starting point of the model fiber is Y = v(M0) + M0.

$$ p\left( t \right) = \omega \frac{1}{j} + \left( {1 - \omega } \right)\sum\nolimits_{x = 1}^{i} {\frac{1}{i}p\left( {t |x} \right)} $$
(3)
$$ p\left( {t |x} \right) = \frac{1}{{\left( {2\pi \sigma^{2} } \right)^{{\frac{\upxi}{2}}} }}exp\left( { - \frac{{\left\| {t - M_{i} } \right\|^{2} }}{{2\sigma^{2} }}} \right) $$
(4)

The formula above is the probability density of GMM. \( \omega \in \left( {0,1} \right) \) is a weighting parameter of spill points, redundant points and noise points. ξ represents the dimension of the point set, which in this article is 1. The p(t) expresses the probability that t, the data point, is generated by M GMM centers (including the influence of noise). After setting the parameters of the model, EM algorithm is used to optimize the calculation to get the transforming relationship between two point sets.

The centroid position of the point set Tj is adjusted by a series of transformation parameters that can be estimated by minimizing the following negative log likelihood function below.

$$ E(\theta ,\upsigma2) = - \sum\nolimits_{n = 1}^{N} {log} \sum\nolimits_{m = 1}^{M} {p\left( x \right)p\left( {t |x} \right) } $$
(5)

EM Algorithm is used to calculate θ and σ2. E step establishes the target function Q which represents the upper bound of the negative log-likelihood function (5). Based on the Gaussian mixture model clustering and EM algorithm, E-step can be derived as foll-ows:

$$ Q\left( W \right) = \sum\nolimits_{m = 1}^{M} {\sum\nolimits_{n = 1}^{N} {p\left( {m |x_{n} } \right) \times \left\| {\frac{{x_{n} - y_{m} - G\left( {m,R} \right)W}}{\sigma }} \right\|^{2} /2^{{\sigma^{2} }} + D} } $$
(6)
$$ P_{mn} = exp\left( { - \frac{1}{2}\left\| {\frac{{y_{m} - x_{n} }}{\sigma }} \right\|^{2} } \right)/\sum\nolimits_{m = 1}^{M} {{ \exp }\left( { - \frac{1}{2}\left\| {\frac{{y_{m} - x_{n} }}{\sigma }} \right\|^{2} } \right)} $$
(7)

The posterior probability P can be obtained by calculating the parameter values, where G is the matrix, m and R are the rows and columns of the matrix. By calculating the partial derivative of W for the above equation, we can obtain Eq. (7) above.

3 Experiment and Results

3.1 Registration Based on CPD

All studies were obtained on 1.5T MR units (Siemens, Sonata, VA25 operating system). Twelve DT imaging data were acquired by using a single-shot, echo-planar imaging sequence with sensitivity encoding and a parallel imaging factor of 2.0. The imaging matrix was 96 × 96 with a field of view of 240 × 240 mm. Diffusion weighting was encoded along 30 independent orientations and the b-value was 700 s/mm2 [19]. Three kinds of fiber bundles are selected to make the registration based on CPD, namely, CC, UNC and SLF. The registration results are showed in Fig. 3. The blue line presents the model set while the red line shows the target. In order to verify the reliability and superiority of the proposed tract-based registration method, we chose three other classical registration methods for comparison, namely Affine [20], Rigid [21], and ICP [22].

Fig. 3.
figure 3

Results of fiber registration based on CPD

Full size image

3.2 Evaluation Based on Distance Between Fiber Tracts

We choose to implement the mean of the closest distances since it provides an estimate that uses all the available data. It is more discriminative than the minimum of the closest distances. We use K-means algorithm to find the nearest point of the target fiber set to the model set which is defined to be the mean nearest distance as the evaluation based on the distance between fiber tracts. Here, Mean Nearest Neighbor Distance (MND) is defined to be the index. The Euclidean distance is used to measure the distance of two samples where Ni represent the nearest point of the registered fibers from the model fiber point set Mi as Eq. (8) shows.

$$ MND = \frac{1}{n}\mathop \sum \limits_{i = 1}^{n} \left\| {N_{i} - M_{i} } \right\| $$
(8)

Figure 4 shows the Mean Nearest Neighbor Distance between model fibers and the reference fibers of twelve subjects. Smaller values of MND show a better registration because it indicates the difference levels after registration between each subject and the model. As we can see, CPD shows the lowest value while the other three methods exhibits a higher result, indicating the CPD registration algorithm stands out among the registration methods evaluated.

Fig. 4.
figure 4

The Mean Nearest Neighbor Distance of CC (a), SLF (b), UNC (c), and the mean of twelve subjects with the Standard deviation (d) on four registration methods

Full size image

3.3 Evaluation Based on FA Profile Along the Fiber Tracts

A key function of registration is for comparison among individual images. So the registration assessment performance also should be the measurement of anatomical structures. A profile metric based on the normative correlation is calculated along each fiber. The FA value is obtained on the corresponding points captured from template for each registered subject. Figure 5 show the FA profile of three selected fibers from 12 subjects using different tract-based registration methods. In these figures, the x-coordinate represents the arc length of the fiber bundles, and the y-coordinate is the value of FA. According to the FA distribution, the mean FA distribution of CPD algorithm has the highest degree of template fitting and almost coincides. However, the distribution of other algorithms is not as high as the FA distribution of template fiber bundles, which indicates that CPD algorithm guarantees the distribution of fiber bundles’ characteristic parameters to the greatest extent while ensuring the registration accuracy. Therefore, the application of CPD algorithm in fiber bundle registration enables more accurate results to be obtained in subsequent analysis.

Fig. 5.
figure 5

FA profiles of the CC (a), SLF (b), UNC (c) for the four registration methods.

Full size image

On the basis of the automatic clustering and marking of fiber bundles, it is possible to find the landmarks whose shape features are stable according to the shape information of the fiber bundles. The feature points are described by maximum curvature along the fiber bundles and the landmark is also defined to be the origin of the FA profile. pl is the landmark point where the cosine value is the maximum. Figure 6 is the schematic diagram of fiber landmark where the deeper the red is, the greater curvature it shows.

Fig. 6.
figure 6

The schematic diagram of fiber landmark.

Full size image
$$ p_{l} = Arg\,max\,cos\left( {p_{i} ,p_{j} } \right) $$
(9)

4 Conclusion

In this paper, we proposed a novel tract-based analysis method based on the CPD non-linear registration algorithm. CPD algorithm is applicable to multi-dimensional point set registration under rigid and non-rigid transformation, it has strong robustness to the influence of noise, out of line point and missing point. To evaluate the performance of the proposed method, the mean nearest neighbor distance and the FA profile were calculated and the comparison with three traditional registration methods was conducted. The results support the superiority of CPD over the existing methods. Thus, the tract-based analysis based on the new CPD non-linear registration poses an exciting research direction for brain mapping.