A PSO-Powell Hybrid Method to Extract Fiber Orientations from ODF

High angular resolution diffusion imaging (HARDI) has opened up new perspectives for the delineation of crossing and branching fiber pathways by orientation distribution function (ODF). The fiber orientations contained in an imaging voxel are the key factor in tractography. To extract real fiber orientations from ODF, a hybrid method is proposed for computing the principal directions of ODF by combining the variation of Particle Swarm Optimization (PSO) algorithm with the modified Powell algorithm. This method is comprised of the global searching ability of PSO and the powerful local optimizing of Powell search. This combination can guarantee finding all the diffusion directions without applying sliding windows and improve the accuracy and efficiency. The proposed approach was evaluated on simulated crossing-fiber datasets, Tractometer, and in vivo datasets. The results show that this method could correctly identify fiber directions under a range of noise levels. This method was compared with the state-of-the-art methods, such as modified Powell, ball-stick model, and diffusion decomposition, showing that it outperformed them. As to the multimodal voxels where different fiber populations exist, the proposed approach allows us to improve the estimation accuracy of fiber orientations from ODF. It can play a significant role in the nerve fiber tracking.


Introduction
At present, tractography based on diffusion-weighted magnetic resonance imaging (DWI) is the only noninvasive tool to obtain information on the neural architecture of the human brain white matter (WM) in vivo. The structural connectivity inferred from tractography is critical for understanding the functional coupling between cortical regions of the brain and for the characterization of neurodegenerative diseases and for medical applications [1][2][3][4]. In deterministic tractography, it is the important step to resolve the fiber orientations populated in each imaging voxel.
In current literature, there are three mathematical models applied to retrieve fiber orientations from DWI raw datasets: apparent diffusion coefficient (ADC), diffusion tensor (DT), and ODF. However, the local maxima of ADC profile do not necessarily coincide with the underlying fiber directions, making the extraction of orientation information difficult [5][6][7]. This is due to the nature of the ADC measurement which is the projection of spin displacements onto the diffusing gradient axis. The limitation of the DT model is the Gaussian diffusion assumption, which implies that there can only be a single-fiber population per voxel [8][9][10]. It is known that many voxels have low diffusion anisotropy due to the crossing, branching, and fanning of multiple fibers [11][12][13][14]. The ODF is defined as the radial projection of the spherical diffusion function, which is a function on the unit sphere describing the probability averaged over the voxel that a particle will diffuse into any solid angle [15,16]. As a spherical function, ODF has its local maxima aligned with the underlying fiber directions at every voxel. Until now, ODF is most widely employed to determine the fiber orientations with high angular resolution.
As the water molecules in WM tissues tend to diffuse along fibers when contained in fiber bundles, the principal directions of ODF agree with the true synthetic fiber directions [17]. The ODF field is promising for the estimation of neuronal fibers. The orientation of a particular fiber population could be estimated by finding the peaks of the corresponding reconstructed ODF. For this reason, a major focus within the DWI community has been directed at developing methods to compute the ODF. Although diffusion spectrum imaging (DSI) was firstly developed to image complex distributions of intravoxel fiber orientations with a more detailed, complete, and accurate view of WM architecture locally, the time-consuming MRI signal sampling restricts its application [18]. Currently, Q-ball imaging (QBI), constant solid angle QBI (CSA-QBI), and diffusion orientation transform (DOT) are constantly used to construct the ODF because they are more time-saving than DSI [15,16,19]. However, there is still an important issue to be solved in the neural fiber tracking from ODF. That is how to extract the fiber directions from ODF accurately with high angular resolution. In general, there are two kinds of methods to extract fiber orientations from ODF, and the first is to resolve the ODF to multicompartment partial volume models, and the other is to directly search the peaks of ODF.
The method of diffusion decomposition obtains fiber orientations by decomposing ODF into standard component ODFs [20]. The advantage of the approach is that a fiber vector can be easily represented by the direction of component ODF but not by spherical harmonics. The decomposition algorithm provides a sparse solution to improve the ability in resolving crossing fibers and to avoid false fibers as encountered in diffusion deconvolution. Ball and sticks model is a model-fitting approach, which decomposes HARDI signals into isotropic and anisotropic diffusion components directly [21]. However, the two models suffer from the shortcomings regarding model selection. We must determine the number of diffusion compartments with a priori structural knowledge of each voxel and must use nonlinear optimization to obtain the fiber orientations. What is more, the two methods are sensitive to noise and to the number of HARDI measurements. Essentially, the two models are ill-posed inverse problems.
The ODF, as a probability distribution function, should be nonnegative. As the principal directions of ODF are consistent with fiber orientations, extracting the fiber orientations from ODF could be boiled down to a multimodal optimization problem. At present, several methods exist to extract ODF's maxima, such as finite difference method, Powell's method, and spherical Newton's method. By searching for local maxima of persistent angular structure (PAS) using Powell's method, the orientations of WM fibers were revealed [22]. Tournier et al. estimated the orientation of a particular fiber population by finding the peak of the corresponding reconstructed ODF using a spherical Newton's method [23]. This method has the merit of high convergence speed, but it is susceptible to the position of the starting point. In the iteration not only gradient vectors and their modulus but also Hessian matrix and its inverse matrix are needed to be computed over and over again. This is quite time-consuming and memory-consuming. Sequential quadratic programming (SQP) is an iterative method for nonlinear optimization. It usually is used on the mathematical problems for which the objective function and the constraints are twice continuously differentiable [24]. But SQP was found to introduce biases in the peak distributions via the constraints.
In order to find all the fiber-along vectors from ODFs and improve the precision of multipeak searching, in this work, we have introduced a novel methodology to estimate the fiber directions directly from ODF based on PSO-Powell hybrid algorithm. In this method, the global search ability of PSO is combined with the strong local search ability of modified Powell algorithm. This combination can not only improve the solution accuracy but also speed the searching at the same time. Only using the function value information without the need to calculate derivatives makes it very useful to solve ODF optimization. It can correctly retrieve the orientations corresponding to underlying intravoxel fibers populations. Results on the simulated datasets, Tractometer, and in vivo HARDI datasets illustrate the effectiveness of the proposed approach.

