Tri-linear interpolation-based cerebral white matter fiber imaging.

Diffusion tensor imaging is a unique method to visualize white matter fibers three-dimensionally, non-invasively and in vivo, and therefore it is an important tool for observing and researching neural regeneration. Different diffusion tensor imaging-based fiber tracking methods have been already investigated, but making the computing faster, fiber tracking longer and smoother and the details shown clearer are needed to be improved for clinical applications. This study proposed a new fiber tracking strategy based on tri-linear interpolation. We selected a patient with acute infarction of the right basal ganglia and designed experiments based on either the tri-linear interpolation algorithm or tensorline algorithm. Fiber tracking in the same regions of interest (genu of the corpus callosum) was performed separately. The validity of the tri-linear interpolation algorithm was verified by quantitative analysis, and its feasibility in clinical diagnosis was confirmed by the contrast between tracking results and the disease condition of the patient as well as the actual brain anatomy. Statistical results showed that the maximum length and average length of the white matter fibers tracked by the tri-linear interpolation algorithm were significantly longer. The tracking images of the fibers indicated that this method can obtain smoother tracked fibers, more obvious orientation and clearer details. Tracking fiber abnormalities are in good agreement with the actual condition of patients, and tracking displayed fibers that passed though the corpus callosum, which was consistent with the anatomical structures of the brain. Therefore, the tri-linear interpolation algorithm can achieve a clear, anatomically correct and reliable tracking result.

Research Highlights

(1) Tri-linear interpolation algorithm for fiber tracking can reduce the noise and partial volume effects, thus obtaining more rapid calculations, more tracked fibers, longer and smoother tracked fibers, more obvious orientation and clearer details.

(2) Comparisons of the tri-linear interpolation algorithm and the tensorline algorithm help to define the theoretic application value. Furthermore, the adaptivity and availability of the tri-linear interpolation algorithm in tracking white matter in human brains have been verified through experimental testing, clinical applicability and comparisons with actual anatomical structures.

INTRODUCTION

Fiber tracking based on diffusion tensor imaging is the predominant way to monitor white matter fibers three dimensionally, non-invasively and in vivo in recent years[12]. Fiber tracking can observe the sparse jostling, disruptions, destructions and other anomalies of white matter fiber tracts. Called tractography[3], these three-dimensional color coded maps can display the complicated directions of white matter fiber tracts more intuitively, and therefore has wider applications and receiving more attention in clinical and scientific research[45678].

At present, fiber tracking algorithms are divided into deterministic tractography[91011] and probabilistic tractography[12131415]. In deterministic algorithms, Mori et al[16] proposed the Fiber Assignment by Continuous Tracking algorithm, which tracks fibers along the main diffusion direction of a voxel. Weinstein and colleagues[17] studied the Tensor Deflection algorithm, which uses the entire diffusion tensor to deflect the estimated fiber trajectory. Subsequently, Weinstein and colleagues[18] proposed the tensorline algorithm, which combines the advantages of the Fiber Assignment by Continuous Tracking and Tensor Deflection algorithms. Moreover, some algorithms have also investigated the reduction of tracking errors by applying energy minimization methods. Parker and colleagues[19] adopted the Fast Marching Method to track white matter fibers. Prados et al[20] applied Riemannian geometry and control theory to track fiber bundles by computing the geodesic distance between seed points and termination points. In probabilistic algorithms, Wu et al[21] proposed the Global Optimization Algorithm and obtained an optimal tracking result by making the fiber tracts pass through two selected regions. In addition, Cheng et al[22] presented Tractography Incorporating A Priori Anatomic Knowledge, which combines known anatomical structures and achieved a more accurate tracking result.

Deterministic algorithms track fibers mainly depending on diffusion direction, which are simple and can lighten the calculation burden, and have been most widely used in clinical diagnosis[23]. However, deterministic algorithms are susceptible to noise and partial volume effects, which result in the accumulation of tracking errors. Probabilistic algorithms can effectively reduce noise and partial volume effects, thus decreasing the accumulated errors and providing more fiber orientations. Unfortunately, their calculations are very complicated, time-consuming and easy to produce additional ambiguous fibers, which make the application of these algorithms difficult[2425]. To reduce the impact of partial volume effects and noise in deterministic algorithms, lighten computational burden, improve tracking fiber length and provide clearer details, this study proposes the tri-linear interpolation algorithm for white matter fiber tracking. Validation of this method demonstrates its suitability for clinical applications.

