Processing pipeline

To call the FDT GUI, either run Fdt (Fdt_gui on a Mac), or run fsl and press the FDT button.

A typical processing pipeline (and approximate time required for each stage, based on an Intel 2.66GHz processor, and a 60-direction whole-brain dataset of dimensions 128x128x70, at 2 mm isotropic resolution) would consist of:

1. Any study or scanner-specific pre-processing (e.g., conversion from DICOM to NIFTI using for example dcm2niix, removal of images affected by large artifacts). This would be done manually by the user.

2. Check your diffusion data (see below)

3. topup - Susceptibility-induced distortion correction (fieldmap estimation) using topup.

4. bet - Brain extraction using BET on the distortion-corrected b0 (output of topup).

5. eddy - Distortion correction using eddy. eddy simultaneously corrects for eddy currents and subject motion. It can also use a fieldmap estimated by topup and correct for susceptibility-induced distortions. (~hours-days on CPU or less if multi-threaded or few hours on GPU).

After running eddy, all the diffusion-weighted images should be in alignment with each other and be undistorted. We can now either run a microstructural analysis or tractography. For a microstructural analysis the further pipeline would be:

1. dtifit - Fitting of diffusion tensors on eddy-corrected data (~1 minute)

2. TBSS - Comparison of the fractional anisotropy maps (or others) between subjects (see TBSS)

Or run tractography we need to:

1. bedpostx - Fitting of the probabilistic diffusion model on corrected data (~15 hours single-threaded on CPU or less if multi-threaded or ~30 minutes run on GPU)

2. Registration to structural reference image and/or standard space - (3-6 minutes)

3. probtrackx - Probabilistic tractography run on the output of bedpostx (execution time depends very much on what the user wishes to do. Generating a connectivity distribution from a single voxel of interest takes about 1 second)

4. Further post-processing of probtrackx outputs can be carried out if required using the command-line utilities

Diffusion data in FSL

Diffusion data is usually stored in a 4D NIFTI file, where each volume is acquired with a different b-value and/or gradient orientation. To run FDT tools we need:

• Diffusion weighted data (data.nii.gz): A 4D series of data volumes. This will include diffusion-weighted volumes and volume(s) with no diffusion weighting.

• Gradient directions (bvecs): An ASCII text file containing a list of gradient directions applied during diffusion weighted volumes. The order of entries in this file must match the order of volumes in the input data series.

The format is

x_1 x_2 x_3 ... x_n
y_1 y_2 y_3 ... y_n
z_1 z_2 z_3 ... z_n

Vectors are normalised to unit length within the dtifit code. For volumes in which there was no diffusion weighting, the entry should still be present, although the direction of the vector does not matter! For technical details see this FAQ entry.

• b values (bvals): An ASCII text file containing a list of b values applied during each volume acquisition. The b values are assumed to be in s/mm^2 units. The order of entries in this file must match the order of volumes in the input data and entries in the gradient directions text file.

The format is

b_1 b_2 b_3 ... b_n

Check your data

To check whether your bvecs and bvals are correct, it is recommended to run a quick check on the raw data. This can be run after preprocessing or even before if the subject did not move too much:

1. Brain extraction using BET.

2. Run dtifit

3. View the principal eigenvector (V1) to check the vectors are correctly oriented with respect to the anatomy. If there is a problem then the bvecs need to be modified (change signs and/or permute order of components) in order to get the V1 vectors correctly oriented (repeat dtifit after modifying bvecs and view V1 vectors again until they are correct) before proceeding with any further processing.

Distortion Correction

Diffusion weighted images tend to be distorted due to susceptibility-induced distortions and eddy currents in the gradient coils. EDDY offers a state of the art method for correcting eddy current distortions and subject motion. It can also use a fieldmap (estimated by TOPUP) to consider and correct susceptibility-induced distortions. It also allows detection of outliers, data points that have been corrupted by motion-induced signal dropout. Finally, it permits slice-wise outlier detection and replacement, and slice-to-volume registration. We refer the user to the detailed documentation page of eddy for further details.

Note: The old eddy_correct tool has now been superseded by eddy that offers substantially better distortion correction.

DTIFIT

DTIFIT fits a diffusion tensor model at each voxel. You would typically run dtifit on data that has been pre-processed and eddy current corrected. Note that dtifit is not necessary in order to use the probabilistic tractography (which depends on the output of BEDPOSTX, not DTIFIT).

In the FDT GUI, use the top left drop down menu to select DTIFIT.

Input: You can specify an input directory containing all the required files with standardized filenames, or alternatively you can specify input files manually by turning on the specify input files manually switch. If an input directory is specified then all files must be named as shown in parentheses below. If input files are specified manually they can have any filename. Required files are:

• Diffusion weighted data (data): A 4D series of data volumes. This will include diffusion-weighted volumes and volume(s) with no diffusion weighting.

• BET binary brain mask (nodif_brain_mask): A single binarised volume in diffusion space containing ones inside the brain and zeroes outside the brain.

• Output basename: User specifies a basename that will be used to name the outputs of dtifit. If the directory input option is used then the basename will be dti.

• Gradient directions (bvecs): An ASCII text file containing a list of gradient directions applied during diffusion weighted volumes.
• b values (bvals): An ASCII text file containing a list of b values applied during each volume acquisition.

Advanced options In addition to the required input above, the user can choose to apply a weighted least-squares regression instead of the default standard linear regression. The user can also choose to save the tensor elements and/or the sum of squared error. This last output can be useful for detecting artefacts.

Outputs of dtifit

• <basename>_V1 - 1st eigenvector

• <basename>_V2 - 2nd eigenvector

• <basename>_V3 - 3rd eigenvector

• <basename>_L1 - 1st eigenvalue

• <basename>_L2 - 2nd eigenvalue

• <basename>_L3 - 3rd eigenvalue

• <basename>_MD - mean diffusivity