ODF Construction.
In the past decades, respectable researchers have tried to extract fiber orientations with high angular resolution from raw DWI datasets. However, image acquisition under clinical conditions with limited measurement time faces the problem of poor spatial and angular resolution and the technique's high susceptibility to noise [25][26][27]. The advent of HARDI has provided the chance to delineate multifiber pathways effectively and efficiently by ODF. Essentially, ODF is a function of two angular variables and as (1). The equation expresses diffusion probability in the direction ( , ) of water molecules contained in WM tissue: where is the displacement radius, is the azimuth angle, and is the elevation angle in spherical coordinate. In this work, we applied the PSO-Powell hybrid method to the ODF fields which were constructed with QBI, CSA-QBI, and DOT. In order to be able to delineate fiber crossings even at low angles and avoid unnecessary loss of angular resolution, we choose a high SH order of 8 for ODF constructions in CSA-QBI and DOT. Higher-order spherical harmonics are necessary to resolve fibers that are separated by small angles but also introduce noise [28][29][30]. We applied a set of 724 directions evenly distributed on a unit sphere to evaluate the ODFs of the testing HARDI datasets.
After ODF fields have been constructed, it is crucial to resolve the principal directions of ODF, which are aligned with the underlying fiber directions. In each imaging voxel, the directions of fiber tracts are parallel to the directions of maximum diffusion that are defined as local maxima of the ODF [18,19,23,24,31,32]. In the next section, the whole process of PSO-Powell hybrid algorithm was described in detail, which was applied to extract the fiber directions through finding the peaks of ODF.