RESULTS

Tri-linear interpolation algorithm achieved a more complete display of white matter fiber than tensorline algorithm

Fiber tracking on the same regions of interest was performed using both the tri-linear interpolation algorithm and the tensorline algorithm, and their comparisons were made in terms of tracking fiber number, length, tract smoothness and image details. The quantitative analysis in Figure 1A indicates that the total number of tracked fibers by the tri-linear interpolation algorithm was slightly greater than the tensorline algorithm. However, the longest fiber length and the average fiber length of the former method were significantly higher. Tracking results in Figure 2 show that the fibers tracked by the tri-linear interpolation algorithm were distinctly longer, smoother, and had more distinguishable details.

Figure 1

Tracking index comparison of the same regions of interest (corpus callosum of the patient with cerebral infarction) obtained by tri-linear interpolation algorithm and tensorline algorithm.

(A) The diffusion-weighted imaging fractional anisotropy map of the patient. Regions of interest are shown in red.

(B) The index comparison using tri-linear interpolation algorithm and tensorline algorithm of the same regions of interest (genu of corpus callosum), where the total number of tracked fibers is one of the important indicators to evaluate the effectiveness of the fiber tracking algorithms. Lsd indicates the length standard deviation, i.e., the representative fluctuating amount of the tracking fibers’ length; the wider range of the length is, the larger standard deviation will be. Min and Max denote, respectively, the minimum and the maximum fiber length.

Figure 2

Tracking display comparison of the same regions of interest (the corpus callosum knee) in the patient with cerebral infarction.

(A) Tracking result of the tri-linear interpolation algorithm.

(B) Tracking result of the tensorline algorithm.

The fibers tracked by the tri-linear interpolation algorithm are longer, and the fibers that passed though the corpus callosum are more concentrated and clearer. Red, green, and blue indicate right-left, anterior-posterior and superior-inferior orientations, respectively.

Tracking results obtained by the tri-linear interpolation algorithm were in agreement with the actual disease and cerebral anatomy

Fibers tracked by the tri-linear interpolation algorithm (Figure 3), which passed though the corpus callosum in the axial and sagittal planes, illustrated that the fibers obtained were smoother and longer. Furthermore, when combined with color-coded mapping, the orientations of the fibers were easily distinguishable. The images in the coronal plane show that fiber tracts in the lower right region of the cerebral infarction of the patient's brain were obviously abnormal compared with those in the lower left regions. Figure 4 shows the diffusion-weighted images of the lesion, as well as the coronal view of tracked fibers, which is in good agreement with the expected tracking results.

Figure 3

Different views of fibers tracked by the tri-linear interpolation algorithm (corpus callosum of the patient with cerebral infarction).

(A) The axial plane view of the fibers that passed though the corpus callosum.

(B) The sagittal plane view of the fibers that passed though corpus callosum.

The nerve fibers that passed though the corpus callosum are divergent, and their orientations are more obvious. Red, green, and blue indicate right-left, anterior-posterior and superior-inferior orientations, respectively.

Figure 4

Coronal plane view of the fibers in the corpus callosum of the patient with cerebral infarction.

(A) Acute infarction of the right basal ganglia.

(B) Fiber distribution of the two sides in coronal plane is asymmetric, and the fibers within the yellow box are clearly abnormal, which is in good agreement with the cerebral infarction disease. Red, green, and blue indicate right-left, anterior-posterior and superior-inferior orientations, respectively.

Compared with the tracking result of the actual structure in the coronal view (Figure 5), the tracked fibers matched well with our anatomical understanding of the brain structure. Thus, we can draw the conclusion that this algorithm can track and display white matter nerve fibers better than the conventional approach.

Figure 5

Comparison between coronal anatomy and fiber tracking results of the patient with cerebral infarction.

(A) Actual coronal anatomical diagrams of the brain.

