# Regularization of DT-MRI Using 3D Median Filtering Methods

Authors: Soondong Kwon; Dongyoun Kim; Bongsoo Han; Kiwoon Kwon
Journal: Journal of Applied Mathematics
 (2014)

---

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


---

## Body

## 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]. Sincein 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 [3–5]. 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 [6–9]. 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 [12–15].Among many regularization techniques we are interested in the median filtering of diffusion tensor data [16–18]. 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 3
×
3
×
3 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. Section2 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 matrixA
=
{
a
j
k
}
j
,
k
=
1
,
…
,
3:
(1)
∥
A
∥
=
∑
j
,
k
=
1
3
a
j
k
2
.
Let B be an admissible set of matrices and let S
=
{
A
1
,
…
,
A
n
} be any set of n matrices. Then the median A
Median of the set S is defined as
(2)
A
Median
=
argmin
X
∈
B
⁡
E
S
(
X
)
,
E
S
(
X
)
=
∑
i
=
1
n
∥
X
-
A
i
∥
.
In general median filtering method, B is chosen as the set of symmetric positive-definite 3
×
3 matrices. The neighboring set S can be chosen arbitrarily in general, but we specify just two neighboring sets S
2 and S
3 for two-dimensional and three-dimensional cases. That is to say, in two-dimensional methods, S
2 has chosen nine neighboring matrices in each axial (or sagital or coronal) slice as in Figure 1(a) and, in three-dimensional cases, S
3 has chosen 27 neighboring matrices which is 9 matrices for every consecutive three axial slices as in Figure 1(b). After the median filtering, A
Median, replaces A
5 in two-dimensional methods and A
14 in three-dimensional methods as in Figure 1, the Frobenius norm matrix median filter A
Median defined in (2) preserves the following properties [16–18].(1)
If all matrices inS are symmetric, scalar-valued, and positive-definite, then A
Median is also symmetric, scalar-valued, and positive-definite, respectively.
(2)
IfA
Median is Frobenius matrix median of S, k
A
Median and B
θ
A
Median are Frobenius matrix median of k
S and B
θ
S, respectively, where k is a given positive constant and B
θ is a matrix rotating all vectors in three dimensions by angle θ for some given orientation.Figure 1
Computations of (a) SM2D and (b) SM3D medians.


### 2.1. Simple Median Filtering Methods (SM2D and SM3D)

