GTRACT V4

GTRACT Overview

This is a manual for the GTRACT program. This program suite now incorporates a full DTI data processing pipeline. In this version of the tools, all of the functional units exist as command line programs that integrate with Slicer3 using the Execution model framework. The standard format for the DTI data in this version is the NRRD format. Anatomical images and label maps can be in any format supported by ITK. The program support several version of fiber tracking including:

1. Free tracking - Basic streamlines with only a starting region
2. Streamlines - Streamlines tracking between two regions of interest
3. Graph Search - Modified streamlines that tries multiple solutions in low anisotropy regions
4. Guided Tracking - Uses an initial guess to the fiber track to guide the tracking
5. Fast Marching - Method based on the method proposed by Geoffrey Parker

Funding for this version of the GTRACT program was provided by NIH/NINDS R01NS050568-01A2S1. The software is hosted by the Neuroimaging Informatics Tools and Resources Clearinghouse (NITRC).

Installing GTRACT

The easiest way to install GTRACT is to obtain pre-compiled binaries from the NITRC website. You should be able to download and install these on most platforms. Currently this is just an archive file for the various platforms that can be unpacked wherever the user has write access.

Installing Into Slicer3

The command line modules do support the Slicer3 execution model and could be placed within the Slicer3 build and will be become active when the user starts Slicer3. The modules are installed under the BRAINS->GTRACT menu in the Slicer3 Module tree. More details will be forthcoming on this. Copyright information for GTRACT can be obtained here.

Building GTRACT

GTRACT depends on only ITK and VTK. It also requires CMake for building. Most recent versions of VTK and ITK should work. Presently GTRACT requires using CMake version 2.8 or higher for the builds.

Conversion from DICOM to NRRD Format

The program for the conversion of DICOM images into Gordan Kindlmann's NRRD format is supported by the Slicer3 DicomToNrrd. This is included within the GTRACT package for convenience. This program should handle most images from GE and Siemens scanners including Mosaic format images. For Siemens scanners, software version above VB13 should be supported. We would like to thank Xiaodong Tao for his great work on this tool. This program should read the diffusion directions from the DICOM header.

DicomToNrrdConverter --inputDicomDirectory DICOM_directory_name --outputVolume NRRD_file_name.nrrd

The output file will contain the Diffusion weighted images in a vector format. The tools within the GTRACT analysis pipeline now expect the images in NRRD format.

GTRACT Programs

The following are links to the full documentation of GTRACT's programs. The same documentation can be accessed on the command line using the following:

${GTRACT_prog_name} --help

• extractNrrdVectorIndex Documentation
• gtractAnisotropyMap Documentation
• gtractAverageBvalues Documentation
• gtractClipAnisotropy Documentation
• gtractCoRegAnatomy Documentation
• gtractConcatDwi Documentation
• gtractCopyImageOrientation Documentation
• gtractCoregBvalues Documentation
• gtractCreateGuideFiber Documentation
• gtractFiberTracking Documentation
• gtractResampleAnisotropy Documentation
• gtractResampleB0 Documentation
• gtractTensor Documentation
• gtractInvertRigidTransform Documentation
• gtractInvertBSplineTransform Documentation
• gtractResampleCodeImage Documentation
• gtractFastMarchingTracking Documentation
• gtractCostFastMarching Documentation

DTI Image Co-Registration

Once the images are in NRRD format the next step in the analysis pipeline is the elimination of motion and eddy current artifacts between the diffusion weighted images. The program for this step of the analysis is gtractCoregBvalues. This program will support registration between two diffusion series in the case where you would like to average this data across runs. It also would allow for co-registration with a T2 weighted image or field map in the same plane as the DTI data.

A mutual information metric cost function is used for the registration because of the differences in the signal intensity as a result of the diffusion gradients. The fixed image for the registration should be a b0 image. The user has the option to perform a rigid or a full affine registration. The full affine allows the registration procedure to correct for eddy current distortions that may exist in the data. This is specified with the eddyCurrentCorrection option. In general the standard parameters work fairly well for a rigid registration, and you probably just need to specify the images for registration. When the eddyCurrent Correction is specified, then we have had to adjust the relaxationFactor and maximumStepSize parameters. We typically use values of around 0.25 for the relaxation factor and 0.1 for the maximum step size when this option is turned on.