(B) The corresponding tracking fibers. The transverse fiber bundle in red basically coincides with the location and direction of the corpus callosum, and the coded colors and the orientation of the fibers are consistent with the actual anatomical structure. Red, green, and blue indicate right-left, anterior-posterior and superior-inferior orientations, respectively.

DISCUSSION

In terms of diffusion tensor imaging based pixel visualization, the main diffusion tensor direction combined with color-coding ability can also be used to analyze fiber orientation. However, it can only be used to observe the trend of fiber orientations roughly and with a two-dimensional flat panel display. Therefore, users cannot clearly visualize the results. Based on diffusion tensor imaging theory, deterministic methods and probabilistic methods have emerged, where each method can display white matter fibers intuitively and comprehensively with color-coding in three-dimensions, facilitating visualization. Unfortunately, both methods fail to reduce the cumulative error caused by noise and partial volume effects, and also fail to simultaneously compute simply and quickly[232425]. Furthermore, no matter which method is employed, the inherent defects of diffusion tensor imaging that can only indicate one direction mean that neither of them can be appropriately used to solve fiber crossing issues per voxel. Therefore, some scholars[262728] put forward the Q-space and Q-Ball imaging methods, where both methods are model-free and can infer the fiber crossing condition within each voxel. Unfortunately, their imaging conditions are precluded by patient safety guidelines that limit the maximum time rate of change of magnetic fields and thus cannot be easily applied to clinical diagnosis. Ozarsla and Mareci[29] developed the diffusion tensor model based on a higher order tensor method and can also infer the fiber crossing within one voxel, but this is still in an early stage of development with little reports on its clinical usage.

This work applied the tri-linear interpolation algorithm based on classical deterministic approaches, which can reduce the influence of noise and partial volume effects effectively. This method not only preserves the advantages of quick and simple calculations but can also obtain better tracking results. Comparison of the tri-linear interpolation algorithm and the tensorline algorithm indicates that in the same regions of interest, tri-linear interpolation tracks 295 fibers, which is slightly more than the 280 fibers tracked by the tensorline algorithm. The longest fiber length and average fiber length are 191 and 121 mm, which are significantly longer than the 110 and 76 mm of the tensorline algorithm. Tracking images also show that the tracked fibers are longer and smoother with clearer displayed details. Finally, we selected data from a female subject with a right cerebral infarction for fiber tracking. The results show that fiber tracts in the lesions are obviously abnormal, which corresponds with the disease very well. At the same time, the coronal view of the tracked fiber shows that the orientation of the fibers in the corpus callosum is consistent with the anatomical structure. However, the correctness of the white matter nerve fiber tracking has no uniform standard, making tracking accuracy difficult to evaluate. Because of limitations of image resolution, the accuracy of fiber tracking itself is not high. Additionally, fiber tracking is more sensitive to noise and partial volume effects; and although the tri-linear interpolation algorithm can reduce these, it still fails to eliminate these effects. Therefore, future research should focus on improving the image resolution and the signal-to-noise ratio. Finding methods that can comprehensively display the diffusion directions of water molecules and better infer the fibers crossing within each voxel will also improve the effectiveness of this model.

MATERIALS AND METHODS

Design

A computer-aided nerve fiber tracking experiment.

Time and setting

The experiment was performed in 2011 in the Medical Image Evaluation Center of Tianjin Huanhu Hospital, China.

Materials

The diffusion tensor imaging data sets were acquired from a female patient with acute infarction of the right basal ganglia, aged 56 years and weighing 77 kg, 15 hours after onset. We applied a Siemens 3.0 T nuclear magnetic resonance imaging system (Erlangen, Bavaria, Germany) with a birdcage head coil, and held the patient's head in place with sponge pads. The device for calculation and visualization was a DELL Vostro 2200 computer (Dell (China) Co., Ltd., Xiamen, Fujian Province, China).

According to the Medical Institution Regulations issued by the State Council of China[30], the experimental program and risk were conveyed to the patient, and an informed consent form was signed before the experiment.

Methods

Magnetic resonance imaging detection