PSO-Powell Hybrid
Optimization. Hybrid strategies for optimization are implemented by combining a heuristic algorithm with a mathematical algorithm. This strategy increases reliability in comparison to mathematical methods Computational and Mathematical Methods in Medicine 3 and increases efficiency in comparison to pure heuristics algorithms [33]. In this work, PSO is coupled with the modified Powell's method in order to obtain a fast and reliable hybrid algorithm. As the heuristic algorithm, PSO is utilized to cover the entire search space, while modified Powell method, as the mathematical algorithm, starting from a point inside this region, quickly reaches local maxima. The combination makes the hybrid algorithm reliable and at the same time maintains properties which lead to rapid convergence. The PSO algorithm is used to cover the entire search space, identifying the region of local maximum, while Powell algorithm quickly reaches the maximum. Another reason for this choice is the fact that this hybrid algorithm does not require the gradients of ODF [34,35].
In the optimization of ODF, in order to make PSO go through the domains where the local maxima locate, the c 2 parameter in (2) is assigned to zero [34]. Since the two random factors of rand( ) in (2) could increase the stochastic motion of the particles, we removed them from (2) so that the hybrid algorithm could converge to all local maxima as soon as possible. The revised particles' velocity update equation is described by (3).
where ( ) is the inertia weight, stands for the velocity of particles, 1 and 2 represent learning factors, rand( ) represents a random value between 0 and 1, is the optimal position of the th iteration, is the global optimal position, and stands for the current position. Because the construction points located on the spherical surface have no structure or order between their relative locations, we interpolated the query points by triangulating the known ones. This step involves traversing of the triangulation data structure to find the triangle that encloses the query point. Once the point is found, the subsequent step is to compute the value of the query point by the nearest-neighbor interpolation method. The detailed procedure of PSO-Powell hybrid searching algorithm is outlined as follows.
Step 1. Initialize the positions for a swarm of particles of size , and initialize the parameters including the maximum number max of PSO iterations, the maximum number of hybrid iterations, the precision of Powell searching, and the precision of PSO searching.
Step 2. Evaluate the fitness of each particle.
Step 3. If the total number of iterations is greater than , we would stop the iteration and output all the local maxima. Otherwise, turn to Step 4.
Step 4. If ≤ max ( is the number of PSO iterations), the speed and position of the particles would be updated according to (3) and (4). And then and should be updated at the same time.
Step 5. If Tolerance PSO < , then we search the extrema for the particles of last PSO iteration using modified Powell method.
where and denote the step length. And is obtained by the linear search.
Step 6. The new extremum is added to the extreme set.
Step 7. Reinitialize particle position and velocity. Then go to Step 2.
In this hybrid method, the term of Tolerance PSO < is the condition of the transformation from PSO to Powell search, in which Tolerance PSO represents the tolerance of PSO optimization and stands for the Powell searching threshold value. In the construction of ODF, some noise would be introduced. We could select a normalized ODF value as the thresholding, and this could avoid selecting small peaks that may appear due to noise and transformation [27].   This hybrid algorithm could search all the peaks in the feasible region through revised PSO and has the ability of global search for all extremum points. With the strong local search ability of modified Powell method, the accuracy and convergence speed of the hybrid algorithm are improved. The algorithm only uses the value of the ODF without the need to calculate derivatives. In this work, the parameters of the hybrid algorithm are set as follows according to [34,35]. PSO population size is 100. The max inertia weight is 0.2, and the min inertia weight is 0.1. The acceleration factor 1 is 0.5. The maximum number of PSO iterations is 120. The searching threshold value is 0.02. The maximum number of PSO-Powell hybrid iterations is 20.

Results
We utilized multitensor simulated datasets, Tractometer datasets, and in vivo datasets to evaluate the methods for extracting fiber orientations, including ball-stick, modified Powell, diffusion decomposition, and PSO-Powell model. The angular deviations for the four methods were compared and the results proved the validity and feasibility of PSO-Powell hybrid algorithm.

Simulation Study.
The synthetic datasets were acquired using the multitensor model [36], which leads to an analytical computation of exact ODF. For a given -value of 3000 s/mm 2 , noise level of SNR = 20 [37], and 64 encoding directions that uniformly are distributed on the unit sphere, we generated DWI raw signals. The simulation parameters of the synthetic datasets about one single-fiber and two and three crossing fibers are shown in Table 1. After the ODF fields were constructed with QBI, CSA-QBI, and DOT, we applied three algorithms including PSO-Powell, modified Powell with sliding windows, and diffusion decomposition to extract the fiber orientations. The constructed ODFs were displayed in Figure 1. The ball-stick model is a simplified model for multiple tensors, and the fiber orientations are directly estimated from DWI raw signals. We directly applied it to the comparison in Figures 2-4. Figures 2, 3, and 4 show the angular deviations of ballstick, PSO-Powell, modified Powell, and diffusion decomposition applied to synthetic ODFs without noise and with noise of SNR = 20. The legends located on the right side denote the angular deviations of the evaluated methods. Ballstick model directly resolves the fiber directions from raw DWI signals. The other three methods were used to extract fiber directions from ODFs constructed with QBI, CSA-QBI, and DOT. The azimuth and elevation deviations were separately grouped for quantitative comparison. As shown in the figures, the diffusion decomposition with QBI ODF has the highest angular deviation for the estimated azimuth and elevation . The results of the modified Powell method show better accuracy for noised and noiseless ODFs. Its angular deviation may be due to ODF construction incompatibility and sliding windows radius, and the error can be reduced in low SNR conditions. The ball-stick model performs best for single-fiber and two-crossing-fiber ODFs (Figures 2 and 3) but worse for three-crossing-fiber ODF (Figure 4). In sum, PSO-Powell method showed substantial better performance than other methods. The overall results suggest that PSO-Powell method can be applied to the fiber orientations extraction from ODF fields constructed with QBI, CSA-QBI, and DOT. From the above quantitative comparisons, the hybrid algorithm was shown to be reliable and to perform better than the ball-stick model, modified Powell  draw the conclusion that PSO-Powell got the better convergence speed than modified Powell and diffusion decomposition. Because ball-stick model resolves the orientations through nonlinear curve fitting in least-squares sense under the condition that the number of diffusion compartments is given, it requires significantly less time. Powell's method takes the longest time because of the repeated application of the sliding window, and the radius is 0.4. In diffusion decomposition, much time must be spent on the iterative regression analysis. These methods were tested on the PC equipped with Intel Core i5-3337U, 4 G RAM, and Windows 10.

