Abstract

DT-MRI (diffusion tensor magnetic resonance imaging) tractography is a method to determine the architecture of axonal fibers in the central nervous system by computing the direction of the principal eigenvectors obtained from tensor matrix, which is different from the conventional isotropic MRI. Tractography based on DT-MRI is known to need many computations and is highly sensitive to noise. Hence, adequate regularization methods, such as image processing techniques, are in demand. Among many regularization methods we are interested in the median filtering method. In this paper, we extended two-dimensional median filters already developed to three-dimensional median filters. We compared four median filtering methods which are two-dimensional simple median method (SM2D), two-dimensional successive Fermat method (SF2D), three-dimensional simple median method (SM3D), and three-dimensional successive Fermat method (SF3D). Three kinds of synthetic data with different altitude angles from axial slices and one kind of human data from MR scanner are considered for numerical implementation by the four filtering methods.

1. Introduction

DT-MRI tractography is a method of noninvasively tracing neuronal fiber bundles. DT-MRI measures anisotropy per pixel and provides the directional information of eigenvectors relevant for fiber tractography. Since DT-MRI data which eigenvectors are computed from usually contain noise, the calculated principal eigenvector direction assumed to be the fiber direction may be different from the real direction in the voxel. As the propagation becomes longer, these differences in the voxels, how small it is, make the whole computed fiber direction deviate far away from the real fiber direction [1]. Since in vivo field maps are usually corrupted by noise, the tractography problem turns out to be a mathematically “ill-posed” problem which means that the tracking results are very sensitive to perturbations by noise. The mathematical attempt to stabilize the solution is known as regularization [2].

Many approaches have been attempted to stabilize the noise problem. The approaches are to stabilize six or more diffusion weighted tensors. Since a diffusion-weighted image is a scalar image, there are many conventional image processing techniques [35]. The diffusion weighted images make one diffusion tensor for each voxel using Stejskal-Tanner formula. There are conflicting opinions for using more than six diffusion-weighed images to make a diffusion tensor for each voxel [69]. Diffusion tensor regularization shows different aspects than conventional scalar image regularization [10, 11]. In fiber tractography, the information about PEV (principal eigenvector) and FA (fractional anisotropy) from the diffusion tensor for each voxel is required. PEV and FA are stabilized using many assumptions or facts including low curvature, small total variation, and orthogonality of eigenpairs [1215].

Among many regularization techniques we are interested in the median filtering of diffusion tensor data [1618]. It is known that the median filtering for anisotropic tensor data preserves the good property of denoising and structure-preserving, which is well known for median filtering of isotropic data [18]. In addition, Kwon et al. devised a successive Fermat filtering for tensor data [19]. Tensor-valued median filter is computed using minimization problem such as gradient descent method or Newton type method in general. Unlike the conventional minimization algorithm, more simple methods, simple median (SM) and successive Fermat (SF) methods, are suggested in [18, 19], respectively. In the papers, only one slice filtering method, that is to say, two-dimensional method for three-dimensional tensor, is considered. In this regard, we call these two-dimensional methods SM2D and SF2D, respectively. But the real DT-MR tensor data are composed of numerous axial, coronal, or sagittal slices. In this paper, we extend SM2D and SF2D to three-dimensional methods SM3D and SF3D, respectively, finding median tensor from neighbor tensors. That is to say, the SF3D/SM3D median tensor is the Fermat/simple median of the three SF2D/SM2D medians from the three consecutive slices.

This paper is organized as follows. Section 2 describes the three-dimensional extension (SM3D and SF3D) of the two-dimensional median algorithms (SM2D and SF2D). In Section 3, we generate the synthetic data and describe MR data for numerical simulations of DT-MR tractography. And we also explained the error measures to test the numerical experiment. In Section 4, a performance characterization of the four-median filtering is presented. Section 5 is devoted to concluding remarks.

2. 3D Median Filtering Methods