% gtractCoregBvalues --fixedImageVolume baseImage.nhdr --movingImageVolume baseImage.nhdr --fixedImageIndex 0 \
 --outputVolume outputImage.nhdr

Tensor Image Calculation

This step takes the diffusion tensor image data and generates a tensor representation of the data based on the signal intensity decay, b values applied and the diffusion directions. The apparent diffusion (ADC) for a given orientation is computed on a pixel-by-pixel basis by fitting the image data (i.e. voxel intensities) to the Stejskal-Tanner equation:

  Sb = S0 exp-b ADC

If at least six orientations are used, then the diffusion tensor can be computed. The measured ADC values is related to the diffusion tensor elements via the following equation:

ADC = M D

where ADC is a vector of the measured apparent diffusion coefficients, M is the matrix describing the diffusion directions applied, and D is a vector of the diffusion tensor elements. This can be solved using standard matrix operations

D = M-1 ADC

In this step we are solving for the diffusion tensor which contains size elements (D_xx, D_yy, D_zz, D_xy, D_xz, D_yz). These can then be used to define the diffusion tensor

|Dxx Dxy Dxz|
|Dyx Dyy Dyz| = D
|Dzx Dzy Dzz|

This is a symmetric matrix therefore D_xy=D_yx, D_xz=D_zx, and D_yz=D_zy.

The program gtractTensor is responsible for the estimation of the Tensor. The program uses itk::DiffusionTensor3DReconstructionImageFilter for this step. There are several options for the analysis that can be tuned by the user.

• Background threshold - This sets a threshold used on the b=0 image to remove background voxels from the processing. Typically a value around 100 works well for Siemens DTI images. We have found that 500 works better for GE datasets. Check your data particularly in the globus pallidus to make sure that brain tissue is not being eliminated with this threshold.

• Median filter - If the filter is enabled, a median filter will be applied to the data. The x, y, and z values specify the size of the median filter in voxels.

• Isotropic Resampling - If this option is enabled, the data will be resample to isotropic resolution specified. This is recommended if fiber tracking is to be performed.

Example:

% gtractTensor --inputVolume Series0003_coreg.nhdr --outputVolume Series0003_Tensor.nhdr \
 --useMedianFilter --medianFilterSize 1,1,1 --backgroundSuppressingThreshold

--Admin 10:44, 2 July 2009 (CDT)

NOTE: We have recently made changes to this program top make it more flexible. There is now a --ignoreIndex option that allows the user to remove specific gradient directions if artifacts exist in a particular gradient direction. A DTI volume can now have multiple B0 and gradient direction images. These are now handled properly without requiring averaging of the diffusion weighted images. Averaging of the B0 image is done internally in this program.

Compute Anisotropy Images

Anisotropy mapping can be used to generate the following diffusion tensor indices:

• Fractional Anisotropy (FA)
• Relative Anisotropy (RA)
• Volume Ratio (VR)
• Lattice Index (LI)
• Coherence Index (CI)
• Mean Diffusivity (ADC)
• Axial Diffusivity (AD)
• Radial Diffusivity (RD)

Anisotropy images are used for fiber tracking, but the anisotropy scalars are not defined along the path. Instead the Tensor representation is included as point data allowing all of these metric to be computed using only the fiber tract point data. The images can be saved in any ITK supported format, but it is suggested that you use an image format that supports the definition of the image origin. This includes NRRD, NifTI, and Meta formats. These images can also be used for scalar analysis including regional anisotropy measures or VBM style analysis.

The command line parameters are defined here gtractAnisotropyMap.

Example:

% gtractAnisotropyMap --inputVolume Series0003_Tensor.nhdr \
  --anisotropyType FA --outputVolume Series0003_FA.nhdr

Alignment with Anatomical Image