• <basename>_FA - fractional anisotropy (isotropic ~ 0; stick-like ~1)

• <basename>_MO - mode of the anisotropy (oblate ~ -1; isotropic ~ 0; prolate ~ 1)

• <basename>_S0 - raw T2 signal with no diffusion weighting

Optional output

• <basename>_sse - Sum of squared error of diffusion tensor fit.

• <basename>_tensor - tensor as a 4D file in this order: Dxx, Dxy, Dxz, Dyy, Dyz, Dzz

command line utility

dtifit

Compulsory arguments (You MUST set one or more of):
        -k,--data       dti data file
        -o,--out        Output basename
        -m,--mask       Bet binary mask file
        -r,--bvecs      b vectors file
        -b,--bvals      b values file

Optional arguments (You may optionally specify one or more of):
        -V,--verbose    switch on diagnostic messages
        -h,--help       display this message
        --cni   Input confound regressors
        --sse   Output sum of squared errors
        -w,--wls        Fit the tensor with weighted least squares
        --kurt  Output mean kurtosis map (for multi-shell data)
        --kurtdir       Output  parallel/perpendicular kurtosis maps (for multi-shell data)
        --littlebit     Only process small area of brain
        --save_tensor   Save the elements of the tensor
        -z,--zmin       min z
        -Z,--zmax       max z
        -y,--ymin       min y
        -Y,--ymax       max y
        -x,--xmin       min x
        -X,--xmax       max x
        --gradnonlin    Gradient Nonlinearity Tensor file

Two notable options only available from the command line are:

• --kurt: Adds a mean kurtosis to the fit (for multi-shell data), which is output in a <basename>_kurt file. Note that fitting the mean kurtosis is different from the full kurtosis tensor model.

• --gradnonlin: If provided, dtifit will consider the variability in the b-value and gradient orientation due to gradient non-linearities.

BEDPOSTX

BEDPOSTX stands for Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques. The X stands for modelling Crossing Fibres. bedpostx runs Markov Chain Monte Carlo sampling to build up distributions on diffusion parameters at each voxel. It creates all the files necessary for running probabilistic tractography. For an overview of the modelling carried out within bedpostx see this technical report.

bedpostx models crossing fibres within each voxel of the brain. Crucially, bedpostx automatically determines the number of crossing fibres per voxel. For details on the basic crossing fibre model see Behrens et al, NeuroImage 2007. By default, bedpostx runs a slightly more complex model Jbabdi et al. MRM 2012, where the diffusion coefficient is modelled using a Gamma distribution.

bedpostx takes about 15 hours to run but will automatically batch if run on an SGE-capable system.

Note that bedpostx is a wrapper script for a command-line tool called xfibres.

In the FDT GUI, use the top left drop down menu to select BEDPOSTX.

Input directory: Use the browse button to select an input directory. That directory must contain the following files:

• data: A 4D series of data volumes. This will include diffusion-weighted volumes and volume(s) with no diffusion weighting.

• nodif_brain_mask: 3D binary brain mask volume derived from running bet on nodif (i.e. on a volume with no diffusion weighting).

• bvecs (with no file extension): An ASCII text file containing a list of gradient directions applied during diffusion weighted volumes.

• bvals (with no file extension): An ASCII text file containing a list of bvalues applied during each volume acquisition.

Tip: Run bedpostx_datacheck in command line to check if your input directory contains the correct files required for bedpostx.

Outputs of BEDPOSTX

bedpostx creates a new directory at the same level as the input directory called <indir>.bedpostX which contains all the files you need for probabilistic tractography. Highlights are (<i> indicates the i-th fibre. It ranges from 1 to the maximum number of fibres set in the advanced options.):

• merged_th<i>samples - 4D volume - Samples from the distribution on theta

• merged_ph<i>samples - 4D volume - Samples from the distribution on phi

  • theta and phi together represent the principal diffusion direction in spherical polar co-ordinates

• merged_f<i>samples - 4D volume - Samples from the distribution on anisotropic volume fraction (see technical report).

• mean_th<i>samples - 3D Volume - Mean of distribution on theta

• mean_ph<i>samples - 3D Volume - Mean of distribution on phi

• mean_f<i>samples - 3D Volume - Mean of distribution on f anisotropy

  • Note that in each voxel, fibres are ordered according to a decreasing mean f-value

• mean_dsamples - 3D Volume - Mean of distribution on diffusivity d

• mean_d_stdsamples - 3D Volume - Mean of distribution on diffusivity variance parameter d_std (not produced if --model=1)

• mean_S0samples - 3D Volume - Mean of distribution on T2w baseline signal intensity S0

• dyads<i> - Mean of PDD distribution in vector form. Note that this file can be loaded into FSLeyes for easy viewing of diffusion directions

• dyads<i>_dispersion - 3D Volume - Uncertainty on the estimated fibre orientation. Characterizes how wide the orientation distribution is around the respective PDD.(how is this calculated?)

• nodif_brain_mask - binary mask created from nodif_brain - copied from input directory

Advanced Options

You may change some options before running bedpostx, depending on the questions you want to ask or the quality of your diffusion data. The default values of these parameters are the ones used in the corresponding paper (Behrens et al, NeuroImage 2007).

• Fibres: Number of fibres modelled per voxel.
• Weight: Multiplicative factor for the prior on the additional modelled fibres. A smaller factor means more weighting for the additional fibres.
• Burnin: Number of iterations before starting the sampling. These might be increased if the data are noisy, and the MCMC needs more iterations to converge.

Additionally, the following alternative models are available in the advanced options:

• Single-Shell Model : Use this option if the data contains only a single non-zero bvalue. This will revert the model to use the one described in Behrens et al, NeuroImage 2007 rather than in Jbabdi et al. MRM 2012 (i.e., the the diffusion coefficient is modelled as a single value rather than using a Gamma distribution). Command line argument: --model=1.