In this section, we describe general median filtering and explain briefly SM2D and SF2D. We also extend the two-dimensional methods to the three-dimensional methods: SM3D and SF3D. In this paper, we use the following Frobenius norm as a matrix norm of a matrix : Let be an admissible set of matrices and let be any set of matrices. Then the median of the set is defined as In general median filtering method, is chosen as the set of symmetric positive-definite matrices. The neighboring set can be chosen arbitrarily in general, but we specify just two neighboring sets and for two-dimensional and three-dimensional cases. That is to say, in two-dimensional methods, has chosen nine neighboring matrices in each axial (or sagital or coronal) slice as in Figure 1(a) and, in three-dimensional cases, has chosen 27 neighboring matrices which is 9 matrices for every consecutive three axial slices as in Figure 1(b). After the median filtering, , replaces in two-dimensional methods and in three-dimensional methods as in Figure 1, the Frobenius norm matrix median filter defined in (2) preserves the following properties [1618].(1)If all matrices in are symmetric, scalar-valued, and positive-definite, then is also symmetric, scalar-valued, and positive-definite, respectively.(2)If is Frobenius matrix median of , and are Frobenius matrix median of and , respectively, where is a given positive constant and is a matrix rotating all vectors in three dimensions by angle for some given orientation.

2.1. Simple Median Filtering Methods (SM2D and SM3D)

The SM2D method is to find a simple approximation of the median as in the following equation: It is well known that is contained in the convex hull of the set neighboring set , whereas is a member of , which is one of the nine vertices , , of the convex hull. We extend the SM filter to three dimensions. While the two-dimensional SM filter uses a mask for filtering, we propose three-dimensional SM (SM3D) filtering using a sized mask that covers above ()th slice and below ()th slice in addition to the current th slice, 9 tensors for each. Consequently, the neighboring matrix set is defined as as in Figure 1(b). More explicitly, SM3D is to find the solution of the following minimization problem: and SM3D filtering returns better results than SM2D filtering, because it uses all the three-dimensional neighboring tensor matrices; however, it is more computationally intensive. The SM3D method is explained in Figure 1.

2.2. Successive Fermat Median Filtering Methods (SF2D and SF3D)

Fermat median filtering algorithm proposed by Kwon et al. [19] is extended to three dimensions. Fermat median filtering method is based on the Fermat-Torricelli problem, which is raised by Pierre de Fermat in 1643 to Evangelista Torricelli:

“Given three points in the two-dimensional plane, find the point having the minimal sum of distances to these three points.”

The above problem is the same as the minimizing solution of (2) for three two-dimensional vectors, since any three points in Euclidean space form a two-dimensional hyperplane in the Euclidean space. Several solutions and generalizations are given in [2022]. The algorithm to find Fermat point is obtained by the following properties (see Figure 2).(1)If one of the three angles in the triangle is greater than or equal to 120°, the Fermat point is the vertex at that angle (Figure 2(a)).(2)If all the three angles in the triangle are less than 120°, the Fermat point is the intersection of the two straight lines joining any vertex of the triangle and its symmetrical point. The symmetrical point of , for example, is computed as follows: the symmetrical point of a point of the triangle is chosen by making equilateral triangle at the outside of the triangle but contacting the opposite side of the point (Figure 2(b)).

We call the algorithm to find Fermat point for three tensors using the above properties, Fermat algorithm. Then, SF2D is a method to find the approximation of for neighboring set by using Fermat algorithms successively as follows: For the illustration of (5), see Figure 3(a). The error between and is analyzed in [19].

To define a three-dimensional SF3D, we divided 27 points in into three groups with nine tensors in three consecutive slices. A tensor computed by SF3D is the Fermat point of the three tensors which are SF2D solutions for each slice. Let us assume that in is contained in th axial slice. The SF3D method is depicted in Figure 2(b) and mathematically it is formulated as follows: See Figure 3(b) for the illustration of SF3D.

3. Numerical Experiments

To evaluate the four filtering methods, we prepared three kinds of synthetic data depending on the angle deviating from each axial slice and human data from the MR scanner. These data are three-dimensional.

The flowchart of our numerical experiment is given in Figure 4. Given synthetic or human data, diffusion tensor is computed. In order to compare the filtering methods with respect to noise level, we added Gaussian noise with mean 0 and standard deviation = 0.001, 0.01, and 0.1 to the original data. The noisy images were reconstructed using the four methods: SM2D, SM3D, SF2D, and SF3D. The eigenvalues (, , ) and eigenvectors were calculated by using the power method [23]. Fractional anisotropy (FA) of the tensor for each voxel is given in relation to the three eigenvalues by

The fiber tracking was based on the fiber assignment continuous tracking (FACT) algorithm and a brute-force reconstruction approach [24]. Following the analysis in [1, 24], the fiber tracking was started at the center of an every voxel with FA value greater than 0.3, and terminated at the voxel with an FA less than 0.3 and a tract turning-angle less than 70° in our numerical implementation. Computing FA value, FACT algorithm, and brute-force approach are done using DTI-Studio (CMRM, Johns Hopkins Medical Institute, USA).