The two registration methods supported for alignment with Anatomical images, Rigid and B-Spline. The rigid registration performs a rigid body registration with the anatomical images and should be done as well to initialize the B-Spline transform. The B-SPline transform is the deformable transform, where the user can control the amount of deformation based on the number of control points as well as the maximum distance that these points can move. This allows for some susceptibility related distortions to be removed from the diffusion weighted images. It is recommended that for image co-registration with the B-Spline Transform that the skull stripped (i.e. image containing only brain with skull removed) image be used.

Copy Reference Image Orientation

Currently, the registration process requires that the diffusion weighted images and the anatomical images have the same image orientation (i.e. Axial, Coronal, Sagittal). It is suggested that you copy the image orientation from the diffusion weighted images and apply this to the anatomical image. This image can be subsequently removed after the registration step is complete. We anticipate that this limitation will be removed in future versions of the registration programs. The program gtractCopyImageOrientation is used for this step. The complete command line arguments for gtractCopyImageOrientation can be generated using the '--help' flag.

% gtractCopyImageOrientation --inputVolume T1_ACPC_clip.nii.gz \
     --inputReferenceVolume  DWI_coreg.nhdr --outputVolume T1_ACPC_Axial.nhdr

Anatomical Rigid

The rigid registration uses a mutual information metric to align the b0 image with the anatomical image. A reasonable set of parameters is specified as the default set of parameters. These were initially tested on a data from both a Siemens 1.5T and 3.0T scanners. Here is the complete set of parameters for gtractCoRegAnatomy.

Example

% gtractCoRegAnatomy --inputVolume DWI_coreg.nhdr --transformType Rigid \
     --vectorIndex 0 --inputAnatomicalVolume T1_ACPC_Axial.nhdr \
     --outputTransformName ACPC_rigid.mat --numberOfIterations 1000

It is recommended that the rigid registration be run before the B-Spline nonlinear registration is run.

T1 Volume

T1 Weighted with B0 image overlay

Anatomical B-Spline

The B-SPline registration places a low dimensional grid in the image that is deformed. The B-Spline registration works best if the anatomical image has been skull stripped and a rigid registration is used as the staring point for the registration. The user can specify the number of control points to be used in the registration and the maximum distance that each control point should move. In general the amount of motion in the slice selection and read-out directions direction should be kept low. The distortion is in the phase encoding direction in the images. A bulk transform can be applied to the B-Spline transform to initialize the transform. This should be the rigid transform generated from the previous step for the optimal results. The program gtractCoRegAnatomy performs this registration step and a detailed list of parameters can be found here, gtractCoRegAnatomy

Example:

% gtractCoRegAnatomy --inputVolume DWI_coreg.nhdr --transformType Bspline \
    --vectorIndex 0 --inputRigidTransform ACPC_rigid.mat --inputAnatomicalVolume T1_ACPC_Axial.nhdr \
    --outputTransformName ACPC_bspline.mat --numberOfIterations 1000

T1 Weighted Volume with B0 image (Rigid Registration) overlay

T1 Weighted Volume with B0 image (B-Spline Registration) overlay.

Notice: Improved registration especially at the anterior portion of the corpus callosum

Invert Rigid Transform

The program gtractInvertRigidTransform is used to invert rigid transforms. Typically this is used to map images in anatomical space to the diffusion weighted image space. This is useful for transforming labels defined on the anatomical image into the space of the diffusion weighted images to be used for fiber tracking. Complete parameters for gtractInvertRigidTransform are available using the '--help' option.

Example:

% gtractInvertRigidTransform --inputTransform ACPC_rigid.xfm \
                              --outputTransform ACPC_Inverse_Rigid.xfm

Invert B-Spline

Since the B-Spline transform cannot be directly inverted, we have implemented a thin plate spline approximation to its inverse. This can be created using the gtractInvertBspline. The X, Y and Z grids specify the grid utilized to sample the B-Spline transform. The parameters for gtractInvertBSplineTransform can be obtained using the '-help' option. The inverse needs to be computed to map fiber tracts from the DTI space to the anatomical image and to map regions of interest from the anatomical image onto the diffusion weighted images using this non-linear transform that accounts in part for susceptibility artifacts.