The magnetic resonance imaging experimental data were acquired using a Siemens 3.0 T nuclear magnetic resonance imaging system with an 8 channel phased birdcage head coil in the Imaging Department of Tianjin Huanhu Hospital, China. The field of vision was 256 mm × 256 mm, dimension of sampling matrix was 128 × 128, 44 slices of continuous scanning, and the thickness of scanning was 4 mm. The echo-planar pulse sequence was applied to the diffusion weighted imaging with the following parameters: time of repetition = 6 000 ms, time of echo = 93 ms,b value = 1 000 s/mm2, 20 diffusion gradients, and 3 repetitions.

Data processing

First, we acquired required data from DICOM files of the patient's diffusion weighted images, and reduced image noise using a Gaussian filter. We achieved the principal eigenvalue and its corresponding eigenvector by calculation of the diffusion tensor matrix of each voxel, and expressed the diffusion gradient using the above principal eigenvalue and its associated eigenvector. To obtain the smoother direction of diffusion gradient, we adopted a tri-linear interpolation algorithm. We then set the step size, obtained the coordinates and diffusion gradient of the next point, and began tracking. The whole data processing flow is shown in Figure 6.

Figure 6

Data processing.

Image processing

All diffusion tensor data processing, fiber tracking and visualization in this experiment were performed by the visualization toolkit with combination of the VC++ software platform. For the sake of displaying the tracking fibers using the broken line form, the smaller step size meant the better smoothness of the fiber, but a heaver calculation burden. The step size was set to be 0.2 mm in this experiment. Fiber direction in each voxel was color-coded: the right-left, anterior-posterior and superior-inferior orientations were coded to be red, green, and blue, respectively. The fiber colors of every step were changing within red, green and blue according to the dominant diffusion direction of each voxel, which yielded a better score for visualizing fiber orientation.

Image preprocessing

To shorten the sampling time, fast image acquisition processes were used in the common clinical diagnosis, resulting in low signal-to-noise ratio and poor quality of images. Moreover, the eddy current of fast imaging and the periodic pulsation of cerebrospinal fluid also led to a lower image quality with white Gaussian noise in diffusion weighted images still present. In addition, diffusion weighted image sequences were susceptible to noise and the least squares method that applied to calculate the diffusion tensor is sensitive to outliers, which also caused the obtained diffusion tensor imaging images to be susceptible to noise[31]. Thus far, with respect to noise elimination methods, mean filtering, median filtering and Gaussian filtering are commonly used[32]. Median filtering is better able to denoise impulse noise, but with a poor ability to denoise Gaussian noise. Although a mean filter yields a better score for reducing Gaussian noise, it is easy to blur the image details; while the Gaussian filter can reduce the Gaussian noise very well, it can also better protect image details. Considering that white Gaussian noise was most prevalent in diffusion weighted imaging, the Gaussian filter was chosen for noise reduction.

Gaussian smoothing is an image smoothing processing method based on the neighborhood weighted average idea, where different locations of the pixels are given different weights[33]. The one-dimensional zero-mean Gaussian function is as follows:

where the Gaussian distribution parameter σ determines the width of Gaussian filter. The two-dimensional zero-mean discrete Gaussian function is always chosen to be the smoothing filter in image processing as follows:

There are two ways to design a Gaussian filter: convolution and the Gaussian template method. According to the separability of the Gaussian function, a two-dimensional Gaussian filter can be achieved by convolving two Gaussian filters of one-dimension successively along the horizontal and vertical directions. The above operations can be accomplished by using a one-dimensional Gaussian function to make two convolutions that transpose to the original image. The Gaussian template method designs a Gaussian filter directly by calculating the template weighted value from the discrete Gaussian distribution. For calculation simplicity, the filter weights are always integers, and filtering template weights must be standardized so that uniform grayscale regions of the image are guaranteed not to be affected. The Gaussian filter was chosen in this study, where an appropriate template size and σ are the key points in the design. Because of the two-dimensional Gaussian function, when r > 3σ, G < 0.01, it is generally preferable to choose the filter whose width is less than 2σ2, that is the template width m=2×2σ2+1. The larger Gaussian model involves a wide weighted range, which can easily blur the image; here a Gaussian template with the σ2=1/2, and m=2×2×(1/2)+1=3 is chosen for equation (3).