• Model Noise Floor : Use this option if (part of) the data is operating close to the noise floor. This will fit one extra parameter f0 that captures the noise floor. This is also described in Jbabdi et al. MRM 2012. This option will output the posterior mean of the noise floor parameter mean_f0samples. Command line argument: --f0 --ardf0.

• Rician Noise : Use this option to replace the default Gaussian noise assumption with Rician noise. This will output the posterior mean of a noise precision (i.e. inverse of variance) parameter mean_tausamples. Command line argument: --rician.

Command line utility

Usage: bedpostx <subject directory> [options]

 expects to find bvals and bvecs in subject directory
 expects to find data and nodif_brain_mask in subject directory
 expects to find grad_dev in subject directory, if -g is set

 options (old syntax)
  -n (number of fibres per voxel, default 3)
  -w (ARD weight, more weight means less secondary fibres per voxel, default 1)
  -b (burnin period, default 1000)
  -j (number of jumps, default 1250)
  -s (sample every, default 25)
  -model (Deconvolution model. 1: with sticks, 2: with sticks with a range of diffusivities (default), 3: with zeppelins)
  -g (consider gradient nonlinearities, default off)

 ALTERNATIVELY: you can pass on xfibres options directly to bedpostx
  For example:  bedpostx <subject directory> --noard --cnonlinear
  Type 'xfibres --help' for a list of available options
  Default options will be bedpostx default (see above), and not xfibres default.

 Note: Use EITHER old OR new syntax.

Registration within FDT

If tractography results are to be stored in any space other than diffusion space then registration must be run.

Registration within FDT uses FLIRT, but there is also an option for using FNIRT nonlinear registration to standard space. When using the GUI, registration can only be applied after bedpostx has been run. Typically, registration will be run between three spaces:

• Diffusion space (using the nodif_brain image in the bedpostX directory.)

• Structural space (using a structural image, e.g., the space of a high-resolution T1-weighted image of the same subject)

• Standard space (by default, the MNI152 brain stored within $FSLDIR/data/standard)

Note that the structural (T1-weighted) image must have had BET applied. The nodif_brain image should be the brain extracted version of the nodif image that is stored in the bedpostX directory. Create this image using fslroi then bet on the data if it does not already exist. (it is important that the user checks the quality of bet results on these images and adjust the settings in bet where appropriate).

Transformation matrices, and their inverses, will be derived from diffusion to structural space and from structural to standard space. Relevant matrices will be concatenated to produce transformation matrices between diffusion and standard space. The resulting matrices are stored within the xfms subdirectory of the bedpostX directory and named as follows:

• diff2str.mat - from diffusion to structural space
• str2diff.mat - from structural to diffusion space
• diff2standard.mat - from diffusion to standard space
• standard2diff.mat - from standard to diffusion space
• str2standard.mat - from structural to standard space
• standard2str.mat - from standard to structural space

By default, transformation matrices between diffusion and structural space are derived using 6 degrees of freedom, the correlation ratio cost function, and normal search; transformation matrices between structural and standard space are derived using 12 degrees of freedom, the correlation ratio cost function and normal search. These parameters may be adjusted if required using the drop down menus in the registration panel.

In the GUI, it is possible to use nonlinear FNIRT registration between structural and standard space. In this case, the user needs to provide a non-brain-extracted structural image, as this will make FNIRT more accurate.

Note that diffusion images will typically have geometric distortions due to the effect of field inhomogeneities that are not affecting the structural image. Therefore, an affine FLIRT registration to structural space may not be accurate enough. We recommend that the user corrects for these distortions either using fieldmap acquisitions or using the topup tool. This must be done in the pre-processing and prior to registration.

PROBTRACKX - probabilistic tracking with crossing fibres

For details about probabilistic tractography as implemented in FDT, see this technical report, and for details about crossing fibre modelling in FDT, see Behrens et al, NeuroImage 2007.

After bedpostx has been applied it is possible to run tractography analyses using probtrackx2. Briefly, probtrackx2 produces sample streamlines, by starting from some seed and then iterate between (1) drawing an orientation from the voxel-wise bedpostX distributions, (2) taking a step in this direction, and (3) checking for any termination criteria. These sample streamlines can then be used to build up a histogram of how many streamlines visited each voxel or the number of streamlines connecting specific brain regions. This streamline distribution can be thought of as the posterior distribution on the streamline location or the connectivity distribution.

The following explains how to run probtrackx. The GUI can be started using Fdt (fdt_gui on the mac). The command-line options (which are slightly more flexible than the GUI) can be viewed using:

probtrackx2 --help

Note that the GUI only runs the new version (probtrackx2). Users of the old version (probtrackx) can only run it from the command line.

PROBTRACKX requires the specification of a bedpostX directory. This directory must contain all the following images:

• merged_th<i>samples

• merged_ph<i>samples

• merged_f<i>samples

• nodif_brain_mask

Here we will first discuss how to use seed, termination, waypoint, and exclusion masks to select a specific subset of all possible streamlines (i.e., to guide the tractograph). Afterwards, we will present the various probtrackX outputs. Finally, we have a section on using surface masks in probtrackX.

Guiding the tractography

A whole range of tools is available to guide the tractography. The most powerful of these is to select the seed region, which defines where the streamlines originate. In addition to this we can set termination masks to stop the streamlines, and waypoint/exclusion masks to filter out those streamlines not relevant for our analysis (e.g., filter out the streamlines not part of our white matter tract of interest).

Very Important: Every target mask must be in the same space as the seed masks (or the reference image in the case of a single voxel or surface seed). Any surface must also use the same convention (i.e. FreeSurfer/Caret/FIRST/voxel - see below).

The results from probtrackX will be binned in the same space as the input masks. More information on selecting this space is provided below.

Seed specification

The seed can be specified in three different ways (selected in the GUI in the Data tab):