Example

% gtractInvertBSplineTransform --inputReferenceVolume T1_ACPC.nhdr --inputTransform ACPC_bspline.xfm \
                             --outputTransform ACPC_Inverse_Bspline.xfm

Resampling Images

Three programs have been written to resample the various types of images that may be encountered in diffusion tensor analysis. All of the programs accept only ITK transforms and in particular Versor Rigid or B-Spline' transforms that are described above.

• gtractResampleCodeImage - Resamples a mask or label - Image pixel type is signed short and nearest neighbor interpolation is used

Example:

% gtractResampleCodeImage --inputTransform ACPC_Inverse_Rigid.xfm \
   --transformType Rigid --inputCodeVolume Thal_axial.nhdr \
   --inputReferenceVolume DWI_coreg.nhdr --outputVolume Thal_DWI.nhdr

• gtractResampleB0 - Resamples the B0 image - Image pixel type is signed short and linear interpolation is used

Example:

% gtractResampleB0 --inputVolume DWI_Coreg.nhdr --inputAnatomicalVolume T1_ACPC.nhdr \
   --inputTransform ACPC_Rigid.txt --vectorIndex 0 \
   --transformType Rigid --outputVolume B0_ACPC_rigid.nii.gz

• gtractResampleAnisotropy - Resamples an anisotropy map - Image pixel type is float and linear interpolation is used

Example:

% gtractResampleAnisotropy --inputAnisotropyVolume FA.nhdr --inputAnatomicalVolume T1_ACPC.nhdr \
   --inputTransform ACPC_Rigid.txt --vectorIndex 0 \
   --transformType Rigid --outputVolume FA_ACPC_rigid.nii.gz

When running these programs the user must also specify a template image that will be used to define the origin, orientation, spacing and size of the resampled image.

You may want to save the aligned B0 image after each of the co-registration steps with the anatomical image to check the registration quality with another tool such as BRAINS or Slicer.

Fiber Tracking

Once an isotropy image, a tensor image, and a spatial object have either been created and loaded fiber tracking can take place. Selecting the Tracking tab will allow the user to specify the type of tracking to be performed. Currently the GTRACT tool supports a number of fiber tracking methods. These include:

1. Free Tracking
2. Streamlines Tracking
3. Graph Search Tracking
4. Guided Tracking
5. Fast Marching Tracking

The output of the fiber tracking is vtkPolyData (i.e. Polylines) that can be loaded into Slicer3 for visualization. The poly data can be saved in either old VTK format files (.vtk) or in the new VTK XML format (.xml). Currently, Slicer3 only accepts the old .vtk file format. The polylines contain point data that defines the Tensor at each point along the fiber tract. This can then be used to be rendered as glyphs in Slicer3 and can be used to define several scalar measures without referencing back to the anisotropy images.

Free Tracking

The Free tracking is a basic streamlines algorithm. This is a direct implementation of the method original proposed by Basser et al. The tracking follows the primary eigenvector. The tracking begins with seed points in the starting region. Only those voxels above the specified anisotropy threshold in the starting region are used as seed points. Tracking terminates either as a result of maximum fiber length, low anisotropy, or large curvature. This is a great way to explore your data, but the tools to appropriately edit these down to a fiber tract of interest do not yet exist. The user should explore using the other approaches described below once they have performed an initial review of their data. A complete set of parameters for gtractFiberTracking can be obtained using the '--help' option.

Example:

% gtractFiberTracking --inputTensorVolume Tensor.nhdr \
   --inputAnisotropyVolume FA.nhdr \
   --trackingMethod Free \
   --inputStartingSeedsLabelMapVolume Thal_DWI.nhdr \
   --startingSeedsLabel 1 --seedThreshold 0.25 \
   --curvatureThreshold 60 --outputTract freeTracking.vtk

Streamline Tracking

The streamlines algorithm is a direct implementation of the method original proposed by Basser et al. The tracking follows the primary eigenvector. The tracking begins with seed points in the starting region. Only those voxels above the specified anisotropy threshold in the starting region are used as seed points. Tracking terminates either by reaching the ending region or reaching some stopping criteria. Stopping criteria are specified using the following parameters:

• Tracking Threshold: Anisotropy values of the next point along the path
• Curvature: Angle of the next tracking direction relative to the current path direction
• Max Length: Length of the track in voxels

Only paths terminating in the ending region are kept in this method. The TEND algorithm proposed by Lazar et al. (Human Brain Mapping 18:306-321, 2003) has been instrumented. This can be enabled using the --useTend option while performing Streamlines tracking. This utilizes the entire diffusion tensor to deflect the incoming vector instead of simply following the primary eigenvector. The TEND parameters are set using the --tendF and --tendG options.

v_out = fe₁ + (1-f) * ((1-g)v_in+gD · v_in)

Where:

• v_out is the outgoing vector direction
• v_in is the incoming vector direction defined in the previous step
• e_in is the primary eigenvector direction
• D is the diffusion Tensor

The complete set of command line arguments for gtractFiberTracking are available with the '--help' option at the command line.

Example:

% gtractFiberTracking --inputTensorVolume Tensor.nhdr \
  --inputAnisotropyVolume FA.nhdr \
  --trackingMethod Streamline \
  --inputStartingSeedsLabelMapVolume startingSeedsLabelMap.nhdr \
  --startingSeedsLabel 1 \
  --inputEndingSeedsLabelMapVolume endingSeedsLabelMap.nhdr \
  --endingSeedsLabel 2 \
  --curvatureThreshold 60 \
  --outputTract streamlineTracking.vtk

Graph Search Tracking

The Graph Search tracking is the first step in the full GTRACT algorithm developed by Cheng et al. (NeuroImage 31(3): 1075-1085, 2006) for finding the tracks in a tensor image. This method was developed to generate fibers in a Tensor representation where crossing fibers occur. The graph search algorithm follows the primary eigenvector in non-ambiguous regions and utilizes branching and a graph search algorithm in ambiguous regions. Ambiguous tracking regions are defined based on two criteria:

• Branching AI Threshold: Anisotropy values below this value and above the tracking threshold
• Curvature Major Eigen: Angles of the primary eigenvector direction and the current tracking direction

In regions that meet this criteria, two or three tracking paths are considered. The first is the standard primary eigenvector direction. The second is the secondary eigenvector direction. This is based on the assumption that these regions may be prolate regions. If the Random Walk option is selected then a third direction is also considered. This direction is defined by a cone pointing from the current position to the centroid of the ending region. The interior angle of the cone is specified by the user with the Branch/Guide Angle parameter. A vector contained inside of the cone is selected at random and used as the third direction. This method can also utilize the TEND option where the primary tracking direction is that specified by the TEND method instead of the primary eigenvector.

The parameter '--maximumBranchPoints' allows the tracking to have this number of branches being considered at a time. If this number of branch points is exceeded at any time, then the algorithm will revert back to a streamline algorithm until the number of branches is reduced. This allows the user to constrain the computational complexity of the algorithm. A complete set of command line arguments for gtractFiberTracking can be obtained using the '--help' option from the command line.

Example:

% gtractFiberTracking --inputTensorVolume Tensor.nhdr \
 --inputAnisotropyVolume FA.nhdr \
 --trackingMethod GraphSearch \
 --inputStartingSeedsLabelMapVolume startingSeedsLabelMap.nhdr \
 --startingSeedsLabel 1 \
 --inputEndingSeedsLabelMapVolume endingSeedsLabelMap.nhdr \
 --endingSeedsLabel 2 \
 --inputTract guideFiber.vtk \
 --curvatureThreshold 60 \
 --branchingThreshold 0.20 \
 --maximumBranchPoints 5 \
 --branchingAngle 45.0 \
 --outputTract graphSearchTracking.vtk

Guided Tracking