Figure 7 illustrates the difference between a fractional anisotropy image processed by Gaussian filtering and that which contains noise. It is obvious that the white Gaussian noise (the randomly distributed noise points in Figure 7A) was basically removed and the image details were well preserved.

Figure 7

Comparison of the Gaussian filtering effect.

(A) Fractional anisotropy image with Gaussian noise, where the condition polluted by the Gaussian noise can be seen.

(B) Fractional anisotropy image after Gaussian filtering, where noise in figure (left) is basically removed and the detail preserved.

Acquisition of each voxel gradient

The tri-linear interpolation algorithm is a method to interpolate the gradient of each voxel represented by the principle eigenvalue and associated eigenvector of the diffusion tensor matrix. For the diffusion-weighted imaging that using the Stejskal-Tanner pulsed gradient spin echo sequence, diffusion signal intensity is controlled by the diffusion sensitive b factor[34]. The diffusion-weighted signal intensity S and non-diffusion-weighted signal intensity S0 have the following relationship[35]:

where D is the apparent diffusion coefficient and b is the diffusion sensitive factor,

where γ is the gyromagnetic ratio, G is the diffusion gradient pulse intensity and δ is the diffusion gradient pulse duration. The diffusion of water molecules is three-dimensional, assuming that the diffusion tensor is a second-order tensor, which is a 3 × 3 matrix[36]. Because the tensor is symmetric, only six unique elements are required to fully characterize the tensor[37] as shown in equation (6):

The vector vi=(xi yi zi) is used to indicate the direction of the diffusion pulse gradient, and the relationship between the non-diffusion weighted S0 and Si under different gradient direction fields will be derived by substituting equation (5) and (6) into equation (4)38]:

The diffusion tensor is a real symmetric matrix with 6 unknown parameters, which theoretically need only to apply 6 noncollinear diffusion gradient directions to obtain the tensor matrix:

Equation (7) can be expressed as:

Because the gradient number n is larger than 6, the coefficient matrix is not square. To solve equation (8), multiple linear regression was applied, that is using sample observations to estimate the parameters D, D, D, D, D, D, making the residual sum between actual observation value and estimated value to be minimized[39].

where =(e1e2e3) is a matrix composed by eigenvectors of diffusion tensor and Λ is a diagonal matrix formed by the eigenvalue. Assign λ1≥λ2≥λ3, and choose the biggest eigenvalue λ1 and its corresponding eigenvector e1 to be the predominant diffusion gradient of this voxel, then the gradient f(x,y,z) of point P(x,y,z) will be λ1e1.

For the purpose of estimating the diffusion anisotropy, some researchers put forward some parameters to evaluate its characteristics, which are diffusion tensor trace, mean diffusivity, relative anisotropy, fractional anisotropy (FA) and volume ratio[39]. Among them, FA is the main evaluation parameter, which can provide good gray matter contrast, and has a high signal-noise ratio[40]. Meanwhile, FA is also an important parameter to determine whether or not to stop fiber tracking, defined as:

Tri-linear interpolation algorithm

The energy minimization and limiting principle based tri-linear interpolation linearly interpolates points within a voxel three-dimensionally according to the given values at the vertices of the cube grid cell where the sampling point is located[41]. Define the gradient of point M(x,y,z) to be the principle water molecular diffusion gradient, namely λ1(x,y,z)e1. Assuming that the gradient of a point is g(x,y,z), then g(x,y,z)=λ1(x,y,z)e1. In fact, because of limitations by the imaging resolution and slice thickness, the voxel acquired signal cannot include every point evenly, so g(x,y,z) are actually a series of discrete water molecular diffusion gradients. Assuming that P(x,y,z) is a known point within a voxel, x,y,z are then its coordinate values, x+, y+, z+ are the coordinates with integer values closest to point P(x,y,z) in the positive direction, and x-, y-, z- are the coordinates closest to point P(x,y,z) in the negative direction. Suppose that x, y, z is the coordinate difference between P(x,y,z) and P(x-,y-,z-)

First, interpolate along the y axis:

which is the corresponding function value of the point P2(x-,y,z-) in Figure 8, similarly we have:

Figure 8

Schematic diagram of tri-linear interpolation algorithm.