• Single voxel seed: Generates streamlines from a single user-specified voxel. In the GUI this voxel is defined by setting the x,y,z coordinates in voxels or millimetres. From the command line this is enabled by providing a text file with the voxel coordinates to the -x,--seed flag and setting the --simple flag.

• Single mask: This is the most common setting, where the seed region is defined by a user-provided mask. Streamlines are seeded from each voxel in the mask. This mask can either be a volumetric (i.e., NIFTI) file or a surface file. On the command line this is achieved by passing on a NIFTI or surface filename to the -x,--seed flag

• Multiple mask: Streamlines are seeded from each of the masks. This can be used to generate streamlines from both a volumetric and surface mask or to produce an ROIxROI connectivity matrix. All the masks have to be in the same space and surface masks need to use the same encoding.

If the seed is a single voxel or a surface, a reference image will also need to be provided (--seedref on the command line).

The number of individual streamlines (or samples) that are drawn in each voxel or surface vertex can be set using the -P/--nsamples flag or in the Options>Basic Options section in the GUI. By default this is set to 5000 as we are confident that convergence is reached with this number of samples. However, reducing this number will speed up processing and can be useful for preliminary or exploratory analyses.

By default each streamline will be initialised in the centre of each voxel (or at the surface vertex location). This can be changed to sample in a sphere around the voxel centre or vertex using the --sampvox flag or in the Options>Advanced Options section in the GUI

Waypoint masks

If set, only streamlines that pass through the waypoint mask will be considered valid. Invalid streamlines are not counted when producing the various probtrackX outputs. These can be set in the GUI Data tab or using the --waypoints flag.

If more than one mask is set, then by default probtrackx will perform an AND operation. I.e. tracts that do not pass through ALL these masks will be discarded from the calculation of the connectivity distribution. You can set this to an OR using the Options>Waypoint option section in the GUI (--waycond on the command line). In this case, streamlines are retained if they intersect at least one of the waypoint masks. In the GUI use the add and remove buttons to make a list of waypoint masks. On the command line provide an ASCII file with the list of waypoint masks to the --waypoints flag.

Waypoint masks can be made more strict by selecting Force waypoint crossing in selected order in the Options>Waypoint option section in the GUI (--wayorder on the command line). With this option only streamlines are accepted that cross the waypoint mask in the order provided in the GUI or text file. Any streamline that reaches a waypoint out of order is terminated and rejected. This is only used when the waypoint condition is set to AND.

Exclusion mask

If an exclusion mask is to be used then check the box in the Data tab of the GUI and use the browse button to locate the mask file (or use the --avoid flag on the command line). Pathways will be discarded if they enter the exclusion mask. For example, an exclusion mask of the midline will remove pathways that cross into the other hemisphere.

Termination mask

If a termination mask is to be used then check the box and use the browse button to locate the mask file (or use the --stop flag on the command line). Pathways will be terminated as soon as they enter the termination mask. The difference between an exclusion and a termination mask is that in the latter case, the tract is stopped at the target mask, but included in the calculation of the connectivity distribution, while in the former case, the tract is completely discarded. (Note that paths are always terminated when they reach the brain surface as defined by nodif_brain_mask)

Streamlines are terminated as soon as they reach this mask. If you want to only terminate the streamlines once they leave the provided mask use the --wtstop mask instead. This can be useful to track streamlines through some termination mask without allowing them to leave. Streamlines that are seeded within the --wtstop mask are allowed to leave once without being terminated immediately.

Termination/exclusion criteria

Streamlines are terminated, when one of the following criteria is fulfilled:

• A termination or exclusion mask is used (see above).

• The streamline leaves the nodif_brain_mask in the input bedpostX directory.

• The maximum number of steps is reached. By default, samples are terminated when they have travelled 2000 steps. Using the default step length of 0.5mm this corresponds to a distance of 1m. Adjust using Options>Advanced Options GUI section or -S,--nsteps on the command line.

• The curvature of the next step is too high compared with the previous one. We limit how sharply pathways can turn in order to exclude implausible pathways. This number is the cosine of the minimum allowable angle between two steps. By default this is set to 0.2 (corresponding to a minimum angle of approximately ±80 degrees). Adjusting this number can enable pathways with sharper angles to be detected (Options>Basic Options GUI section or -c,--cthr on the command line).

• The streamline loops back on itself -i.e paths that travel to a point where they have already been. In the GUI this is enabled by default (see Options>Basic Options). It can be activated on the command line using -l,--loopcheck.

• The fractional anisotropic volumes (stored in merged_f<i>samples) can be used to influence the tractography. If activated (Options>Advanced Options in the GUI; -f,--usef on the command line) the tracts stop if the anisotropy is lower than a random variable between 0 (low anisotropy) and 1 (high anisotropy).

When terminated, the next step is to determine whether the streamline is valid (i.e., whether it will be counted in the probtrackX output). Streamlines are invalid if one of the following is true:

• The streamline reached the exclusion mask.

• The waypoint constraints were not met (see above for the various options)

• The total streamline length is shorter than some user-defined threshold. This minimum threshold is disabled by default (i.e., all streamlines pass this test). It can be set in the Options>Advanced Options GUI section or using the --disttrhresh flag on the command line.

• None of the target regions were reached when using the --network flag to compute an ROIxROI connectivity matrix (see below).

Masks in non-diffusion space

All masks must be provided in the same space (i.e., the NIFTI files must have the same size and be registered with each other). The output files will be produced in this same space. By default, probtrackX will assume that this space is diffusion space, however the masks can also be provided in a different space by selecting Seed space is not Diffusion on the GUI Data tab or using the --xfm and/or --invxfm.

If the transformation to the diffusion space is linear (e.g., for masks from a subject-specific structural image) only the transformation matrix from seed space to diffusion space has to be provided (using the --xfm flag on the command line). For non-linear transforms (e.g., for masks in standard space) the FNIRT non-linear transform has be provided both from seed space to diffusion space (--xfm) and from diffusion space to seed space (--invxfm).