The second phase of the GTRACT algorithm is Guided Tracking. This method incorporates anatomical information about the track orientation using an initial guess of the fiber track. In the originally proposed GTRACT method, this would be created from the fibers resulting from the Graph Search tracking. However, in practice this can be created using any method and could be defined manually. To create the guide fiber the program gtractCreateGuideFiber can be used. This program will load a fiber tract that has been generated and create a centerline representation of the fiber tract (i.e. a single fiber). The complete set of parameters for gtractCreateGuideFiber can be obtained using the '--help' option.

Example:

%  gtractCreateGuideFiber --inputFiber graph.vtk --outputFiber guide.vtk --numberOfPoints 100

In this method, the fiber tracking follows the primary eigenvector direction unless it deviates from the guide fiber track by a angle greater than that specified by the '--guidedCurvatureThreshold' parameter. The program gtractFiberTracking uses this algorithm. The user must specify the guide fiber when running this program.

Example:

%  gtractFiberTracking --inputTensorVolume Tensor.nhdr \
 --inputAnisotropyVolume FA.nhdr \
 --inputStartingSeedsLabelMapVolume --startingSeed.nhdr \
 --startingSeedsLabel 1 \
 --inputEndingSeedsLabelMapVolume --endingSeed.nhdr \
 --endingSeedsLable 2 \
 --inputTract guide.vtk \
 --trackingMethod Guided \
 --guidedCurvatureThreshold 30.0 \
 --maximumGuideDistance 12.0 \
 --outputTract guidedTracking.vtk

Fast Marching Tracking

This algorithm implements a fast marching fiber tracking algorithm to identify fiber tracts from a tensor image. This algorithm is based on the work by G. Parker et al. (IEEE Transactions On Medical Imaging, 21(5): 505-512, 2002) An additional feature of this version is the inclusion of anisotropy in the cost function calculation. The user has the ability to weight the contribution of the anisotropy weight.

This algorithm is broken into two parts. The first phase of the algorithm generates a cost image. This starts with a seed region specified by the user. The cost of the path tracking from the seed region to every voxel in the brain is computed. The cost is based on the direction of the front, the primary eigenvector direction, and the anisotropy value. The user can adjust the weight of the anisotropy value using the '--anisotropyWeight' flag. The program gtractCostFastMarching is responsible for the generation of the cost image. This may take several minutes to generate and increases with the size of the seed region. Currently, both a cost and speed image are generated. The speed image is not used in the fiber tracking phase and will become an optional parameter in future releases. Currently, it is a required parameter.

Example:

% gtractCostFastMarching --inputTensorVolume Tensor.nhdr \
    --inputAnisotropyVolume FA.nhdr \
    --inputStartingSeedsLabelMapVolume seeds.nhdr \
    --startingSeedsLabel 1 --outputCostVolume cost.nhdr \
    --outputSpeedVolume speed.nhdr --anisotropyWeight 0.5

Once the cost image is generated, the user can then generate tracts between the any region of the brain and the seed region specified in the previous step. Fiber tracking is fairly quick and uses a gradient descent solution from the defined regions back to the seed points specified in gtractCostFastMarching. The parameters for gtractFastMarchingTracking can be obtained using the '-help' option from the command line. The resulting paths can be constrained using the anisotropy threshold and the number of iterations. This will help to eliminate unrealistic fiber tracts.

Example:

% gtractFastMarchingTracking --inputTensorVolume Tensor.nhdr \
    --inputAnisotropyVolume FA.nhdr \
    --inputStartingSeedsLabelMapVolume seedsEnd.nhdr \
    --startingSeedsLabel 1 --inputCostVolume cost.nhdr \
    --seedThreshold 0.3 --trackingThreshold 0.2 \
    --outputTract fastMarching.vtk

Visualizing Fiber Tracts

The resulting fiber tracts can be loaded into Slicer3 or Paraview. After starting Slicer3, the loading of fiber tracts can be found under the Tractography->DisplayLoadSave module. Selecting this module will bring up the interface to load the fiber tracts. Select the Load Tractography Button and select the fiber tract of interest. The user can also select the type of glyph to be use for display, including, lines, tubes, or ellipsoids. The fiber tracts will be visible in the 3D view and can be overlaid on the image data as well.

File:Slicer3-Display-Tracks.jpg Slicer3 DisplayLoadSave Tab