How to build a tract specific white matter Atlas

Previous steps

What do you have to get installed?

Creating a tract specific white matter atlas is the result of the combination of three main open-source software developed in PICSL. The first main part relies on DTI-TK. The last steps rely on the continuous medial representation model (cm-rep) (see this page for installation and this page for documentation). Another companion tool that will be needed here is convert3D (see here).

Two visualization softwares can help you during the course of this tutorial:

1. ITK-SNAP
2. paraview.

The former is a good segmentation and visualization tools for nifti images, whereas the latter is convenient to deal with .vtk format images.

Building the whole brain white matter Atlas using DTITK

The first task is to build a whole brain white matter atlas (see tutorial page. If you are already familiar with DTI-TK, you can skip the first points of the tutorial and directly go to Registration and Spatial normalization of DTI volumes. At the end of his, you will have created a white matter template (see Figure 1), which is necessary to build the tract specific template.

Figure 1

Main Goal

Now the whole brain DT atlas has been created from your population, the main point is to delineate the tracts of interest and then create a medial (surface-based) representation for each one. As of December 2010, these tracts of interest that have been considered so far are

• the corpus callosum (CC)
• the inferior longitudinal fasciculus (ILF)
• the uncinate fasciculus (UNC)
• the superior longitudinal fasciculus (SLF)
• the cortico-spinal tract (CST)
• the inferior fronto-occipital fasciculus (IFO)

Such an atlas was already created using an aging population (65-83 yo) obtained from the ixi database (see in the example folder). This study was published in \cite{zhang2009}. Another one is in the process of being created. The underlying database is the same but we selected subjects belonging to an adult population from age 18 to 60.

Following the generation of this template, several applications are possible, such as the detection of white matter differences between two populations in these particular tracts (an example is in the aforementioned Zhang paper, where WM changes are identified between ALS and healthy subjects, based on the fractional anisotropy measure) or such as longitudinal studies.

Creating the Surface Representation- 1st part- starting with DTI-TK

What to start with?

The preliminary steps consist of creating an FA images, a white matter mask as well as a whole brain tractography from the DT image.

TVResample -in tensor.nii.gz -xs 224 -ys 224 -zs 144 -xv 1 -yv 1 -zv 1})

TVtool -in tensor.nii.gz -fa. This will output tensor_fa.nii.gz

BinaryThresholdImageFilter tensor_fa.nii.gz wm.nii.gz 0 0.2 0 1

A whole brain tractography will be needed throughout this specific tract atlas construction. To do so, use

SingleTensorFT -in tensor.nii.gz -seed wm.nii.gz -out tensor.vtk

Other factors can be added, such as