3.1. Error Measures

Error measures used in the analysis of numerical examples are also given in this section. Let us define the followings:: The number of voxels: The set of original neighboring tensors: The set of tensors disturbed by Gaussian noise with mean and variance : The set of tensors whose error from original tensor set should be computed.In this paper, can be reconstructed tensor sets using one of the four filtering method or . The error measures used in this paper are the average angular error (AAE) and the average fractional anisotropy (AFA) defined as follows:

3.2. Synthetic Data Generation

The tensor data for each voxel was constructed using the following diagonalization with a diagonal matrix having three eigenvalues and an orthogonal matrix parameterized by azimuth angle and altitude angle as follows: We designed PEV, the first eigenvector, of to have azimuth angle and altitude angle , since the first, second, and third row of corresponds to the corresponding eigenvector for the first, second, and third diagonal elements, respectively. The second and third row of is chosen to make orthogonal matrix.

The number of voxels for the synthesized data was chosen just to verify the efficiency of the filtering in the middle slice for three consecutive slices. The generalization to more slices is not so difficult. Let be the position of the voxel in the xy-plane.

We have considered three kinds of synthetic data, with different altitude angles . We set azimuth angle and make synthetic data zero matrices except 19 continuous diagonal lines as follows: where O is the zero matrix.

The synthesized data and their projection to xy- and yz-plane are depicted in Figure 5.

3.3. Human Data from MR Scanner

The whole procedure obtaining fiber tracts for DT-MR images from the MR scanner data is as follows. The single-shot spin echo is used in the image acquisition and data preprocessing from MR scanner. Echo planar imaging (SE-EPI) pulse sequence with two diffusion sensitizing gradients placed on both sides of the 180° refocusing pulse. Fifty contiguous DT-MR images were obtained at 1.5 T Philips Gyroscan MR scanner with the following imaging parameters: field of view = 224 224 mm2, slice thickness = 3 mm, acquisition matrix = 96, reconstruction matrix = 128, TR = 10,000 ms, TE = 76 ms, factor = 1000 , the number of diffusion sensitizing gradients = 6.

To correct subject head motion and the image distortion due to eddy current, every DTI 3D volume image was realigned to image in FSL (Analysis Group, FMRIB, Oxford, UK). The diffusion tensor was calculated from the 7- volume images (six diffusion weighted images and one image with ). To assess the effects of median filtering on the regularization of the noise in DT-MRI tractography, a noisy data was generated by adding the Gaussian noise with 0 mean and 0.0001 standard deviation to the original human DTI data. Fiber tracking for the original data and noisy data with no filtering, SM3D, SF3D, SM2D and SF2D was performed using the DTI-Studio (CMRM, Johns Hopkins Medical Institute, USA). To select the corticospinal tract (CST), a ROI (Region Of Interest) was drawn on the known anatomical CST area in the pons as shown in Figure 6 and the fibers passing through the ROI were considered as the CST.

4. Numerical Results and Discussion

The numerical results about the three synthetic data in Section 3.2, are given in Figures 7, 8, and 9. Gaussian noise with 0 mean and various standard deviations = 0.1, 0.01, 0.001 are added to the synthesized data in (11). And we compared the performance of the four median filtering regularization algorithms (SM2D, SF2D, SM3D, and SF3D) with respect to AAE and AFA errors in Figure 7. As shown in the figure, three-dimensional filters are superior to two-dimensional filters and SF3D/SF2D is superior to SM3D/SM2D. It is remarkable that SF2D is superior to SM3D in the figure. This shows that the simple median is not so close the median point, whereas the successive median point is much close to the median point even in two-dimension. The xy-plane projection as in Figure 5(a) of the original data , the disturbed image with Gaussian noise having 0 mean and standard deviation 0.01, and reconstructed images by the four methods are pictured in Figure 8. The reconstructed image in Figure 8(c) from SM2D is slightly improved from the disturbed image in Figure 8(b), but still many distorted and broken lines are observed. The distortion and breakage of lines from the original image Figure 8(a) is more improved in Figure 8(d) (SF2D) or Figure 8(e) (SM3D) than the image from Figure 8(b) (SM2D). However, most improved image is found in Figure 8(f), the reconstructed image using SF3D.