Phantom Study.
The proposed method was also evaluated on the phantom of Tractometer, which is popular in fiber tracking test. The DWI datasets of Tractometer were acquired on the 3T Tim Trio MRI systems, equipped with a whole body gradient coil, a whole body transmit coil, and a 12channel head receive coil. A single-shot diffusion-weighted twice refocused spin echo-planar pulse sequence was used to perform the acquisitions in order to compensate for the first order Eddy currents. For each acquisition, there were 64 uniformly distributed diffusion-weighted measurements and one = 0 image, with two repetitions. The spatial resolution for the datasets is 3mm × 3mm × 3mm and 3 slices were   [38]. Figure 5 shows visual comparison of the fiber orientations extracted from ODF fields constructed with QBI, CSA-QBI, and DOT. The ROI regions were marked out by the red squares in Figure 5(a). The top-right ROI contains twocrossing fibers, and the bottom-left ROI contains single fibers. From Figure 5(c), we could see that the marked voxels lost one fiber orientation. A possible cause to this error can be due to the discrepancy between QBI and the actual diffusion pattern in the phantom. In Figure 5(e), the voxels marked by black polygons got some more false directions. In Figure 5(g), the four marked voxels have false orientations, and this is due to the fact that DOT would lead to lager false peaks. In Figure 5(i), the single-fiber directions were correctly estimated. In Figure 5(k), seven voxels have false directions. In Figure 5(q), three voxels have false directions. Overall, in crossing-fiber regions, PSO-Powell method has best performance on the ODF field constructed with DOT. In single-fiber regions, the method stands out on the ODF field constructed with QBI.

In Vivo Study.
In order to further illustrate the effectiveness of PSO-Powell hybrid algorithm on neural fiber orientation extraction, we applied the algorithm to in vivo DWI datasets. The whole-brain HARDI was performed on an adult male volunteers (25 years old) using a 3T Signa EXCITE scanner equipped with an eight-channel phased-array head 8 Computational and Mathematical Methods in Medicine coil. A multislice single-shot echo-planar spin echo pulse sequence was employed to obtain attenuated signals at a diffusion weighting of = 3000 s/mm 2 , where the diffusionencoding directions were distributed uniformly over the surface of a sphere using electrostatic repulsion. An additional acquisition without diffusion weighting at = 0 s/mm 2 was also obtained. The total scan time for whole-brain acquisition of 64 diffusion-encoding directions was 19.3 min with TR/TE 18 s/84 ms and isotropic 2 mm voxel resolution (FOV 260 × 260 mm, matrix 140 × 140, 92 interleaved slices with 2mm slice thickness and no gap). Figure 6 visually shows the fiber orientations estimated from ODF fields constructed with QBI, CSA-QBI, and DOT in the region of the corpus callosum (CC) that is a structure that connects the left and right cerebral hemispheres. After the ODFs were computed, the condition of FA ≥ 0.20 was applied to segment WM tissue out to calculate fiber directions. The results of Figures 6(c) and 6(e) are consistent with known anatomy. In Figure 6(g), some voxels have twocrossing directions, and this is mainly because that DOT has evoked false peaks into ODF due to the impact of MRI Rician noise. Figure 7 shows the fiber orientations estimated from ODF fields constructed with QBI, CSA-QBI, and DOT in the ROI region through which cuneus fibers, vertical occipital fasciculi, and superior longitudinal fasciculi pass. This region contains many multimodal ODFs. In Figure 7(g), crossingfiber directions can all be clearly identified. By contrast, Figures 7(c) and 7(e) only reflect one fiber direction mainly. In the calculation, the directions along which the peak of ODF was less than the mean value were discarded [23], just like the ODFs marked out by red rectangles in Figures 7(d) and 7(f).