where, a1, a2, b1, b2 are the four gradients that after interpolation respectively, and their location are the red points as shown in Figure 8.

Second, interpolate along the x axis:

where c1, c2 are the gradients after the second interpolation, and their location are the blue points as shown in Figure 8.

Finally, interpolate along the z axis:

The gradient g(x,y,z) of point P(x,y,z) is derived by the tri-linear interpolation algorithm, which can be applied to tracking white matter fiber. The principle eigenvalues and eigenvectors can be obtained by processing the diffusion tensor imaging data, and the coordinates of next point can be found by substituting the interpolation results into the step size formula:

where P(x,y,z) is the known point, P(x,y,z) is the next tracking point, s is the step size, g(x,y,z) is the present interpolation gradient.

Calculation according to the above method, and, combined with the tracking termination criteria, such as fiber length range and the largest defection angle, white mater fiber tracking can be realized.

Tensorline algorithm

Tensorline algorithm calculates the next advancing direction by weighting the incoming gradient g, the present gradient g1 and the outgoing gradient g that was detected by the voxel tensor[18]. This process can be written as:

where c1 is the linear factor, and its value is:

and w is the weighted factor, whose value is from 0–1, and its empirical value 0.2 is selected in this experiment. The coordinates of next point can be derived by substituting the result into the step formula. Its tracking procession is the same as the tri-linear interpolation algorithm.

Regions of interest selection and tracking termination criteria

The selection of seed points and tracking termination criteria are very important in the fiber tracking procession. Considering human brain anatomical structure characteristics, the corpus callosum was selected as the regions of interest in this experiment[4243]. The majority of this region is cerebral white matter, and connects the left and the right hemispheres of the brain. Furthermore, this region has a rich distribution of white matter fiber-bundles, which is beneficial to display the tracking results. Select the points located in the sagittal plane as the seed points, since the white matter fibers that connected the left and right sides of the brain all pass though the sagittal plane; the fibers in this area are divergent, and with more obvious fiber orientations. In the fiber tracking processing, because of structural differences between individuals, different diseases have different effects on the patients’ nerve fiber structure, and there is no unified standard about how to control the tracking process accurately at present. In this study, the maximum deflection angle between two adjacent voxels and the minimum fractional anisotropy value are chosen as the tracking termination criteria, that is the maximum deflection angle is 90°, minimum fractional anisotropy value is 0.2, the minimum and maximum fiber length are 20 and 300 mm. The tracking diagram is shown in Figure 9.

Figure 9

Illustration of fiber tracking algorithm.

(A) Fiber tracking terminating conditions diagram. (a) Reach the region whose fractional anisotropy value is too low (shown in dark gray), stop tracking (shown in red point); (b) deflection angle between two adjacent voxels exceeded the threshold, stop tracking; (c) normal tracking circumstance. (B) Tri-linear interpolation tracking diagram. Black arrows indicate the voxel gradient around interpolation points, red arrows represent the gradient calculated by tri-linear interpolation, and blue arrows indicate the tracking trajectory.

White matter fiber tracking procession

We applied Siemens 3.0 T nuclear magnetic resonance imaging system with a birdcage head coil, and fixed the patient's head with sponge pads. We first predefined the regions of interest and selected seed points within the regions of interest. Next, we initiated from the seed points and calculated the diffusion direction using tri-linear interpolation. We defined the step size and calculated the coordinates of the next point, and judged whether to stop tracking. Finally, fiber tracking was performed and the tracking results were displayed. Figure 10 is the tracking flow diagram of the tri-linear interpolation algorithm.

Figure 10

Tri-linear interpolation algorithm tracking flow diagram.

ROI: Region of interest; Y: yes; N: no.

Fiber tracking of the same regions of interest by the tri-linear interpolation algorithm and the tensorline algorithm

Fibers of the same regions of interest were tracked by tri-linear interpolation algorithm and the tensorline algorithm, the tracking indexes were compared and analyzed.

Clinical application of the tri-linear interpolation algorithm

We selected the diffusion weighted imaging data set of a female patient with acute infarction of the right basal ganglia, and chose the corpus callosum as regions of interest. We next compared the coronal view of the tracking result with the actual anatomical structure of the brain.