Principal eigenvector field segmentation for reproducible diffusion tensor tractography of white matter structures

Available online at www.sciencedirect.com

Magnetic Resonance Imaging 29 (2011) 1088 – 1100

Principal eigenvector field segmentation for reproducible diffusion tensor tractography of white matter structures Ram K.S. Rathore a,⁎, Rakesh K. Gupta b , Shruti Agarwal a , Richa Trivedi b,c , Rajendra P. Tripathi c , Rishi Awasthi b b

a Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Delhi, India Department of Radiodiagnosis, Sanjay Gandhi Postgraduate Institute of Medical Sciences, Lucknow, Delhi, India c NMR Research Centre, Institute of Nuclear Medicine and Allied Sciences, Delhi, India Received 13 September 2010; revised 16 April 2011; accepted 22 April 2011

Abstract The study was aimed to test the feasibility of utilizing an algorithmically determinable stable fiber mass (SFM) map obtained by an unsupervised principal eigenvector field segmentation (PEVFS) for automatic delineation of 18 white matter (WM) tracts: (1) corpus callosum (CC), (2) tapetum (TP), (3) inferior longitudinal fasciculus (ILF), (4) uncinate fasciculus (UNC), (5) inferior fronto-occipital fasciculus (IFO), (6) optic pathways (OP), (7) superior longitudinal fasciculus (SLF), (8) arcuate fasciculus (AF), (9) fornix (FX), (10) cingulum (CG), (11) anterior thalamic radiation (ATR), (12) superior thalamic radiation (STR), (13) posterior thalamic radiation (PTR), (14) corticospinal/corticopontine tract (CST/CPT), (15) medial lemniscus (ML), (16) superior cerebellar peduncle (SCP), (17) middle cerebellar peduncle (MCP) and (18) inferior cerebellar peduncle (ICP). Diffusion tensor imaging (DTI)-derived fractional anisotropy (FA) and the principal eigenvector field have been used to create the SFM consisting of a collection of linear voxel structures which are grouped together by color-coding them into seven natural classes to provide PEVFS signature segments which greatly facilitate the selection of regions of interest (ROIs) for fiber tractography using just a single mouse click, as compared with a manual drawing of ROIs in the classical approach. All the 18 fiber bundles have been successfully reconstructed, in all the subjects, using the single ROIs provided by the SFM approach, with their reproducibility characterized by the fact that the ROI selection is user independent. The essentially automatic PEVFS method is robust, efficient and compares favorably with the classical ROI methods for diffusion tensor tractography (DTT). © 2011 Elsevier Inc. All rights reserved. Keywords: Diffusion tensor imaging; Tractography; Reproducibility; White matter

1. Introduction The widely used diffusion tensor imaging (DTI)-based diffusion tensor tractography (DTT) technique to study human white matter (WM) anatomy [1–9] is impeded by several shortcomings. Factors like noise, partial volume effects and complex fiber architectures within a pixel can lead to production of false-positive and false-negative

⁎ Corresponding author. Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur, UP 208016, India. Tel.: +91 512 2598741; fax: +91 512 2597500. E-mail addresses: [email protected], [email protected] (R.K.S. Rathore). 0730-725X/$ – see front matter © 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2011.04.014

results. A further challenge is placement of reference regions of interest (ROIs). Most common methods of fiber tracking (FT) rely on a knowledge-based selection of ROIs on appropriate slices and then generate the fibers from there followed by cleaning and refining the bundle obtained [6,10–12] which is normally quite time consuming and involves a considerable amount of trial and error. The task of marking the ROIs is significantly more difficult if the structure lies very close to another region like the tapetum. Even in the absence of these difficulties the values and the significance of the parameters used in the selection of ROIs depend heavily on the size, shape, number, etc., of the ROIs marked out [13]. In some semiautomatic methods and in multi-ROI techniques [6,10–12], the accuracy of the results can be increased if the ROIs are made both on the fractional

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

1089

Table 1 Summary of the mean±S.D. of FA values of 18 major fiber tracts of healthy controls (n=20) Fiber tracts

Mean±S.D.

Fiber tracts

Mean±S.D.

Fiber tracts

Mean±S.D.

CC Right TP Left TP Right STR Left STR Right PTR Left PTR Right UNC Left UNC Right CST Left CST Right ML

0.49±0.02 0.47±0.02 0.47±0.02 0.44±0.03 0.44±0.03 0.44±0.02 0.44±0.02 0.37±0.03 0.36±0.03 0.51±0.02 0.52±0.03 0.44±0.02

Left ML Right OP Left OP Crux of FX Right tail of FX Left tail of FX Right AF Left AF Right ATR Left ATR Right IFO Left IFO

0.43±0.03 0.42±0.03 0.40±0.04 0.34±0.03 0.43±0.05 0.41±0.05 0.43±0.03 0.43±0.02 0.39±0.03 0.40±0.03 0.44±0.02 0.43±0.03

Right MCP Left MCP Right SCP Left SCP Right ICP Left ICP Right SLF Left SLF Right ILF Left ILF Right CG Left CG

0.44±0.02 0.43±0.02 0.43±0.04 0.42±0.03 0.50±0.05 0.49±0.05 0.43±0.03 0.42±0.02 0.41±0.03 0.41±0.02 0.39±0.02 0.40±0.02