For example, for analyses in structural space:

• Linear: xfms/str2diff.mat

• Nonlinear: xfms/str2diff_warp and xfms/diff2str_warp

The actual registration step to get these affine and non-linear transforms is discussed above.

Adjust streamlines path

Here we list some options available to adjust the path streamlines take (rather than seeding, terminating, and rejecting them). These are available in the Advanced Options section in the Options tab of the GUI or with the command line flags listed:

• Use modified Euler streamlining (--modeuler): Use modified Euler integration as opposed to simple Euler for computing probabilistic streamlines. More accurate but slower.

• Step length (-S,--nsteps; default 0.5mm): This determines the length of each step. This setting may be adjusted from default e.g., depending on the voxel size being used, or if tracking is being performed on different sized brains (e.g., infants or animals).

• Subsidiary fibre volume threshold (--fibthresh; default=0.01): Threshold on the volume fraction of subsidiary fibres. Fibres with a lower volume fraction than this threshold (as given in the mean_f<i>samples file in the bedpostx directory) are discarded during tractography.

ProbtrackX outputs

The seed, termination, waypoint, and exclusion masks defined above determine which streamlines will be counted in the probtrackX output files. Here we will discuss the various output files available.

A successful run of probtrackX will produce the following files in the output directory:

• probtrackx.log - a text record of the command that was run. This file can be useful to rerun the same or a slightly altered version of probtrackX.

• fdt.log - a log of the setup of the FDT GUI when the analysis was run. To recover this GUI setup, type Fdt fdt.log. This file is only produced if run from the GUI.

• waytotal - a text file containing a single number corresponding to the total number of generated tracts that have not been rejected by inclusion/exclusion mask criteria.

All probtrackX outputs count the number of streamlines in some way. The option Use distance correction (Options>Advanced Options in GUI or --pd flag) corrects for the fact that connectivity distribution (i.e., number of streamlines) drops with distance from the seed mask. If this option is checked, the connectivity distribution is the expected length of the pathways that cross each voxel times the number of samples that cross it.

In addition to counting the number of streamlines, the average length of the streamlines up to that point can also be stored by setting the --ompl flag in addition to setting at least one of the output listed below.

Streamline density map

This is the default output when running from the GUI (activate using the --opd flag in the probtrackX output). The result is a 3D image file (called fdt_paths) containing the number of streamlines reaching each voxel. This can be thought of as quantifying the connectivity from the seed region.

By carefully selecting the seed region, termination mask, and inclusion/exclusions masks, this can be used to reconstruct specific white matter tracts. XTRACT leverages probtrackX and a set of carefully curated masks to reconstruct many of the major white matter tracts.

ROI by ROI connectivity matrix

The output here quantifies the number of streamlines seeded in one region of interest (ROI) that reached another ROI. For N ROIs, the result will be an NxN matrix (called fdt_network_matrix), where each row quantifies how many streamlines seeded in a certain ROI reached the other ROIs.

In the GUI this output can be produced by providing the N ROIs as multiple seed masks (see above for instructions). An ASCII text file listing the ROIs can also be provided on the command line. On the command line the --network flag should also be added to produce the NxN matrix. Note that adding the --network flag will cause any streamlines that are seeded in one ROI and do not reach any other ROIs to be marked as invalid. This will alter any other probtrackX output produced.

voxel by ROI connectivity

To estimate the voxelxROI connectivity probtrackx quantifies connectivity values between each voxel in a seed mask and any number of user-specified target masks. In the example on the right, seed voxels in the thalamus are classified according to the probability of connection to different cortical target masks.

In the GUI this output is created by selecting Classification Targets on the Data tab (note that this option is only available if the seed space is set to Single Mask). Use the add button to locate each target mask. Targets must be in the same space as the seed mask. When all targets are loaded you can press the save list button to save the list of targets as a text file. If you already have a text file list of required targets (including their path) then you can load it with the load list button or set the --targetmasks and --os2t flags on the command line.

The output directory will contain a single volume for each target mask, named seeds_to_{target} where {target} is replaced by the file name of the relevant target mask. In these output images, the value of each voxel within the seed mask is the number of samples seeded from that voxel reaching the relevant target mask. The value of all voxels outside the seed mask will be zero. Note that the files seeds_to_{target} will have the same format (either volumes or surfaces) as the seed mask.

There are command line utilities that can be run on these outputs:

• proj_thresh - for thresholding some outputs of probtrackx
• find_the_bigggest - for performing hard segmentation on the outputs of probtrackx, see example on the right

Voxel by voxel connectivity matrix

Voxel by voxel connectivity matrices are available in the GUI through the Matrix Options in the Options tab or using --matrix1/2/3 and associated flags in the command line (see probtrackx2 --help). These can be used to store a connectivity matrix between (1) all seed points and all other seed points, or (2) all seed points and all points in a target mask, or (3) all pairs of points in a target mask (or a pair of target masks). All masks used here MUST BE IN SEED SPACE. These options can be used separately or in conjunction.

Matrix1 and Matrix2 store the number of samples (potentially modulated by distance when distance correction is used) between the rows (seed points) and the columns of the matrix.

Matrix3 stores the number of sample streamlines (potentially modulated by distance when distance correction is used) between each pair of points in a pair of target masks. In this case, when a streamline is sent from each seed point in both directions, if the streamline hits the two target masks at two locations along either sides of the streamline, the corresponding row and column of matrix3 is filled.

Note on the ouput format: matrix1,2,3 are stored as 3 column ASCII coding of sparse matrices. These files can be loaded into matlab using e.g.:

    x=load('fdt_matrix1.dot');
    M=spconvert(x);

Typical uses of Matrix options:

• Matrix1: The seed mask(s) can represent grey matter, so this would be all GM to GM connectivity.
• Matrix2: The seed can be a grey matter region, and mask2 the rest of the brain (can be low res). The results can then be used for blind classification of the seed mask.
• Matrix3: The 'target' masks can be the whole of grey matter, and the seed mask white matter. This option allows more sensitive reconstruction of grey-to-grey connectivity as pathways are seeded from all their constituting locations, rather than just their end-points as in Matrix1.

Example use of Matrix2 for clustering:

A typical use of Matrix2 is blind (i.e. hypothesis-free) classification. Say you have run probtrackx with a seed roi (e.g. thalamus) and the --omatrix2 option by setting the Matrix2 target mask to a whole brain mask (typically this mask would be lower resolution than the seed mask). The resulting Matrix2 can be used in Matlab to perform 'kmeans' classification. Below is an example Matlab code to do this (this example is only valid if the seed mask is a NIFTI volume, not a surface):

   % Load Matrix2
   x=load('fdt_matrix2.dot');
   M=full(spconvert(x));
   % Calculate cross-correlation
   CC  = 1+corrcoef(M');
   % Do kmeans with k clusters
   idx = kmeans(CC,k);   % k is the number of clusters
   % Load coordinate information to save results
   addpath([getenv('FSLDIR') '/etc/matlab']);
   [mask,~,scales] = read_avw('fdt_paths');
   mask = 0*mask;
   coord = load('coords_for_fdt_matrix2')+1;
   ind   = sub2ind(size(mask),coord(:,1),coord(:,2),coord(:,3));
   [~,~,j] = unique(idx);
   mask(ind) = j;
   save_avw(mask,'clusters','i',scales);
   !fslcpgeom fdt_paths clusters

Using surfaces

It is possible to use surface files in probtrackx. The preferred format is GIFTI.

Surface files typically describe vertex coordinates in mm. Annoyingly, different softwares use different conventions to transform mm to voxel coordinates. Probtrackx can deal with the following conventions: freesurfer, caret, first and voxel. The latter convention simply means that mm and voxel coordinates are the same.

Switching conventions

When using surfaces in probtrackx, all surfaces MUST use the same convention. However, we provide a command-line tool (surf2surf) to transform between different conventions. This tool can also be used to convert between different surface file formats.

surf2surf -i <inputSurface> -o <outputSurface> [options]

Compulsory arguments (You MUST set one or more of):
        -i,--surfin     input surface
        -o,--surfout    output surface

Optional arguments (You may optionally specify one or more of):
        --convin        input convention [default=caret] - only used if output convention is different
        --convout       output convention [default=same as input]
        --volin         input ref volume - Must set this if changing conventions
        --volout        output ref volume [default=same as input]
        --xfm           in-to-out ascii matrix or out-to-in warpfield [default=identity]
        --outputtype    output type: ASCII, VTK, GIFTI_ASCII, GIFTI_BIN, GIFTI_BIN_GZ (default)
        --values        set output scalar values (e.g. --values=mysurface.func.gii or --values=1)

Projecting data onto the surface

It is often useful to project 3D data onto the cortical surface. We provide the following command-line tool (surf_proj) to do exactly that:

Usage: surf_proj [options]

Compulsory arguments (You MUST set one or more of):
        --data  data to project onto surface
        --surf  surface file
        --out   output file

Optional arguments (You may optionally specify one or more of):
        --meshref       surface volume ref (default=same as data)
        --xfm   data2surf transform (default=Identity)
        --meshspace     meshspace (default='caret')
        --step  average over step (mm - default=1)
        --direction     if>0 goes towards brain (default=0 ie both directions)
        --operation     what to do with values: 'mean' (default), 'max', 'median', 'last'
        --surfout       output surface file, not ascii matrix (valid only for scalars)

Using FreeSurfer surfaces

In order to use FreeSurfer-generated surfaces as masks in probtrackx, you need to add the following command line options (which can also be set in the GUI):

   --meshspace=freesurfer --seedref=orig.nii.gz

Where the orig.nii.gz file can be generated by converting the orig.mgz file (found in the FreeSurfer output directories) using FreeSurfer's mri_convert.

Additionally, using FreeSurfer with probtrackx requires a few specific steps that we describe below.

FreeSurfer Registration

FreeSurfer operates in a conformed space, which is different from the original structural image space that it has received as an input. When tracking from FreeSurfer surfaces, it is necessary to provide a transformation between conformed space and diffusion space. Below are a few steps that show how to achieve this.

We assume that you have ran dtifit on your diffusion data with an FA map called dti_FA.nii.gz (we recommend using an FA map to register to T1 structural images), and also that you have a file called struct.nii.gz that you have used as an input to FreeSurfer recon_all program.

We will carry on a few steps that aim at calculating the following transformations: fa<->struct<->freesurfer. Then we will concatenate these transformations to get fa<->freesurfer.

Let us start with struct<->freesufer (assuming john is the subject's name):

tkregister2 --mov $SUBJECTS_DIR/john/mri/orig.mgz --targ $SUBJECTS_DIR/john/mri/rawavg.mgz --regheader --reg junk --fslregout freesurfer2struct.mat --noedit
convert_xfm -omat struct2freesurfer.mat -inverse freesurfer2struct.mat

Now transforming FA to struct:

flirt -in dti_FA -ref struct_brain -omat fa2struct.mat
convert_xfm -omat struct2fa.mat -inverse fa2struct.mat

The final stage is to concatenate these transformations:

convert_xfm -omat fa2freesurfer.mat -concat struct2freesurfer.mat fa2struct.mat
convert_xfm -omat freesurfer2fa.mat -inverse fa2freesurfer.mat

Label files

Label files from FreeSurfer are useful as cortical ROIs for tractography. In order to use a label file (or collection of labels) in probtrackx, you must first transform it (them) into a surface file, using the label2surf command.

For instance, say you want to create a surface ROI for Brodmann areas 44 and 45. You must first decide which surface to use (pial, white) and which hemisphere (lh or rh). Say you want BA44/45 on the lh.white surface, then you must first transform the relevant surface file into an FSL-compatible format (ASCII or VTK or GIFTI), then use label2surf to transform the label files into a surface file. Below, the first command is a FreeSurfer command:

mris_convert lh.white lh.white.gii
echo lh.BA44.label lh.BA45.label > listOfAreas.txt
label2surf -s lh.white.gii -o lh.BA44.gii -l listOfAreas.txt

The full set of arguments for label2surf is shown below:

label2surf
         Transforms a group of labels into a surface

Usage:
label2surf -s <surface> -o <outputsurface> -l <labels>

Compulsory arguments (You MUST set one or more of):
        -s,--surf       input surface
        -o,--out        output surface
        -l,--labels     ascii list of label files

Optional arguments (You may optionally specify one or more of):
        -v,--verbose    switch on diagnostic messages
        -h,--help       display this message

Utilities

proj_thresh

proj_thresh is a command line utility that provides an alternative way of expressing connection probability in connectivity-based segmentation. It is run on the output of probtrackx when classification targets are used.

The output of Connectivity-based seed classification is a single volume for each target mask, named seeds_to_{target} where {target} is replaced by the file name of the relevant target mask. In these output images, the value of each voxel within the seed mask is the number of samples seeded from that voxel reaching the target mask. proj_thresh is run as follows:

proj_thresh <list of volumes/surfaces> threshold

Where the list of volumes is the outputs of Connectivity-based seed classification (i.e., files named seeds_to_target1 etc etc) and threshold is expressed as a number of samples For each voxel in the seeds mask that has a value above threshold for at least one target mask, proj_thresh calculates the number of samples reaching each target mask as a proportion of the total number of samples reaching any target mask. The output of proj_thresh is a single volume for each target mask.

find_the_biggest

find_the_biggest is a command line utility that performs hard segmentation of a seed region on the basis of outputs of probtrackx when classification targets are being used.

The output of Connectivity-based seed classification is a single volume for each target mask, named seeds_to_{target} where {target} is replaced by the file name of the relevant target mask. In these output images, the value of each voxel within the seed mask is the number of samples seeded from that voxel reaching the target mask. find_the_biggest classifies seed voxels according to the target mask with which they show the highest probability of connection. It is run as follows:

find_the_biggest <list of volumes/surfaces> <output>

Where the list of volumes is the outputs of Connectivity-based seed classification (i.e., files named seeds_to_target1 etc etc).

The example below uses probtrackx and find_the_biggest to perform hard segmentation of the thalamus on the basis of its connections to cortex.

vecreg - Registration of vector images

After running dtifit or bedpostx, it is often useful to register vector data to another space. For example, one might want to represent V1 for different subjects in standard space. vecreg is a command line tool that allows to perform such registration.

Vector images cannot be registered by simply applying a transformation (as calculated by, say, FLIRT) to every voxel's coordinates. The corresponding vectors have to be reoriented accordingly (see D. Alexander 2001, IEEE-TMI 20:1131-39). vecreg performs this operation for you. The image on the right shows the effect of applying vecreg (right) to the V1 image on the left, compared to simply applying voxelwise transformation (e.g. using applyxfm4D) to the vectors (centre).

Important: vecreg does not calculate a transformation, but simply applies a given transformation to the input vector field. vecreg can apply a linear transformation calculated with FLIRT, or a non-linear transformation calculated by FNIRT.

types of input that may be used for vecreg from DTIFIT: V1,V2,V3,tensor from BEDPOSTX: dyads1, dyads2, etc.

Command line options

vecreg -i <input4D> -o <output4D> -r <refvol> [-t <transform>]


Compulsory arguments (You MUST set one or more of):
        -i,--input      filename for input vector or tensor field
        -o,--output     filename for output registered vector or tensor field
        -r,--ref        filename for reference (target) volume

Optional arguments (You may optionally specify one or more of):
        -v,--verbose    switch on diagnostic messages
        -h,--help       display this message
        -t,--affine     filename for affine transformation matrix
        -w,--warpfield  filename for 4D warp field for nonlinear registration
        --rotmat        filename for secondary affine matrix
                        if set, this will be used for the rotation of the vector/tensor field
        --rotwarp       filename for secondary warp field
                        if set, this will be used for the rotation of the vector/tensor field
        --interp        interpolation method : nearestneighbour, trilinear (default), sinc or spline
        -m,--mask       brain mask in input space
        --refmask       brain mask in output space (useful for speed up of nonlinear reg)

qboot - Estimation of fibre orientations using q-ball ODFs and residual bootstrap

qboot is a command line tool that allows estimation of diffusion ODFs and fibre orientations from them. Its output can be used as an input for probtrackX in order to perform probabilistic tractography.

ODF estimation is performed using a real spherical harmonics basis. Fibre orientations are estimated as the local maxima of the ODFs. Both deterministic and probabilistic estimation can be performed. For the latter, residual bootstrap is performed to infer on the ODF shape and obtain a distribution of fibre orientations. For more details on the implementation see Sotiropoulos2011 (S.N. Sotiropoulos, I. Aganj, S. Jbabdi, G. Sapiro, C. Lenglet and T.E. Behrens, "Inference on Constant Solid Angle Orientation Distribution Functions from Diffusion-Weighted MRI", p.609, Quebec, Canada, OHBM, 2011).

qboot allows reconstruction of q-ball ODFs (Tuch DS, MRM 2004), CSA ODFs (Aganj I et al, MRM, 2010) and variants of them, obtained via Laplacian sharpening and Laplace-Beltrami regularization (Descoteaux et al, MRM, 2007). Both spherical harmonic coefficients of the reconstructed ODFs and fibre orientation estimates may be returned as output. A real spherical harmonic basis is employed (Aganj I et al, MRM, 2010).

Input files for qboot : Similar to dtifit and bedpostx, qboot needs a 4D data file, a binary mask_file, a bvecs and a bvals file.

Command-line options

qboot -k data_file -m nodif_brain_mask -r bvecs -b bvals

Compulsory arguments (You MUST set one or more of):
        -k,--data       Data file
        -m,--mask       Mask file
        -r,--bvecs      b vectors file
        -b,--bvals      b values file

Optional arguments (You may optionally specify one or more of):
        --ld,--logdir   Output directory (default is logdir)
        --forcedir      Use the actual directory name given - i.e. don't add + to make a new directory
        --q             File provided with multi-shell data. Indicates the number of directions for each shell
        --model         Which model to use. 1=Tuch's ODFs, 2=CSA ODFs (default), 3=multi-shell CSA ODFs
        --lmax          Maximum spherical harmonic order employed (must be even, default=4)
        --npeaks        Maximum number of ODF peaks to be detected (default 2)
        --thr           Minimum threshold for a local maxima to be considered an ODF peak.
                        Expressed as a fraction of the maximum ODF value (default 0.4)
        --ns,--nsamples Number of bootstrap samples (default is 50)
        --lambda        Laplace-Beltrami regularization parameter (default is 0)
        --delta         Signal attenuation regularization parameter for model=2 (default is 0.01)
        --alpha         Laplacian sharpening parameter for model=1 (default is 0, should be smaller than 1)
        --seed          Seed for pseudo-random number generator
        --gfa           Compute a generalised FA, using the mean ODF in each voxel
        --savecoeff     Save the ODF coefficients instead of the peaks.
        --savemeancoeff Save the mean ODF coefficients across all samples
        -V,--verbose    Switch on diagnostic messages

Possible Outputs of qboot

• qboot will by default return distributions of ODF peaks. Peaks are detected as local ODF maxima using a finite-difference scheme. Setting --ns=1 , will perform deterministic estimation of the ODFs and their peaks.

• If --savecoeff is set, then only samples of the ODF coefficients will be saved, without any peaks.

• If --savecoeff --savemeancoeff are set, then only the mean ODF coefficients (across the samples) will be saved for each voxel and again no peaks.

• If --savemeancoeff is set, then the mean ODF coefficients (across the samples) will be saved for each voxel, along with the samples of the ODF peaks.

• If none of the above is set, samples of the ODF peaks are saved.

Running qboot on a cluster

Similar to bedpostX qboot can be parallelised if run on an SGE-capable system. The qboot_parallel script can be employed for this purpose.

Multi-shell data assumptions

The current implementation of qboot can estimate multi-shell ODFs, assuming the following for the data:

• Three non-zero shell data are assumed. The bvalues should form an arithmetic progression, e.g. 0, 1000, 2000, 3000 or 0, 2000, 4000, 6000.

It is assumed that all data from each shell are grouped together and shells are one after the other in data, bvecs and bvals. So, if for example 3 directions are available for 3 shells these should appear as: dir1_shell1, dir2_shell1, dir3_shell1, dir1_shell2,dir2_shell2, dir3_shell2, dir1_shell3,dir2_shell3, dir3_shell3.

• There is no assumption on the number of directions in each shell. Each shell can have its own directions and number of directions. The minimum number of directions in a shell is dictated by the number of coefficients estimated, i.e. if lmax=4 => Num_of_coeff=15 => All shells should have at least 15 directions. This is because interpolation of the data occurs, when: i) different number of directions exist in different shells or ii) same number of directions, but corresponding directions between shells (e.g. dir1_shell1 vs dir1_shell2) differ more than 1 degree. In this case each shell is fitted with SH of the same order, which is by default 10, unless the number of directions in the shell with the less points is not enough (then order is decreased).