Discussion
The Powell-PSO hybrid algorithm overcomes the shortcomings of PSO algorithm, such as poor accuracy and slow convergence speed. At the same time, it keeps the merits of Powell method, such as computational accuracy and convergence rate. By searching all the peaks in the feasible region, the hybrid algorithm has the ability of searching all the local maximums of ODFs. Since the algorithm only uses the value of the function information without the need to calculate derivatives, it is suitable for solving the problem of differentiable and nondifferentiable multifield optimization (MFO) [35].
The comparisons with ball-stick model, modified Powell, and diffusion decomposition methods show that this hybrid algorithm is effective in searching the diffusion directions from ODFs. But when ODFs have more spurious local maximums due to noise or computational models, more hybrid iterations would be carried out.
In the proposed method, the c 2 parameter is set to zero to increase the stochastic motion of PSO particles. This would make the hybrid algorithm go through the domains where local maxima locate. The other parameters were determined according to [34,35], such as inertia weight , acceleration factor 1 , searching threshold , and the number of hybrid and PSO iterations. In simulated experiments, the proposed hybrid method has been tested on single-fiber, two-fiber-crossing, and three-fiber-crossing ODFs. Through quantitative comparisons (Figures 2-4, Tables 2-4), it outperformed other methods. In Tractometer and real experiments, PSO-Powell method has also been verified using single-fiber and crossing-fiber ODFs constructed with QBI, CSA-QBI, and DOT. The results (Figures 5-7) conformed with the ground truth satisfactorily.
The diameter of bundles of axons considered in fiber tractography are on the order of 1 mm, and individual physical fibers on the order of 1-30 m. At the current resolution of DWI, there are between one-third and two-thirds of imaging voxels in the human brain WM that contain fiber crossing bundles [1,8]. But higher resolution of fiber directions comes at the cost of higher susceptibility to noise. The low detection rate and accuracy of local maximums of ODFs at critical crossing angles are well-known problems.
Theoretically, there are two types of parameters that can be extracted from ODFs: the orientation of each fiber population and its volume fractions. Having estimated the orientations of the various fiber populations, the corresponding volume fractions could be computed by finding the set of weighting factors providing the best fit of the weighted sum of their respective signal attenuation profiles to the actual measured signal attenuation. Further study should be conducted to propose a more mature model for computing the fiber volume fractions, meanwhile identifying the fiber orientations.
HARDI has become a common tool for the reconstruction of WM architecture by ODF. In recent years, numerous tracking methods based upon HARDI have emerged in order to overcome the shortcomings of ADC and DT which has no the ability to resolve complex WM architecture such as fiber crossing, branching, and kissing. Furthermore, HARDI is more time-saving than DSI because it just needs to sample on a spherical surface in -space. Due to the fact that higher sensitivity to crossing fibers results in higher sensitivity to noise, the false fiber orientations created by noise can cause substantial error in fiber tracking [39]. In order to increase the accuracy of local maximum detection of ODFs and avoid leaving out nondominant diffusion peaks, we should take into account information from neighboring voxels in further research. In order to better estimate the fiber orientations, we may present a straightforward yet effective method for the smoothing, regularization, and sharpening of diffusion profiles in crossing areas of ODF fields obtained from DWI datasets.

Conclusion
In this work, a novel hybrid method for determining the principal directions of the ODF is proposed by combining the PSO and the modified Powell search. The method tends to extract fiber orientations from ODF fields, leading to no less resolution, but with fewer and smaller spurious peaks, particularly at low SNR. The method can improve the estimation of fiber orientations in heterogeneous WM regions and boost the reliability of fiber tracking. On the basis of the quantitative analyses with the synthetic and phantom datasets, we can conclude that this method is promising to improve the estimation accuracy of fiber orientations from ODF. The experiments with data from a healthy human subject, acquired under clinical imaging conditions, show the method's potential to notably optimize the reconstruction of noncrossing and crossing-fiber orientations. This method may improve the tracking for the construction of brain structural connectivity network.

Conflicts of Interest
The authors declare that they have no conflicts of interest.