Computational and Mathematical Methods in Medicine

Volume 2018 (2018), Article ID 7680164, 12 pages

https://doi.org/10.1155/2018/7680164

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

^{1}School of Electronic Information, Hangzhou Dianzi University, Hangzhou, China^{2}Department of Systems Medicine & Bioengineering, Houston Methodist Hospital, Houston, TX, USA^{3}Department of Biomedical Engineering, University of Houston, Houston, TX, USA

Correspondence should be addressed to Xiaohui Yu; gro.tsidohtemnotsuoh@2uyx

Received 4 August 2017; Revised 20 December 2017; Accepted 26 December 2017; Published 21 January 2018

Academic Editor: Chuangyin Dang

Copyright © 2018 Zhanxiong Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. 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–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–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–10]. It is known that many voxels have low diffusion anisotropy due to the crossing, branching, and fanning of multiple fibers [11–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.

#### 2. Methods

##### 2.1. 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–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–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.

##### 2.2. 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 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 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, and represent learning factors, represents a random value between 0 and 1, is the optimal position of the 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 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 ( 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 , then we search the extrema for the particles of last PSO iteration using modified Powell method.(5.1)The particles of the last PSO iteration are considered as the initial points, . The directions are linearly independent. Generally, are set along the directions of axes. Let .(5.2)Starting from , are obtained by the linear search along the directions of : where and denote the step length. And is obtained by the linear search.(5.3)Let . If , is the solution. Otherwise, starting from , a new solution would be found along by the linear search.(5.4)The parameter corresponding to the maximum drop is determined through(5.5)If would be still linearly independent, and they are still the search directions of the next iteration. Let , then go to (5.2).(5.6)If (7) is false, are linearly dependent. Let (, and ; then go to Step (5.2).

*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 is the condition of the transformation from PSO to Powell search, in which 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 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.

#### 3. 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.

##### 3.1. 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 [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.