• By default, same number of datapoints (including b=0's) are assumed for each shell acquisition. If this is not the case, a --q=qshells.txt file should be provided. The information in the qshells file is of the form N1 N2 N3 , each number indicating the number of datapoints (including b=0's) that correspond to each shell acquisition.

FSLeyes

Displaying DWI images in FSLeyes

Outputs of bedpostx or dtifit can be conveniently displayed in FSLeyes. If you open an image of diffusion vectors (e.g., dtifit_V1 output of dtifit or dyads<i> output of bedpostx) then it is possible to display these vectors using RGB coding (where the colours red,green and blue represent diffusion in the x,y,z axes respectively) or using lines where a line at each voxel represents the principle diffusion direction at that voxel. To do this, set the Overlay type to 3-direction vector image (RGB) or 3-direction vector image (Line).

It is also possible to display multiple lines per voxel - simply load the vector images you wish to view, and set them all to be displayed as line vectors. You can change the colours used to display each image via the overlay display panel (the gear button on the top toolbar). In the example above, the red lines represent the main fibre orientation, and the blue lines represent the secondary fibre orientations (thresholded at an f-value of 0.05), as calculated by bedpostx.

If you only want to visualise fibres within voxels where multiple fibres are supported (e.g. above a certain threshold for the corresponding mean_f<i>samples), you can "clip" the vector image so that only voxels above a certain threshold are shown. Load both the dyads<i> and mean_f<i>samples images, then select the dyads image and open the overlay display panel (the gear button). Set the Clip by option to the mean_f<i>samples image, and then adjust the Clipping range to the threshold you want.

Another option is to create a vector file where the voxels below a certain f-threshold are zeroed. To do this you can use maskdyads:

maskdyads dyads<i> mean_f<i>samples 0.05

FSLeyes can also display the tensor model fit from output of dtifit - choose File -> Add from directory, and select the dtifit output directory. Or, if you have used the --save_tensor option to dtifit, you can load the dtifit_tensor file and set the Overlay type to Diffusion tensor.

For more details see the FSLeyes documentation.

CategoryFDT