The SM2D method is to find a simple approximationA
SM
2
D of the median A
Median as in the following equation:
(3)
A
SM
2
D
=
argmin
X
∈
S
2
⁡
E
S
2
.
It is well known that A
Median is contained in the convex hull of the set neighboring set S
2, whereas A
SM
2
D is a member of S
2, which is one of the nine vertices A
i, i
=
1,2
,
…
,
9, of the convex hull. We extend the SM filter to three dimensions. While the two-dimensional SM filter uses a 3
×
3 mask for filtering, we propose three-dimensional SM (SM3D) filtering using a 3
×
3
×
3 sized mask that covers above (n
-
1)th slice and below (n
+
1)th slice in addition to the current nth slice, 9 tensors for each. Consequently, the neighboring matrix set S
3 is defined as {
A
1
,
A
2
,
A
3
,
…
,
A
27
} as in Figure 1(b). More explicitly, SM3D is to find the solution of the following minimization problem:
(4)
A
SM
3
D
=
argmin
X
∈
S
3
⁡
E
S
3
(
X
)
,
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 [20–22]. 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 (Figure2(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 ofA
1, for example, is computed as follows: the symmetrical point (
A
1
′
) of a point (
A
1
) of the triangle (
Δ
A
1
A
2
A
3
) is chosen by making equilateral triangle (
Δ
A
1
′
A
2
A
3
) at the outside of the triangle but contacting the opposite side (
A
2
A
3
-
) of the point (
A
1
) (Figure 2(b)).Computing Fermat pointF (a) when one of the three angles is greater than 120
° and (b) when all the three angles are less than 120
°. In the latter case, the Fermat point is the intersection of the two straight lines A
1
A
1
′
- and A
3
A
3
′
-.
(a)
(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 ofA
Median for neighboring set S
2 by using Fermat algorithms successively as follows:
(5)
A
SF
2
D
=
argmin
X
∈
B
⁡
E
SF
(
X
)
,
SF
=
{
SF
1
,
S
F
2
,
S
F
3
}
,
SF
i
=
argmin
X
∈
B
⁡
E
F
i
(
X
)
,
F
i
=
{
A
3
(
i
-
1
)
+
1
,
A
3
(
i
-
1
)
+
2
,
A
3
(
i
-
1
)
+
2
}
,
i
=
1,2
,
3
.
For the illustration of (5), see Figure 3(a). The error between A
Median and A
SF
2
D is analyzed in [19].Figure 3
Computations of (a) SF2D and (b) SF3D medians.To define a three-dimensional SF3D, we divided 27 points inS
3 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 A
14 in S
3 is contained in nth axial slice. The SF3D method is depicted in Figure 2(b) and mathematically it is formulated as follows:
(6)
A
SF
3
D
=
argmin
X
∈
B
⁡
E
SF
(
X
)
,
SF
=
{
SF
n
-
1
,
SF
n
,
SF
n
+
1
}
,
SF
m
=
argmin
X
∈
B
⁡
E
F
m
(
X
)
,
m
=
n
-
1
,
n
,
n
+
1
,
F
m
=
{
A
9
(
m
-
n
+
1
)
+
i
}
i
=
1,2
,
…
,
9
.
See Figure 3(b) for the illustration of SF3D.


## 2.1. Simple Median Filtering Methods (SM2D and SM3D)

The SM2D method is to find a simple approximationA
SM
2
D of the median A
Median as in the following equation:
(3)
A
SM
2
D
=
argmin
X
∈
S
2
⁡
E
S
2
.
It is well known that A
Median is contained in the convex hull of the set neighboring set S
2, whereas A
SM
2
D is a member of S
2, which is one of the nine vertices A
i, i
=
1,2
,
…
,
9, of the convex hull. We extend the SM filter to three dimensions. While the two-dimensional SM filter uses a 3
×
3 mask for filtering, we propose three-dimensional SM (SM3D) filtering using a 3
×
3
×
3 sized mask that covers above (n
-
1)th slice and below (n
+
1)th slice in addition to the current nth slice, 9 tensors for each. Consequently, the neighboring matrix set S
3 is defined as {
A
1
,
A
2
,
A
3
,
…
,
A
27
} as in Figure 1(b). More explicitly, SM3D is to find the solution of the following minimization problem:
(4)
A
SM
3
D
=
argmin
X
∈
S
3
⁡
E
S
3
(
X
)
,
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 [20–22]. 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 (Figure2(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 ofA
1, for example, is computed as follows: the symmetrical point (
A
1
′
) of a point (
A
1
) of the triangle (
Δ
A
1
A
2
A
3
) is chosen by making equilateral triangle (
Δ
A
1
′
A
2
A
3
) at the outside of the triangle but contacting the opposite side (
A
2
A
3
-
) of the point (
A
1
) (Figure 2(b)).Computing Fermat pointF (a) when one of the three angles is greater than 120
° and (b) when all the three angles are less than 120
°. In the latter case, the Fermat point is the intersection of the two straight lines A
1
A
1
′
- and A
3
A
3
′
-.
(a)
(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 ofA
Median for neighboring set S
2 by using Fermat algorithms successively as follows:
(5)
A
SF
2
D
=
argmin
X
∈
B
⁡
E
SF
(
X
)
,
SF
=
{
SF
1
,
S
F
2
,
S
F
3
}
,
SF
i
=
argmin
X
∈
B
⁡
E
F
i
(
X
)
,
F
i
=
{
A
3
(
i
-
1
)
+
1
,
A
3
(
i
-
1
)
+
2
,
A
3
(
i
-
1
)
+
2
}
,
i
=
1,2
,
3
.
For the illustration of (5), see Figure 3(a). The error between A
Median and A
SF
2
D is analyzed in [19].Figure 3
Computations of (a) SF2D and (b) SF3D medians.To define a three-dimensional SF3D, we divided 27 points inS
3 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 A
14 in S
3 is contained in nth axial slice. The SF3D method is depicted in Figure 2(b) and mathematically it is formulated as follows:
(6)
A
SF
3
D
=
argmin
X
∈
B
⁡
E
SF
(
X
)
,
SF
=
{
SF
n
-
1
,
SF
n
,
SF
n
+
1
}
,
SF
m
=
argmin
X
∈
B
⁡
E
F
m
(
X
)
,
m
=
n
-
1
,
n
,
n
+
1
,
F
m
=
{
A
9
(
m
-
n
+
1
)
+
i
}
i
=
1,2
,
…
,
9
.
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 Figure4. 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 (λ
1, λ
2, λ
3) and eigenvectors were calculated by using the power method [23]. Fractional anisotropy (FA) of the tensor A for each voxel is given in relation to the three eigenvalues by
(7)
FA
(
A
)
=
3
(
λ
1
-
m
)
2
+
(
λ
2
-
m
)
2
+
(
λ
3
-
m
)
2
2
(
λ
1
2
+
λ
2
2
+
λ
3
2
)
,
m
=
λ
1
+
λ
2
+
λ
3
3
.Figure 4
Flowchart of the fiber tractography loaded from synthetic data (Section3.2) or DTI data (Section 3.3).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:N
: The number of voxels
A
=
{
A
n
}
n
=
1
,
…
,
N: The set of original neighboring tensors
A
σ: The set of tensors disturbed by Gaussian noise with mean 0 and variance σ
2
X: The set of tensors whose error from original tensor set A should be computed.In this paper, X can be reconstructed tensor sets using one of the four filtering method or A
σ. The error measures used in this paper are the average angular error (AAE) and the average fractional anisotropy (AFA) defined as follows:
(8)
AAE
(
X
:
A
)
=
1
N
∑
n
=
1
N
cos
⁡
-
1
⁡
(
PEV
(
X
n
)
·
PEV
(
A
n
)
)
,
FA
(
X
:
A
)
=
1
N
∑
n
=
1
N
|
FA
(
X
n
)
-
FA
(
A
n
)
|
.


### 3.2. Synthetic Data Generation

The tensor data for each voxel was constructed using the following diagonalization with a diagonal matrix having three eigenvalues(
1,0.2,0.2
) and an orthogonal matrix parameterized by azimuth angle θ and altitude angle φ as follows:
(9)
A
θ
,
φ
=
T
θ
,
φ
t
(
1
0
0
0
0.2
0
0
0
0.2
)
T
θ
,
φ
,
T
θ
,
φ
=
(
cos
⁡
ϕ
cos
⁡
θ
cos
⁡
ϕ
sin
θ
sin
ϕ
sin
ϕ
cos
⁡
θ
sin
ϕ
sin
θ
-
cos
⁡
ϕ
sin
θ
-
cos
⁡
θ
0
)
.
We designed PEV, the first eigenvector, of A
θ
,
ϕ to have azimuth angle θ and altitude angle ϕ, since the first, second, and third row of T
θ
,
ϕ corresponds to the corresponding eigenvector for the first, second, and third diagonal elements, respectively. The second and third row of T
θ
,
ϕ is chosen to make T
θ
,
ϕorthogonal matrix.The number of voxels for the synthesized data was chosen128
×
128
×
3 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 (
i
x
,
i
y
) be the position of the voxel in thexy-plane.We have considered three kinds of synthetic data, with different altitude anglesϕ
=
30
°
,
45
°
,
60
°. We set azimuth angle θ
=
45
° and make synthetic data zero matrices except 19 continuous diagonal lines as follows:
(10)
A
I
=
{
A
π
/
4
,
π
/
6
=
(
0.500
0.300
0.245
0.300
0.500
0.245
0.245
0.245
0.400
)
,
i
y
-
10
<
i
x
<
i
y
+
10
,
O
,
otherwise
,
(11)
A
I
I
=
{
A
π
/
4
,
π
/
4
=
(
0.400
0.200
0.283
0.200
0.400
0.283
0.283
0.283
0.600
)
,
i
y
-
10
<
i
x
<
i
y
+
10
,
O
,
otherwise
,
(12)
A
I
I
I
=
{
A
π
/
4
,
π
/
3
=
(
0.500
0.300
0.245
0.300
0.500
0.245
0.245
0.400
0.400
)
,
i
y
-
10
<
i
x
<
i
y
+
10
,
O
,
otherwise
,
where O is the zero matrix.The synthesized data and their projection toxy- andyz-plane are depicted in Figure 5.Synthesized dataA
I, A
I
I, and A
I
I
I given in (10), (11), and (12): (a) three-dimensional representation of A
I and its projection toxy-plane, (b) another three-dimensional representation of A
I and its projection toyz-plane, and resulting two-dimensional representation of the synthetic data (c) A
I, (d) A
I
I, and (e) A
I
I
I.
(a)
(b)
(c)
(d)
(e)


### 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, b factor = 1000 s
/
mm
2, 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 tob
0 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 b
=
0). 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.MRI data for the pons in the (a) sagittal and (b) axial planes. The rectangle in (b) was chosen as ROI to select the corticospinal tract (CST).
(a)
(b)


## 3.1. Error Measures

Error measures used in the analysis of numerical examples are also given in this section. Let us define the followings:N
: The number of voxels
A
=
{
A
n
}
n
=
1
,
…
,
N: The set of original neighboring tensors
A
σ: The set of tensors disturbed by Gaussian noise with mean 0 and variance σ
2
X: The set of tensors whose error from original tensor set A should be computed.In this paper, X can be reconstructed tensor sets using one of the four filtering method or A
σ. The error measures used in this paper are the average angular error (AAE) and the average fractional anisotropy (AFA) defined as follows:
(8)
AAE
(
X
:
A
)
=
1
N
∑
n
=
1
N
cos
⁡
-
1
⁡
(
PEV
(
X
n
)
·
PEV
(
A
n
)
)
,
FA
(
X
:
A
)
=
1
N
∑
n
=
1
N
|
FA
(
X
n
)
-
FA
(
A
n
)
|
.


## 3.2. Synthetic Data Generation

The tensor data for each voxel was constructed using the following diagonalization with a diagonal matrix having three eigenvalues(
1,0.2,0.2
) and an orthogonal matrix parameterized by azimuth angle θ and altitude angle φ as follows:
(9)
A
θ
,
φ
=
T
θ
,
φ
t
(
1
0
0
0
0.2
0
0
0
0.2
)
T
θ
,
φ
,
T
θ
,
φ
=
(
cos
⁡
ϕ
cos
⁡
θ
cos
⁡
ϕ
sin
θ
sin
ϕ
sin
ϕ
cos
⁡
θ
sin
ϕ
sin
θ
-
cos
⁡
ϕ
sin
θ
-
cos
⁡
θ
0
)
.
We designed PEV, the first eigenvector, of A
θ
,
ϕ to have azimuth angle θ and altitude angle ϕ, since the first, second, and third row of T
θ
,
ϕ corresponds to the corresponding eigenvector for the first, second, and third diagonal elements, respectively. The second and third row of T
θ
,
ϕ is chosen to make T
θ
,
ϕorthogonal matrix.The number of voxels for the synthesized data was chosen128
×
128
×
3 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 (
i
x
,
i
y
) be the position of the voxel in thexy-plane.We have considered three kinds of synthetic data, with different altitude anglesϕ
=
30
°
,
45
°
,
60
°. We set azimuth angle θ
=
45
° and make synthetic data zero matrices except 19 continuous diagonal lines as follows:
(10)
A
I
=
{
A
π
/
4
,
π
/
6
=
(
0.500
0.300
0.245
0.300
0.500
0.245
0.245
0.245
0.400
)
,
i
y
-
10
<
i
x
<
i
y
+
10
,
O
,
otherwise
,
(11)
A
I
I
=
{
A
π
/
4
,
π
/
4
=
(
0.400
0.200
0.283
0.200
0.400
0.283
0.283
0.283
0.600
)
,
i
y
-
10
<
i
x
<
i
y
+
10
,
O
,
otherwise
,
(12)
A
I
I
I
=
{
A
π
/
4
,
π
/
3
=
(
0.500
0.300
0.245
0.300
0.500
0.245
0.245
0.400
0.400
)
,
i
y
-
10
<
i
x
<
i
y
+
10
,
O
,
otherwise
,
where O is the zero matrix.The synthesized data and their projection toxy- andyz-plane are depicted in Figure 5.Synthesized dataA
I, A
I
I, and A
I
I
I given in (10), (11), and (12): (a) three-dimensional representation of A
I and its projection toxy-plane, (b) another three-dimensional representation of A
I and its projection toyz-plane, and resulting two-dimensional representation of the synthetic data (c) A
I, (d) A
I
I, and (e) A
I
I
I.
(a)
(b)
(c)
(d)
(e)


## 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, b factor = 1000 s
/
mm
2, 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 tob
0 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 b
=
0). 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.MRI data for the pons in the (a) sagittal and (b) axial planes. The rectangle in (b) was chosen as ROI to select the corticospinal tract (CST).
(a)
(b)


## 4. Numerical Results and Discussion

The numerical results about the three synthetic data in Section3.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 A
I
Iin (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. Thexy-plane projection as in Figure 5(a) of the original data A
I
I, 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.Errors with respect toσ = 0.1, 0.01, and 0.001 when we recover the synthetic data A
I
I using the four median filtering methods: (a) Log (AAE) error and (b) Log (AFA) error.
(a)
(b)Filtering of synthetic dataA
I
I corrupted with Gaussian noise having zero mean and σ = 0.01: (a) the original image, (b) the disturbed image, and the reconstructed images using (c) SM2D, (d) SF2D, (e) SM3D, and (f) SF3D.
(a)
(b)
(c)
(d)
(e)
(f)Errors with respect toA
I, A
I
I, and A
I
I
I with φ
=
30
°
,
45
°
,
60
°, respectively, when we recover the synthetic data with Gaussian noise having mean 0 and σ = 0.01 using the four median filtering methods: (a) Log (AAE) error and (b) Log (AFA) error.
(a)
(b)In Section3.2, we considered three kinds of synthetic data A
I, A
I
I, A
I
I
I with different altitude angles ϕ
=
30
°
,
45
°
,
60
°. In Figure 9, the reconstructed image using the four methods are compared for the three kinds of synthetic data. In the x-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 σ
=
0.0
1. 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 ofxy-plane and propagates in the direction close to the z-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 Figures10, 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 (
σ
=
0
.
000
1
) 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.Errors when we recover human data from MR scanner corrupted with Gaussian noise having zero mean andσ = 0.0001 by using the four median filtering methods: (a) Log (AAE) error and (b) Log (AFA) error.
(a)
(b)CSTs passing through the pons of a human subject: corticospinal tracts obtained from the original data with (a) no filter, (b) SM3D, and (c) SF3D and the CSTs obtained from the disturbed data (σ = 0.0001) with (d) no filter, (e) SM3D, (f) SF3D, (h) SM2D, and (i) SF2D. The images in the bottom left corners show the CSTs in the axial planes represented by the horizontal line in (a). And the left figure and the right one in (g) are magnified images of the green rectangular regions in (b) and (h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 12
Computation time for the four methods to reconstruct the DT-MR data: the blue solid bar and small vertical bar at the top end of the blue bar represent the average and the standard deviation of computational times for the four methods, respectively, after testing the four median filters five times for each.The calculation time average and standard deviation for the four methods are graphed in Figure12. 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 thexy-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.



---

*Source: 285367-2014-07-07.xml*