CC, Corpus callosum; TP, tapetum; ILF, inferior longitudinal fasciculus; UNC, uncinate fasciculus; IFO, inferior fronto-occipital fasciculus; OP, optic pathways; SLF, superior longitudinal fasciculus; AF, arcuate fasciculus; FX, fornix; CG, cingulum; ATR, anterior thalamic radiation; STR, superior thalamic radiation; PTR, posterior thalamic radiation; CST, corticospinal/corticopontine tract; ML, medial lemniscus; SCP, superior cerebellar peduncle; MCP, middle cerebellar peduncle; ICP, inferior cerebellar peduncle.

anisotropy (FA) and on the structural images. Nevertheless, there exists a tradeoff between the time taken to complete the task and the accuracy of the results. Hence, more as a rule, most ROI-based methods are tedious, time-consuming and error prone. Although template-based methods have been tried, they do not have acceptable levels of accuracy [14]. Some recent articles have considered region-based segmentation in DT-MRI using the level set methods [15–18]. The method presented by Wang and Vemuri [16] uses the whole tensor information to define a region-based force as proposed by Paragios and Deriche [19]. WM structures have also been identified using directional correlation-based region growing (DCRG) algorithms [20,21] that use the inner product (IP) of e1 of adjacent voxels, called the directional correlation (DC), to investigate the similarity of fiber tract direction in WM tracts. A major concern of regiongrowing algorithms such as DCRG is that they are iterative in nature and are usually computationally very intensive. The resulting segmentations are usually sensitive to the choice of parameters. Since there is no objective measure of the exact standard for the brain images, it is difficult to determine the optimum threshold exactly. In the absence of any such standard, the establishment of reproducible tracking protocols becomes very important. Principal eigenvector field segmentation (PEVFS) methodology was developed exactly with this need in mind. It utilizes FA and e1 -field maps to construct a presegmented map, called stable fiber mass (SFM). SFM is automatically segmented, color coded and resliced into axial, coronal and sagittal maps containing signature segments of different fiber bundles. These maps are used to automatically delineate the white matter structures by choosing any pixel within a maximal signature segment associated with the structure. Following a click at the pixel, the boundary of the ROI is drawn and the ROI is tracked across slices using fiber assignment by continuous tracking (FACT), followed by cleaning and refinement operations

using a similarity measure to provide the voxel coordinates of fiber bundles. These coordinates constitute the output of the approach and provide masks for the fiber bundles for computing precise statistics on their shape, size, location and distribution of different parameters of clinical interest there, an application of which is provided in Table 1 listing the means and standard deviations (S.D.) of FA values of 18 major fiber tracts of the cohort of the 20 healthy controls of our study.

2. Materials and methods 2.1. Imaging and healthy controls Conventional MRI and DTI were performed on 20 normal human subjects (15 men, 5 women; mean age, 25 years; range, 21–40 years). Imaging was performed on a 3-T MRI scanner (GE Healthcare Technologies, Milwaukee, WI, USA) using a 12-channel head coil. MR imaging parameters included T2-weighted axial fast spin-echo images [repetition time (TR)=4500 ms, echo time (TE)=70 ms, number of excitations (NEX)=1], T1-weighted axial spinecho fluid-attenuated inversion recovery (FLAIR) images (TR=1900 ms, TE=8 ms, NEX=1) and T2-weighted FLAIR images (TR=10,000 ms, TE=140 ms, inversion time=2250 ms, NEX=1). A total of 46 contiguous 3-mm-thick axial sections were acquired with a 240×240-mm2 field of view (FOV) and image matrix of 256×256. DTI data sets were acquired using a single-shot echo-planar dual SE sequence with ramp sampling and acquisition parameters as: b value=0 s/mm2, 1000 s/mm2; slice thickness=3 mm with no gap; FOV=240×240 mm2; TR=12.5 s; TE=86.4 ms; and NEX=1. A total of 46 axial sections were acquired with an image matrix of 256×256 (no zero-filling). The diffusion tensor encoding used was balanced rotationally invariant with 21 noncollinear directions over the unit hemisphere.

1090

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

The DTI data was distortion corrected [22], descalped, interpolated to get isotropic cubic voxel data and then processed to obtain eigenvalues, eigenvectors, FA and other DTI metrics by using in-house-developed methodology at ADISL, IIT Kanpur, India.

2.2. Methodology The PEVFS methodology utilizes four modules: (I) PEVFS; (II) fiber tracking; (III) cleaning; and (IV) refinement, integrated in a reliable, time-efficient, user-

Fig. 1. Flowchart for the PEVFS procedure.

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

Fig. 1 (continued).

1091

1092

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