In Section 3.2, we considered three kinds of synthetic data , , with different altitude angles . In Figure 9, the reconstructed image using the four methods are compared for the three kinds of synthetic data. In the -axis label, the altitude angle in degree of the three synthetic data is considered. Added Gaussin noise in the disturbed image has zero mean and . The two-dimensional filters show no change of errors with respect to altitude angle. However, three-dimensional filters decreased when altitude angle increases: this phenomenon is observed not so outstandingly for SM3D, but very remarkably for SF3D. The increase of altitude angle means that fiber tract strays out of xy-plane and propagates in the direction close to the -axis. Due to this fact, three dimensional filters shows more sensitive to the altitude angle whereas two-dimensional filters shows no change of the errors with respect to the altitude angle.

In Figures 10, 11, and 12, the reconstruction for human data is considered. Gaussian noise with 0 mean and 0.0001 standard deviation are added. As shown in Figure 10, SF methods show better results than SM methods and three-dimensional methods are better than two-dimensional methods. It is remarkable that SF2D make smaller error than SM3D as in Figure 7. In Figures 11(a)11(c) show the CSTs from the original DTI data and Figures  11(d)11(f), 11(h) and 11(i) from noisy ones with no median filtering, SM3D and SF3D, SM2D and SF2D, respectively. And the green rectangular regions in Figures 11(b) and 11(h) were magnified and arranged on the left and the right in Figure 11(g), respectively. The images in the bottom left corners in Figure 11 show the CSTs in the axial planes represented by the horizontal line in Figure 11(a). Not knowing the true CST, we couldn’t tell which one among Figures 11(a)11(c) was closest to the true CSTs. Therefore, Figures 11(a)11(c) were used as the standards to assess the effect of the median filtering in the CST fiber tracking from noisy data: the CSTs in Figures 11(d)11(f) were compared with Figures 11(a)11(c), respectively. In case of no filtering, the number of fibers in the CST in Figure 11(d) is considerably decreased, while the CST in Figure 11(d) shows a similar shape to the one in Figure 11(a). For the SM3D filtering, the CST in Figure 11(e) shows distinct differences in the region enclosed by the blue circles: In the axial image, the fiber tracts indicated by two yellow arrows in Figure 11(b) disappear in Figure 11(e). For the SF3D filtering, although the CSTs in Figures 11(c) and 11(f) show a few differences, they are similar in shape and the fiber tracts in the axial images show little difference. In comparison with 2-dimensional filters, SF3D Figure 11(f) is closer to Figure 11(c) than SF2D Figure 11(i). On the other hand, SM2D Figure 11(h) appears to be closer to Figure 11(b) than SM3D Figure 11(e), however, as one can see in Figure 11(g), the fibers in SM2D (right) show complicate fiber connections by comparison with the left image, which was not realistic and caused from the poor regularization of the noise in the region. Therefore, the SF3D filter shows best regularization of the noise in CST tracking.

The calculation time average and standard deviation for the four methods are graphed in Figure 12. Intel Pentium with 3 GHz CPU and 2 GB RAM was used for the simulation. The difference of calculation times for SM2D and SF2D is not so remarkable, as we know how exactly SF2D recover fiber tracts than SM2D from Figures 7 and 11. And the difference of computation times between SF2D and SF3D is not so outstanding, since we also know the efficiency of SF3D from previous figures. That is to say, even though SF3D needs a little more computation time than SF2D, SF3D recover DT-MR image far better than SF2D. However, SM3D needs too much time than the other three methods and recovering error is larger than SF2D and slightly better than SM2D.

5. Conclusion

In this study, we developed three-dimensional median filters SM3D and SF3D, extending previously developed SM2D and SF2D in [19]. We implemented three kinds of synthetic data with different altitude angle deviating from the axial slices and one human data form MR scanner. As altitude angle increase from the xy-plane we observed that three dimensional median filters more efficient than two-dimensional filters. For all four synthetic and human data, the four median filters (SM2D, SF2D, SM3D, and SF3D) have same tendency with respect to AAE and AFA errors: three-dimensional filters show superior results to corresponding two-dimensional median filters and SF filters show better reconstruction than SM filters. With respect to computation time, SF2D or SF3D filters need not so much time compared to SM2D filter but SM3D filter needs far more computation time than the other three filters. Therefore, SF3D is proved to be the most efficient median filters and is supposed to be one of powerful regularization methods.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research was supported by the Leading Foreign Research Institute Recruitment Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (2010-00757), and by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2013R1A1A2010624).