• -label (default 1) - to be modified if the label of the seed (i.e., the white matter mask) is not 1
• -minFA (default: 0.15) and `-maxAngle' (default: 45) which are stopping criteria
• -stepSize (default: 0.5 mm)
• -upsample (default 1). This number can be changed when using DTIs with very low resolution. To get denser tracts, the value can be set to 2 (or other puissance of 2). It inserts additional point within a

voxel.

First Step: delineating the fibers corresponding to each tract

If you are not familiar with brain white matter tracts in general, you can gather some information in Mori et al., 2002. This paper has been revealed to be particularly helpful to delineate the aforementioned structures. Other resources can be helpful, such as Wakana et al., 2003 and Mori et al., 2002.

The current aim is to delineate each tract. This might be the most important step because everything else results from this.

1. In ITK-Snap, open the FA image (figure 2)
2. Open the white matter mask as an overlay (figure 2)
3. Use the polygon tool to manually trace the ROIs through which the specific tract you are focusing on is passing
4. Once you are satisfied, go to Segmentation -> Save as Image -> Save your ROI as a *.nii.gz
5. To remember to which tract each ROI corresponds, you can use the label editor and name each ROI, then save the label descriptions -> Segmentation -> Save Label Descriptions
6. You will realize very frequently (ok, all the time!) that you not only have to draw ROIs through which the fibers are supposed to go through, but you also have to draw ROIs to exclude some fibers. An example of inclusion ROIs is given for the corpus callosum in figure 2. An example of exclusion ROI is given in figure 3.

Figure 2: An example with the aging atlas. Top left: FA - Top right: FA with white matter mask Bottom: ROIs delineating the corpus callosum

Figure 3: Exclusion ROIs for the corpus callosum

Second Step: Defining each tract

For each tract, you should have drawn a few ROIs through which you expect the fibers to go (the inclusion ROIs) and a few (the exclusion ROIs) through which you expect your fibers not to go, which are here to get rid of the extra annoying fibers. The main tool that is used to generate each tract is called Tracttool. Tracttool will be used to 1/filter the tracts, i.e., create one .vtk file per fiber tract -2/ convert the *.vtk into a .nii file (grey level image of the tract). Obtaining each fiber tract is an iterative procedure. The whole brain tractography (called tensor.vtk see paragraph first step) is filtered again and again until a reasonable fiber tract is obtained. Here is the way to do it:

TractTool -in tensor.vtk -filter filtered_tract1.vtk -roi ROI_inclusion1.nii.gz -label label_for_roi1

(this step is repeated for all the inclusion ROIs. -label is not required is the value of the label used in ITK-SNAP is 1)

 
TractTool -in filtered_tract1.vtk -filter filtered_tract2.vtk -roi ROI_inclusion2.nii.gz -label
label_for_exclusion_roi -inclusion 0

 TractTool -in filtered_tract2.vtk -nifti output.nii.gz -target target.nii.gz

The output output.nii.gz will have the same dimensions as -target, which can be a volume of your choice.

From there, if you notice that your fiber tract is not correct, you can refine the process by creating more exclusion ROIs (see the First step) and re-creating your output (see Second Step). An example of what you can achieve for the CST is given in the example directory.

Third Step: From the fiber tracts to the medial rep

The model behind this representation relies on Yushkevich et al., 2008, which focuses on white matter tracts that are thin sheet-like structures. After segmentation of these tracts, the cm-rep model is fitted to them that describes the skeleton and the boundary of geometrical objects as parametric digital surfaces with predefined topology. This step relies on the continuous medial representation tool developed by Yushkevich et al. The main documentation and binaries can be found here. It also utilizes the convert3D tools that you can find and download here. Now you have your segmentations ready (they are grey-level images called, let's say, tract_greylevel.nii.gz), a few preprocessing steps are needed before building the cm-rep. These use convert3D. The main goal is to obtain a levelset image from the grey level segmentation. As stated in section levelset n_iter, creating a levelset segmentation requires a speed function and an initialization function. The adapted commands are as follows (the recommended parameters are indicated here)

 
c3d tract_greylevel.nii.gz -threshold 0.5 inf 1.0 0.0 -o tract_binary.nii.gz 

(see note at the bottom of this paragraph)

 
c3d tract_binary.nii.gz -shift -0.1 -o tract_speed.nii.gz

 c3d tract_binary.nii.gz -threshold 0.25 inf -0.5 0.5 -o tract_init.nii.gz

 
c3d tract_speed.nii.gz tract_init.nii.gz -levelset-curvature 0.02 -levelset 100 -o tract_levelset.nii.gz

Congratulations! You have made it through the first part. There is still quite a way to go, but we are almost there.

Note: In some cases, especially for the CST, the filtering step is not enough to delineate the portion of the structure that is sheet-like. Consequently, I strongly suggest to act on the binary structure once it is created in the first step (i.e., step that converts the grey level segmentation into a binary image). At this point, you can open the binary image in ITK-SNAP (opening it as the image and as the segmentation is fine). Go to the 3D window -> click update mesh -> click 3D scalpel tool in the 3D Toolbox -> Change the label number to clear label in the Segmentation options and cut the stem of the CST -> Accept and update mesh. Save the new segmentation.

Figure 4 gives you an example.

Figure 4: An example of initial binary files for the cortico-spinal tract. In the bottom image, the stem was cut. The middle and right column illustrate how these two binaries are transformed into skeletons and boundaries during the next steps (see paragraph Building the initial template -1). Modeling the stem using the cmrep model would not be adapted.

Creating the Surface Representation- 2nd part- the Continuous medial representation}

We are going to switch to and focus on the cmp-rep model properly speaking, and use the continuous medical rep tools only. The principle is to build the boundary and the corresponding skeleton for each tract.

Figure 5: Left: An example of initial boundary for the corpus callosum - Right: the corresponding skeleton

Building the initial template -1

In order to initiate the cmrep model, you need a primary template. The construction of this initial template starts with the generation of a skeleton and boundary from the previous levelset file.

 vtklevelset tract_levelset.nii.gz tract_boundary.vtk 0.0

The last and only parameter is the threshold value. The higher the value, the slimmer the boundary.

  cmrep_vskel -Q /PATH_TO_QVORONOI/qvoronoi -p ''p'' -e ''e'' -q ''q'' tract_boundary.vtk tract_skel.vtk

The main factor here is the pruning factor p. If p=0, there is no pruning at all and the skeleton is very rough. The value of $p=1$ is the point, where you can notice pruning. For thinnest structures, it is advised to choose p=1.3 or 1.4. For the structures that contains less detail, p=1.8 to p=2 might be a good estimate. Modifying $p$ gives reasonable results for most of the tracts, however, some issues might happen when the tracts are at the limit of being sheet-like. At this point, you can also tweak the parameter e that is the minimal number of mesh edges separating two generator points of a VD face for it to be considered (e = 3, otherwise) . The last parameter, q is the number of bins in each dimension (q has a default value, so you don't need to specify it). It defines the speed with which the program runs. A good value is q=20-50. Reducing q reduces the number of vertices in the skeleton. For more information on each of the parameters and additional ones, you can type the command cmrep_vskel on the command line. The results should look like the ones of Figure 5 for the corpus callosum or of Figure 4 bottom line for the corticospinal tract.

Building the initial template -2

The next action is to build the true initial template, the one that will serve as an input in the final optimization. This is done with a matlab function, called tempinit.m. This function calls a few other ones, such as vtk_polydata_read.m, vtk_get_point_data.m as well as a package written by Kilian Weinberger that you can find here. tempinit uses this maximum variable unfolding (MVU) algorithm, the general purpose of which is to reduce dimensionality. This function also depends on a library called Csdp, which you can download here. This library contains an executable (found in /bin/csdp), a copy of which should be put in the user's bin directory of the machine where you are running tempinit. The same applies to the c function called triangle. Once compiled, the executable should similarly go to your local bin directory. For more information about the different parameters of this matlab function, you can have a look at the .m file located in the example directory, where defaults and explanations are given for each of them. This program should output two files: medialtemplate.cmrep and medialtemplate.vtk, which are be the starting point for the very last step.

Building the initial template -3

The very last step of this process uses the function cmrep_fit from the cm-rep package that we previously installed. Its inputs consist of the medialtemlate.cmrep files obtained above, the binary image as well as a parameter file and the output directory. The main command is

cmrep_fit param.txt medialtemplate.cmrep tract_binary.nii.gz cmrep_output

The explanation for most of the parameters can be found on this page. The first part of this parameter file contains the default parameters that control the values that you will see appearing on your terminal after launching this command (see Figure 6 and the param_output.txt in the demo folder provided with this tutorial). A few important facts should be retained. During the different iterations you should control that the TOTAL (see Figure 6) is decreasing. Each of the values VOLOV to BNDTRI contribute to this TOTAL and are driven by the default parameters. So you can adjust them by modifying the term weights of the parameter file.

1. VOLOV translates the overlap between the model that is iteratively made and the binary
2. BNDJAC is driven by DefaultParameters.BoundaryJacobianEnergyTerm.Weight
3. DefaultParameters.BoundaryJacobianEnergyTerm.PenaltyA and
4. DefaultParameters.BoundaryJacobianEnergyTerm.PenaltyB
5. GRAD-R corresponds to DefaultParameters.BoundaryGradRPenaltyTerm.Weight
6. RADIUS to DefaultParameters.RadiusPenaltyTerm.Weight
7. LOOPTN to DefaultParameters.LoopTangentSchemeValidityPenaltyTerm.Weight
8. MEDTRI and BNDTRI express the quality of the triangulation, hence how much your model deform

Figure 6: Example of output after submitting the cmrep_fit command. The total value decreases and each of the terms VOLOV to BNDTRI contribute to this decrease.

The second part of the parameter file controls each step of the run, from the alignment stage (translation and rotation) to the deformation stage (fitting of the model). All the parameters are straightforward to understand and thus to modify. One comment could be done about the blurring parameter, which can be tuned down if there are some finer structures to be preserved in the tract.

The output of this final step will be located in the cmrep_output directory, under three subdirectories, called cmrep, image and mesh. The optimized structure are in cmrep_output/mesh/, and called def3.med.vtk (the medial representation) and def3.bnd.vtk (boundary). An example can be found in Figure 7.

Figure 7: Pipeline illustrated on the CST

Congratulations, you have obtained your first white matter tract medial representation!

Checking your template

In order to check the accuracy of the template, it is worth comparing it to the original binary file, through a Dice overlap for instance. cmrep can be used to generate a binary file from the def3.med.vtk

cmrep_fillmesh def3.med.vtk ref.nii.gz def3.nii.gz 

• ref.nii.gz can be the binary image itself.
• The binary def3.nii.gz obtained from the information of the medial representation is then compared to the initial binary model through c3d -verbose def3.nii.gz tract_binary.nii.gz -overlap 1 (if the foreground of the binary image is 1). This command results in the computation of the Dice overlap. A value of 0.88 is satisfying.