friendly [23] automatic program through which the DTT is only a few mouse clicks away from a reproducible output, without a need of any considerable expertise in the field. 2.3. Principal eigenvector field segmentation The PEVFS to delineate WM structures is unsupervised. It is based primarily on the following three observations: (1) adjacent regions within a WM tract possess high directional similarity in water diffusion; (2) many WM structures of the brain have their axes of greatest diffusivity along one of the three cardinal directions (i.e., superior–inferior, anterior– posterior or left–right) or along the diagonals of the unit cube (±1, ±1, ±1) or the diagonals of its faces (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1); and (3) the boundaries of these WM structures consist mostly of gray matter (GM), cerebral spinal fluid (CSF) or other WM structures that have their greatest diffusivity oriented in an almost randomly different direction. The basic principle underlying PEVFS methodology is that, for trajectories defined by a direction field, the best points to start a bidirectional tracking are the ones where the curvature of the trajectory is the least: at the start of a trajectory, any method has to utilize the basic relation r(s+h) −[r(s)+hr′(s)] ≅ 1/2 h2r″ (s), h →0, explicitly (e.g., in Euler's method and FACT) or implicitly (e.g., in Runge– Kutta methods for intermediate points). Here r(s+h) (new point) and r(s) (previous point) are the voxel positions on a fiber being tracked, r′(s) is the direction vector (the principal eigenvector e1) and r″(s) is the rate of change of tangential direction vector whose magnitude κ=|r″| is the curvature at a point associated with the arc length variable s of a trajectory. Hence, to minimize the trajectory deviation in a bi-directional tracking the choice of the seed point on the trajectory should correspond to one of a minimal curvature |r″|. The stable voxels in the PEVFS are precisely the voxels having a practically zero curvature. Thus PEVFS basically helps define the most promising starting or seed regions from which the fiber tracking ought to be initiated, leaving the subsequent tracking and addition (refinement) and deletion of fibers in a bundle to any standard practice in the literature (or the FACT algorithm and a similarity measure-based refinement and deletion procedures in our use, described in later sections). 2.4. Center-Euler-center algorithm and stable voxels The center-Euler-center (CEC) algorithm first obtains the e1 of the diffusion tensor at a seed voxel P (i, j, k) with an FA greater than a set threshold (default=0 or 0.15). The next voxel Q (x, y, z) is computed by
 Q = ROUND P + e1P + 0:5u

ð1Þ

Here and in the sequel, eP1 is the PEV of P, u=(1, 1, 1) and the function ROUND stands for the componentwise “integral part” operation. The tracking continues from the center of the

current voxel and reaches the center of the new voxel using Euler's algorithm with step size h=1, as in Eq. (1). The CEC algorithm follows a track from the seed point only in the positive (forward) direction. The process is repeated while FA remains above the threshold, and the number of voxels covered is under a given maximum. We use CEC to extract stable fiber mass, for which the tracking needs be performed only in one direction, in view of which the stopping condition involving the inner product of the eigenvectors of consecutive voxels, as in FACT, is dropped. To check the stability property of voxel P, proceed in 1 of voxel Q, if the backward direction along the PEV eQ the FA value at Q is greater than the threshold to obtain a voxel P′ as:   P V= ROUND Q − e1Q + 0:5u ð2Þ The voxels P and Q are said to form a stable pair iff P′=P. Mathematically, Eq. (1) can be used to define a function G(P) as: Gð P Þ = Q

ð3Þ

Similarly, Eq. (2) could be used to define a function F(Q) by: F ðQ Þ = P

ð4Þ

Combining Eqs. (3) and (4), the stability of (x, y, z) and (i, j, k) translates to the relations: Gð F ð x; y; zÞÞ = ð x; y; zÞ F ðGði; j; k ÞÞ = ði; j; k Þ

ð5Þ ð6Þ

A satisfaction of any of Eqs. (5) and (6) makes the corresponding point (x, y, z) or (i, j, k) a stable voxel. It is clear from the above that a determination of stable voxels in a volume needs only forward CEC tracking as mentioned earlier. 2.5. Modules of PEVFS PEVFS requires four unsupervised modules to segment WM structures: (a) labeling: assigning an integer value to all voxels in the data set such that the voxels associated with the same positive integer value label constitute a string with any two consecutive members forming a stable pair of voxels; (b) coloring: construction of stable fiber mass pattern (SFMP) map by replacing each assigned label in the labeling module a random color showing strings of stable fibers as pencils (Fig. 2I); (c) grouping: construction of SFM map grouping the fiber pencils from the coloring module with similar orientation by the use of an eight-color scheme providing a segmentation of the volume stack (Fig. 2II); and (d) clicking: pointing to a segment of interest in an axial, sagittal or coronal section in the SFM map to

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

automatically delineate the boundary of the signature segment structure and initiate automatic tracking, cleaning and refinement (Fig. 4I–III). The flowchart of PEVFS procedure is shown in Fig. 1. FA is used to distinguish GM from WM and CSF in PEVFS, as the first step. However, the FA threshold is not mandatory for segmentation using PEVFS and can be inactivated by putting the threshold as zero. As GM and CSF voxels, more as a rule, do not exhibit the stability property, even without applying any FA threshold the GM

1093

and CSF voxels will normally be winnowed out. As described in Ref. [24], the water molecules diffuse with a high degree of directional similarity in the vicinity of the same WM tract, making PEVFS inherently capable of distinguishing one WM tract from the other because of the directional similarity of water diffusion. Starting from a seed point within the WM tract of interest, PEVFS progressively examines the ability to form a stable pair with the closest neighboring voxel lying on the same fiber tract. Stable voxels lying on the same track are given the same label. PEVFS allows identification of the selected WM tract of interest without applying any threshold for the FA and it does not use the DC of the PEVs of the adjacent voxels at all. Furthermore, a threshold smaller than the one used in subsequent tracking, e.g., by FACT or another method, will anyway render the corresponding voxels in the ROI incapable of starting a fiber tracking from there. The outcome of the segmented WM tract is therefore essentially independent of any of the moderately selected threshold values of FA; in fact, with no threshold the segmented ROIs look slightly enlarged, making their identification possibly easier (Fig. 2II). 2.5.1. Labeling The following steps are executed for each seed voxel in the dataset with FA greater than the prescribed threshold value (e.g., 0.0, 0.15, etc.): Step 1. In the beginning, each voxel in the data set is assigned a label of value 0. Let label counter=1 Step 2. Choose a voxel P0 having label zero and with FA greater than the prescribed threshold, i.e., Let label (P0)=0 and FA (P0)Ncutoff Step 3. Assign label of P0 as label counter, i.e., label (P0)=label counter, and increase the label counter by 1. Step 4. Start tracking from the seed voxel P0. If A= {P0, P1…Pn} is the fiber path with seed point P0 obtained using the CEC algorithm and if ei1 denotes the PEV at Pi:
 Pi = ROUND Pi − 1 + e1i − 1 + 0:5u ð7Þ Step 5. Check whether Pj and Pj−1 form a stable pair. Starting from j=1,
 ð8Þ Qj = ROUND Pj − e1i + 0:5u If Qj=Pj−1, then Pj and Pj−1 form a stable pair. Here 1≤j≤n. Step 6. If Pj is not a stable point, then start from Step 2 with another voxel from the dataset. Step 7. If Pj is a stable point with label 0, then assign Pj the label of Pj−1:

 ð9Þ label Pj = label Pj−1

Fig. 2. (I) The various geometrical shapes that voxel groups take in SFMP map. (II) Axial, coronal and sagittal sections of normal brain under the SFM map created without any FA threshold (upper row) and corresponding slices in SFM map which is created with FA N0.15 (lower row).

Repeat Steps 5, 6 and 7 for all 1≤j≤n. When all values of j are exhausted, then restart from Step 2 with another voxel of the dataset.

1094

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

Step 8. If Pj is a stable point with nonzero label, then assign Pk the label of Pj
 labelðPk Þ = label Pj ; where 0VkVj−1 ð10Þ Decrease the label counter by 1 and start from Step 2 with another voxel of the data set. Step 9. Repeat Step 2 to Step 8 for all the voxels in the dataset. In this way, each voxel gets a label. The label zero of a voxel indicates low FA at that point (i.e., FA did not exceed the prescribed cutoff). 2.5.2. Coloring: SFMP For viewing purposes, a color image volume is constructed by assigning colors to all the labels according to their values. The colors are generated randomly and voxels having the same value of labels are colored uniformly. This segmented fiber mass pattern color map is named as SFMP map. Sections of the SFMP map show up in a specific characteristic pattern of colors. These patterns are the geometrical shapes of a horizontal bar, a vertical bar, a diagonal bar and a nonhomogenous group of voxels within a slice (Fig. 2I). The bar for each region is displayed with a different solid color and a different pattern. The nonhomogenous patterns in the image slices in the stack are associated with those fiber bundles which cut the image plane at an pffiffiffiffiffiffiffiffiffiffiffiffiffi angle measuring cos − 1 ð2 = 3Þ or more. 2.5.3. Grouping: SFM Voxels having nonzero labels are grouped and are given distinct colors resembling the standard color-coded FA maps. For grouping, a cube of size 3×3×3 is considered around each labeled voxel P (having nonzero label). The closeness of orientation (l, m, n) of the e1 of the central voxel P to the e1 of neighboring voxel (with the label the same as that of P) is determined and P is assigned its color accordingly. The color scheme is red (l=±1, m=0, n=0), green (l=0, m=±1, n=0), blue (l=0, m=0, n=±1), yellow (l=±1, m =±1, n=0), cyan (l=0, m=±1, n=±1), magenta (l=±1,m=0, n=±1) and white (l=±1, m=±1, n=±1). If no neighboring voxel has the same label as the central voxel, the voxel is given the mean color (128, 128, 128). Such grouping produces a color-coded map resembling the FA color map with the advantage that the structured white matter in this map has clear boundaries. This map is named as the SFM map. The proximity between the central voxel and its neighboring voxels is determined as follows: I. Consider a voxel P0 (i, j, k) having label (P0)N0. II. Search for a voxel P1 in the 26-neighborhood of the voxel P0 such that P1 has its label the same as that of P0 labelðP1 Þ = labelðP0 Þ

ð11Þ

III. Then, the voxel P0 is assigned a color in accordance with the direction from the central voxel P 0 to neighboring voxel P1 (i.e., according to the position of P1 with respect to P0; (l, m, n) is the difference of the coordinates of P1 and P0).

IV. Repeat I to III for each voxel with nonzero labels in the entire dataset. Fig. 2II shows an SFM map without any FA threshold and an SFM color map with FA N0.15 correspondingly. Fig. 3I–III shows different slices of a normal brain in the axial, coronal and sagittal sections under the FA color map and the SFM map, respectively. 2.5.4. Fiber tracking Three-dimensional tract reconstruction was done by using the FACT method, which involves straightforward linear line propagation in the direction of the principal eigenvector [9,25]. An anisotropy threshold of FA N0.15 is used to restrict the fiber to only white matter. The angle of progression is also restricted; the track is terminated upon obtaining an inner product of less than 0.5 between consecutive eigenvectors. The fiber tracking reported in this work uses ROIs to generate fibers from there as compared with the conventional approach that uses ROIs to retain fibers of interest from the entire fiber mass generated by offering each voxel of the volume stack as a seed voxel to a fiber generation routine. The major point of our approach is indeed the choice of ROIs by utilizing the PEVFS. 2.5.5. Cleaning For more accurate estimations in DTI quantifications along a fiber bundle, unwanted fibers can be deleted using the cleaning procedure developed in this study which is described below. These unwanted fibers may manifest themselves for a variety of reasons such as noise and image imperfections that are present in DTI images and limited resolution of voxels that results in branching, kissing, crossing, etc. To delete a fiber, we try to select it uniquely on its projection on any of the three standard orthogonal planes (axial, sagittal and coronal). Reconstructed white matter tracts can be color-coded according to the principal eigenvector direction at a voxel contained in the tract to aid in the identification of the outlier fibers in any of the planes. The outlier fiber is chosen by using the Point ROI tool of the ROI toolbox of our in-house-developed software. After the selection has been made, a kernel (usually a 3×3 array of all 1's) is imposed on the point and any fiber which passes through a voxel which has two of its coordinates the same as those of the point (or the 2-D region obtained after passing the point through the kernel) is deleted. The seed voxels of the deleted fibers are highlighted in gray color (Fig. 4I (c–e)). They can be used to regenerate and/or save the deleted fibers if required. 2.5.6. Refinement Since the number of generated fibers using pixels from the ROI does not exceed the number of pixels in the ROI, it is clear that the remaining fiber mass would have a wiry or a thinner appearance as compared with the fiber mass passing through the ROI generated by the conventional

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

approach. To compensate for the same and to make the reconstructions time efficient, we use fiber bundle similarity measure to generate similar fibers in a spatial neighborhood

1095

of the seed voxels and continue the same by recursively generating similar fibers from the spatial neighborhoods of the seed voxels considered previously.

Fig. 3. Panels (I)–(III) show axial, coronal and sagittal slices, respectively, of the same section in FA color (upper row) and SFM (bottom row) map. Panel (IV) shows the trajectories intersecting two ROIs [e.g., splenium (A1, A2)] as obtained in the classical approach [29]. Panel (V) shows the splenium (Spl) and genu (Ge) fibers obtained using the same slices as used in Panels (IV, A) and (IV, B), respectively, in the SFM map with the same ROIs automatically drawn this time. [Note the stray fiber (indicated by a white arrow) present in Panel IV (A1, A2) which is absent in Panel V (A). Same is the case in Panel IV (B1, B2).]

1096

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

The similarity measure S(A, B) between a pair (A, B) of fiber pathways is defined [26] as S ðA; BÞ = Rcs e − D=C

medial lemniscus (ML), (16) superior cerebellar peduncle (SCP), (17) middle cerebellar peduncle (MCP) and (18) inferior cerebellar peduncle (ICP).

ð12Þ

where Rcs is the corresponding segment ratio, D is the mean Euclidean distance and C regulates a trade-off between D and Rcs; if LA′B′ is the length of the corresponding segment and LA and LB refer to the lengths of fibers A and B, respectively, the corresponding segment ratio Rcs is defined as the ratio of the length of the corresponding segment to the overall length of the pair of fibers

3.2. Comparison with conventional multiple ROI approach Fibers of genu, splenium, UNC, CG, SLF, ATR, IFO and ILF (Fig. 3IV) were also reconstructed with the conventional multiple ROI methods as suggested in the literature [1,27–29] to compare our method with the conventional multiple ROI methods. 3.3. A comparison model

LA VB V Rcs = LA + LB − LA VB V

ð13Þ

Since for the purposes of cleaning and refinement we also generate the fibers from seed points other than the ones from the initial planar ROI, we have modified the procedure for computing S(A, B) by choosing P0 and Q0 to be the mid-points of the respective fibers. This yields a better score for similarity as the seed points for the fibers may technically be at topographically very different parts of the brain and the distance between them should not play a major role in whether two fibers are similar or not. Let A={P−m,…, P−1, P0, P1, …, Pn} and B={Q−s, …, Q−1, Q0, Q1,…,Qr}, where m=n or n+1, s=r or r+1, and P0 and Q0 are the corresponding mid points of the paths A and B, and m, n, s and r are positive integers. Let u=min {n, r} and v=min {m, s}. Then A′={P−v, P−v+1… P−1, P0, P1… Pu} and B′={Q−v, Q−v+1… Q−1, Q0, Q1…Q u} are said to be corresponding segments. The mean Euclidean distance is the average of the distances P0Q0, P−1Q−1, …, P−vQ−v, P1Q1, P2Q2,…, PuQu. We have thus found the average distance between the two fibers. After calculating the distance, the similarity S(A, B) is calculated as suggested in Eq. (12). The value of similarity between two fibers remains between 0 and 1, and, hence, thresholds can be set by giving a percentage value, e.g., 80%, 90%, etc. 3. Results 3.1. Reproducible WM quantitative tractography The proposed methodology was used for a reproducible, reliable and accurate tractography of all the 18 standard fiber bundles (Fig. 4) of interest in the clinical work and studied in the literature dealing with the identification and contextual anatomical studies of brain fibers: (1) corpus callosum (CC), (2) tapetum (TP), (3) inferior longitudinal fasciculus (ILF), (4) uncinate fasciculus (UNC), (5) inferior fronto-occipital fasciculus (IFO), (6) optic pathways (OP), (7) superior longitudinal fasciculus (SLF), (8) arcuate fasciculus (AF), (9) fornix (FX), (10) cingulum (CG), (11) anterior thalamic radiation (ATR), (12) superior thalamic radiation (STR), (13) posterior thalamic radiation (PTR), (14) corticospinal/corticopontine tract (CST/CPT), (15)

The determination of an ideal true ROI being a major challenge is an issue that is constantly under research in the literature. This model contributes to the same. A manual ROI was drawn on a mid-sagittal color-coded FA image and generated fibers (fiber of interest+stray fibers). The area under ROI was termed as manual ROI. After generation of fiber from manual ROI, stray fibers were deleted. The area in manual ROI which gives stray fibers is termed as stray ROI. A true ROI (manual/automatic) is obtained after the removal of unwanted seed voxels. An ideal true ROI model can be developed in two stages, i.e., true basic ROI and leftover ROI. True basic ROI consists of three types of pixels: (1) those pixels which are in true manual ROI but not in true automatic ROI; (2) those pixels which are in true automatic ROI but not in true manual ROI; (3) those pixels which are common to both true manual ROI and true automatic ROI. The leftover ROI consists of those pixels which are not in either of true manual ROI or of true automatic ROI and is obtained by applying the similarity measure between the true basic fibers and fibers from the surroundings of the true basic ROI (Fig. 5). The similarity threshold is fixed at 75%. Manual ROI−Stray ROI=True Manual ROI Automatic (SFM) ROI−Stray ROI=True (SFM) Automatic ROI True Automatic (SFM) ROI+True Manual ROI=True Basic ROI True Basic ROI+Leftover ROI=Ideal True ROI Using one mid-sagittal section of the CC of each of 20 healthy brains, we quantitatively compared the stray fibers that were obtained using both the FA color map and the SFM map (Table 2). We found that our method (SFM map) has a far lower percentage error as compared to the standard method using manual ROI on FA color map. This further implies that the seed pixels in an ROI that contribute to the generation of the associated fiber mass are predominantly stable voxels. Also, when the voxels on the fibers so generated were counted, it was found that the fibers generated by nonstable seed pixels were relatively much shorter as compared with those generated by stable seed pixels, thereby suggesting that the nonstable seed pixels contribute only

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

1097

Fig. 4. Signature segments and associated fibers for the 18 major fiber bundles. (The signature segment consists of a white boundary and its interior.) (I) For each fiber bundle, Panel (a) shows the section of SFM map from which the ROI (indicated by a white arrow) is chosen by a single click of the mouse; Panels (b) and (d) are segment magnified, (d) shows the seed pixels of the stray fibers colored pink (indicated by a yellow arrow); Panel (c) shows the fibers reconstructed from a segment with the stray fibers having been removed in Panel (e). (II) Signature segments and associated fibers for (A) CST, ML; (B) OP; (C) FX, ALF; (D) ATR, IFO; (E) SCP, ICP, MCP; (F) SLF, ILF; and (G) CG are shown. (III) The final reconstructed 18 fibers, after refinement, are shown.

1098

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

Fig. 5. The manual, SFM, true manual, true SFM, true basic and ideal true ROIs and the generated callosal fibers.

peripherally to the fiber mass. This strongly justifies our choice of generating fiber mass by using PEVFS. 4. Discussion This, to the best of our knowledge, is the first report using an entirely DTI-based automated and validated approach to segment WM structures using the property of stability of a

WM voxel. For the first time, the concept of a stable voxel is introduced and the stability property of diffusion ellipsoid is explored and reported. The classical DTT and PEVFS DTT both utilize the principal diffusion direction for tracking WM. However, their approaches for identifying WM tracts differ. DTT traces eigenvectors to reveal the most likely diffusion pathways taken by the water molecules within a WM tract. Instead of

Table 2 A comparison model for ideal true ROI placed on the corpus callosum No. of pixels in ROI C1 SFM C2 Manual C3 True SFM C4 True manual C5 True basic C6 Ideal true C7—%Error (True SFM) C8—%Error (True manual) 100×(C6−C3)/C6 100×(C6−C4)/C6 Brain 1 Brain 2 Brain 3 Brain 4 Brain 5 Brain 6 Brain 7 Brain 8 Brain 9 Brain 10 Brain 11 Brain 12 Brain 13 Brain 14 Brain 15 Brain 16 Brain 17 Brain 18 Brain 19 Brain 20

591 613 626 743 686 739 816 882 754 699 577 542 465 579 695 560 791 679 704 697

558 637 735 762 711 764 800 875 745 712 618 555 483 587 753 553 750 680 695 738

568 567 597 721 641 707 703 744 666 614 571 519 388 485 680 526 687 630 616 631

533 543 578 710 621 692 693 726 658 625 581 521 378 491 679 486 614 619 608 612

576 578 621 732 657 787 729 780 717 625 612 553 418 518 682 539 712 655 671 682

584 613 680 740 665 830 823 835 815 725 637 563 469 575 705 559 740 679 700 718

2.74 7.50 12.21 2.57 3.61 14.82 14.58 10.90 18.28 15.31 10.36 7.82 17.27 15.65 3.55 5.90 7.16 7.22 12.00 12.12

8.73 11.42 15.00 4.05 6.62 16.63 15.80 13.05 19.26 13.79 8.79 7.46 19.40 14.61 3.69 13.06 17.03 8.84 13.14 14.76

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

tracking the possible pathways, the PEVFS method correlates voxels with the stability property of their PEVs to directly group the voxels of the entire WM tract. The most distinguishing feature between classical DTT and PEVFS is that, while the former yields possible threads of water diffusion, PEVFS directly groups voxels corresponding to the WM tract to produce instant WM segmentation. PEVFS scheme uses fiber tractography for segmentation, but the threshold values (such as FA threshold and the inner product threshold) are not being used during segmentation because it might lead to a loss of data and flawed tracking afterwards. Also, because FA is a quantity for analysis too, using it for segmentation could, in principle, lead to systematic error. PEVFS does not use anisotropy as much, as it relies on PEVs instead. Hence, potentially, it may not only successfully differentiate between WM, GM and CSF but also allow finer classifications by being able to differentiate between adjacent fibers with different orientations. PEVFS provides a SFM map consisting of only true colors which are much less in number than those of the FA color map set (2563?23) and hence makes very small structures clearly visible. Keeping the coloring scheme of SFM as that of the FA color map gives it in an obvious extra advantage over the FA color map. Through its precise description of the contours of a WM structure, it allows us to compare its shape across a selected population, like in the study of the corpus callosum's sexual dimorphism [30,31]. In the classical ROI approach, investigators draw polygons using a mouse over one of the many potential 2-D MR images. However, in SFM mapping the investigator need only click the mouse once on any pixel of the colored segment of interest. Also, implementation of PEVFS permits a local comparison between the voxels involved. When a fiber tract takes turns, bends and changes direction, this local comparison is important, since the directional properties of diffusion will not be similar in different locations along the fiber tract. PEVFS tries to tackle the problems associated with manual ROI selection on all fronts simultaneously, by reducing the time for ROI specification, improving its accuracy and reproducibility by lowering inter- and intrarater variations to effectively zero. In the classical approach, editing of the reconstruction results is done [29] using AND, NOT and CUT operation which require the specification of another ROI whose placement is generally not easily inferable and increases the probability of error. However, in our procedure we do not need all these operations. The cleaning procedure devised is user friendly and could be performed in an automatic mode using the lack of similarity. The current approach of refinement of tracking is more time efficient as compared with the brute-force entire fiber mass approach. It allows for branching as there is no restriction on fibers generated from the new seed voxels being attached to the original bundle. It can include fibers which do not penetrate the marked ROI but belong to the same fiber bundle. The proposed protocol has been employed for the major WM fiber bundles. Note that the possible reproducibility of the conventional multiple

1099

ROI method is basically due to the second ROI which does the essential cleaning. However, there is no fixed methodology for the second ROI whose choice depends upon the fiber bundle of interest. For example, genu, splenium and UNC fiber bundles make the second ROI selection in the same coronal slice; ILF, IFO, CG, SLF and ATR fiber bundles make the second ROI selection using altogether different slices (Fig. 3IV). A number of fiber bundles (CC, TP, MCP, ICP, SCP, ML, PTR, STR and FX) (Fig. 4) for which the classical two-ROI approach is said to be unsuitable to reconstruct them reproducibly get successfully reconstructed reproducibly using a single segmental ROI with our approach. Automatic segmentation of the CC has been a difficult problem to date. Wakana et al. [29] go on to say that it is difficult to devise a protocol for the reconstruction of CC fibers reproducibly. The present approach is fully automated and validated approach to segment the CC on a gamut of healthy subjects [32]. Fig. 4 schematically shows the sagittal sections of the refined fibers collectively, i.e., of all the types that we have discussed. To summarize, the PEVFS segmental method is essentially user independent, robust, efficient and reproducible compared with the classical ROI methods which have resulted in its use in clinical studies [32–36].

References [1] Mori S, Wakana S, Nagae-Poetscher LM, Van Zijl PC. MRI atlas of human white matter. Amsterdam: Elsevier; 2005. [2] Basser PJ, Pajevic S, Pierpaoli C, Duda J, Aldroubi A. In vivo fiber tractography using DT-MRI data. Magn Reson Med 2000;44:625–32. [3] Lazar M, Weinstein DM, Tsuruda JS, Hasan KM, Arfanakis K, Meyer E, et al. White matter tractography using diffusion tensor deflection. Hum Brain Mapp 2003;18:306–21. [4] Parker GJ, Stephan KE, Barker GJ, Rowe JB, MacManus DG, Wheeler-Kingshott CA, et al. Initial demonstration of in vivo tracing of axonal projections in the macaque brain and comparison with the human brain using diffusion tensor imaging and fast marching tractography. Neuroimage 2002;15:797–809. [5] Poupon C, Clark CA, Frouin V, Régis J, Bloch I, Le Bihan D, et al. Regularization of diffusion-based direction maps for the tracking of brain white matter fascicles. Neuroimage 2000;12:184–95. [6] Conturo TE, Lori NF, Cull TS, Akbudak E, Snyder AZ, Shimony JS, et al. Tracking neuronal fiber pathways in the living human brain. Proc Natl Acad Sci U S A 1999;96:10422–7. [7] Jones DK, Simmons A, Williams SC, Horsfield MA. Non-invasive assessment of axonal fiber connectivity in the human brain via diffusion tensor MRI. Magn Reson Med 1999;42:37–41. [8] Mori S, Crain BJ, Chacko VP, van Zijl PCM. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann Neurol 1999;45(2):265–9. [9] Poupon C, Mangin JF, Frouin V, Regis J, Poupon F, Pachot-Clouard M, et al. Regularization of MR diffusion tensor maps for tracking brain white matter bundles. New York: Springer-Verlag; 1998. p. 489–98. [10] Huang H, Zhang J, Van Zijl PC, Mori S. Analysis of noise effects on DTI-based tractography using the brute-force and multi-ROI approach. Magn Reson Med 2004;52:559–65. [11] Pierpaoli C, Barnett A, Pajevic S, Chen R, Penix LR, Virta A, et al. Water diffusion changes in Wallerian degeneration and their dependence on white matter architecture. Neuroimage 2001;13:1174–85.

1100

R.K.S. Rathore et al. / Magnetic Resonance Imaging 29 (2011) 1088–1100

[12] Wiegell MR, Larsson HBW, Wedeen VJ. Fiber crossing in human brain depicted with diffusion tensor MR imaging. Radiology 2000;217 (3):897–903. [13] Kanaan RA, Shergill SS, Barker GJ, Catani M, Ng VW, Howard R, et al. Tract-specific anisotropy measurements in diffusion tensor imaging. Psychiatry Res 2006;146(1):73–82. [14] Jones DK, Symms MR, Cercignani M, Howard RJ. The effect of filter size on VBM analyses of DT-MRI data. Neuroimage 2005;26:546–54. [15] Jonassona L, Bressona X, Hagmann P, Cuisenairea O, Meulib R, Thiran JP. White matter fiber tract segmentation in DT-MRI using geometric flows. Med Image Anal 2005;9:223–36. [16] Wang Z, Vemuri B. Tensor field segmentation using region based active contour model. ECCV 2004;(4):304–15. [17] Rousson M, Lenglet C, Deriche R. Level set and region based surface propagation for diffusion tensor MRI segmentation. ECCV Workshops CVAMIA and MMBIA 2004:123–34. [18] Feddern C, Weickert J, Burgeth B. Level-set methods for tensor valued images. Proceedings of the 2nd IEEE Workshop on Variational, Geometric and Level Set Methods in Computer Vision; 2003. p. 65–72. [19] Paragios N, Deriche R. Coupled geodesic active regions for image segmentation. Technical report. Sophia Antipolis, Cedex, France: INRIA; 1999. p. 3783. [20] Lin C, Sun S, Hong C, Chang C. Unsupervised identification of white matter tracts in a mouse brain using a directional correlation-based region growing (DCRG) algorithm. Neuroimage 2005;28(2):380–8. [21] Niogi SN, Mukherjee P, McCandliss BD. Diffusion tensor imaging segmentation of white matter structures using a Reproducible Objective Quantification Scheme (ROQS). Neuroimage 2007;35: 166–74. [22] Woods RP, Grafton ST, Holmes CJ, Cherry SR, Mazziotta JC. Automated image registration: I General methods and intrasubject, intramodality validation. J Comput Assist Tomogr 1998;22:139–52. [23] Agarwal S, Verma S, Trivedi R, Verma M, Sahoo P, Gupta RK, et al. An interface for DTI tractography. Proceedings of the 17th Annual Meeting of ISMRM, Hawaii, USA; 2009. (abstract 3540). [24] Sun SW, Song SK, Hong CY, Chu WC, Chang C. Directional correlation characterization and classification of white matter tracts. Magn Reson Med 2003;49:271–5. [25] Rathore RKS, Bayu G, Gupta RK, Rathore DK, Purwar A, Saksena S, et al. A comparison of fiber tracking by different numerical integration

[26]

[27] [28] [29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

methods. Proceedings of 14th Annual Meeting of ISMRM, Seattle, USA; 2006. (abstract 528). Zhaohua D, Gore JC, Anderson AW. Classification and quantification of neuronal fiber pathways using diffusion tensor MRI. Magn Reson Med 2003;49:716–21. Catani M, Jones DK, Donato R, Ffytche DH. Occipito-temporal connections in the human brain. Brain 2003;126:2093–107. Concha L, Gross DW, Beaulieu C. Diffusion tensor tractography of the limbic system. AJNR 2005;26:2267–74. Wakana S, Caprihan A, Panzenboeck MM, Fallon JH, Perry M, Gollub RL, et al. Reproducibility of quantitative tractography methods applied to cerebral white matter. Neuroimage 2007;36:630–44. Clarke S, Kraftsik R, van der Loos H, Innocenti G. Forms and measures of adult and developing human corpus callosum: is there sexual dimorphism. J Comp Neurol 1989;280:213–30. Pettey D, Gee J. Sexual dimorphism in the corpus callosum: a characteristic of local size variations and a classification driven approach to morphometry. Neuroimage 2002;17:1504–11. Trivedi R, Agarwal S, Rathore RKS, Saksena S, Tripathi RP, Malik GK, et al. Understanding development and lateralization of major cerebral fiber bundles in pediatric population through quantitative diffusion tensor tractography (DTT). Pediatric Res 2009;66(6): 636–41. Agarwal S, Trivedi R, Gupta RK, Rathore RKS. A principal eigenvector based segmental approach for reproducible white matter quantitative tractography. Proceedings of Annual Meeting of ISMRMESMRMB, Stockholm, Sweden; 2010. (abstract 4011). Kumar M, Pal D, Rathore RKS, Ojha BK, Chandra A, Kumar R, et al. Quantification of DTT metrics in various fiber bundle in patients with frontal lobe injury and its correlation with Neuropsychological tests. Proceedings of Annual Meeting of ISMRM-ESMRMB, Stockholm, Sweden; 2010. (abstract 2118). Trivedi R, Anuradha H, Agarwal A, Rathore RKS, Prasad KN, Tripathi RP, et al. Correlation of quantitative diffusion tensor tractography with clinical grades of subacute sclerosing panencephalitis. Proceedings of Annual Meeting of ISMRM-ESMRMB, Stockholm, Sweden; 2010. (abstract 4347). Kumar R, Saksena S, Hussain M, Srivastav A, Rathore RKS, Agarwal S, et al. Prospective longitudinal study of patients who had sustained moderate TBI and a comparison group of age/sex matched healthy controls. J Head Trauma Rehab 2010;25(1):31–42.