Fibre density and cross-section - Multi-tissue CSD


This tutorial explains how to perform fixel-based analysis of fibre density and cross-section with fibre orientation distributions (FODs) computing using multi-tissue CSD using single-shell data or multi-shell data. We note that high b-value (>2000s/mm2) data is recommended to aid the interpretation of AFD being related to the intra-axonal space. See the original paper for more details.

All steps in this tutorial have written as if the commands are being run on a cohort of images, and make extensive use of the foreach script to simplify batch processing. This tutorial also assumes that the imaging dataset is organised with one directory identifying the subject, and all files within identifying the image type. For example:



All commands in this tutorial are run from the subjects path up until step 18, where we change directory to the template path

For all MRtrix scripts and commands, additional information on the command usage and available command-line options can be found by invoking the command with the -help option. Please post any questions or issues on the MRtrix community forum.

Pre-processsing steps

1. DWI denoising

The effective SNR of diffusion data can be improved considerably by exploiting the redundancy in the data to reduce the effects of thermal noise. This functionality is provided in the command dwidenoise:

foreach * : dwidenoise IN/dwi.mif IN/dwi_denoised.mif

Note that this denoising step must be performed prior to any other image pre-processing: any form of image interpolation (e.g. re-gridding images following motion correction) will invalidate the statistical properties of the image data that are exploited by dwidenoise, and make the denoising process prone to errors. Therefore this process is applied as the very first step.

2. DWI general pre-processing

The dwipreproc script is provided for performing general pre-processing of diffusion image data - this includes eddy current-induced distortion correction, motion correction, and (possibly) susceptibility-induced distortion correction. Commands for performing this pre-processing are not yet implemented in MRtrix3; the dwipreproc script in its current form is in fact an interface to the relevant commands that are provided as part of the FSL package. Installation of FSL (including eddy) is therefore required to use this script, and citation of the relevant articles is also required (see the dwipreproc help page).

Usage of this script varies depending on the specific nature of the DWI acquisition with respect to EPI phase encoding - full details are available within the DWI distortion correction using ``dwipreproc` ` page, and the dwipreproc help file.

Here, only a simple example is provided, where a single DWI series is acquired where all volumes have an anterior-posterior (A>>P) phase encoding direction:

foreach * : dwipreproc IN/dwi_denoised.mif IN/dwi_denoised_preproc.mif -rpe_none -pe_dir AP

Fixel-based analysis steps

3. Computing group average tissue response functions

As described here, using the same response function when estimating FOD images for all subjects enables differences in the intra-axonal volume (and therefore DW signal) across subjects to be detected as differences in the FOD amplitude (the AFD). To ensure the response function is representative of your study population, a group average response function can be computed by first estimating a response function per subject, then averaging with the script:

foreach * : dwi2response dhollander IN/dwi_denoised_preproc.mif IN/response_wm.txt IN/response_gm.txt IN/response_csf.txt
average_response */response_wm.txt ../group_average_response_wm.txt
average_response */response_gm.txt ../group_average_response_gm.txt
average_response */response_csf.txt ../group_average_response_csf.txt

4. Upsampling DW images

Upsampling DWI data before computing FODs can increase anatomical contrast and improve downstream spatial normalisation and statistics. We recommend upsampling to a voxel size of 1.25mm (for human brains). If you have data that has smaller voxels than 1.25mm, then we recommend you can skip this step:

foreach * : mrresize IN/dwi_denoised_preproc.mif -vox 1.25 IN/dwi_denoised_preproc_upsampled.mif

5. Compute upsampled brain mask images

Compute a whole brain mask from the upsampled DW images:

foreach * : dwi2mask IN/dwi_denoised_preproc_upsampled.mif IN/dwi_mask_upsampled.mif

Depending on your data, you may find that computing masks on native resolution DWIs gives superiour masks (with less holes), then upsampling. This can be performed using:

foreach * : dwi2mask IN/dwi_denoised_preproc.mif - \| mrresize - -scale 2.0 -interp nearest IN/dwi_mask_upsampled.mif

6. Fibre Orientation Distribution estimation

When performing analysis of AFD, Constrained Spherical Deconvolution (CSD) should be performed using the group average response functions computed at step 3:

foreach * : dwi2fod msmt_csd IN/dwi_denoised_preproc_upsampled.mif ../group_average_response_wm.txt IN/fod.mif ../group_average_response_gm.txt IN/gm.mif  ../group_average_response_csf.txt IN/csf.mif -mask IN/dwi_mask_upsampled.mif

7. Perform simultaneous bias field correction and intensity normalisation

This step performs global intensity normalisation by scaling all tissue types based on a single scale factor. A single multiplicative bias field is also estimated and applied to correct the output:

foreach * : mtbin IN/fod.mif IN/fod_bias_norm.mif IN/gm.mif IN/gm_bias_norm.mif IN/csf.mif IN/csf_bias_norm.mif


We strongly recommend you that you check the scale factors applied during intensity normalisation are not influenced by the variable of interest in your study. For example if one group contains global changes in white matter T2 then this may directly influence the intensity normalisation and therefore bias downstream AFD analysis. To check this we recommend you perform an equivalence test to ensure mean scale factors are the same between groups. To output the scale factor applied for all subjects use foreach * : mrinfo IN/fod_bias_norm.mif -property normalisation_scale_factor.

8. Generate a study-specific unbiased FOD template

Population template creation is the most time consuming step in a fixel-based analysis. If you have a large number of subjects in your study, we recommend building the template from a subset of 20-40 individuals. Subjects should be chosen to ensure the generated template is representative of your population (i.e. equal number of patients and controls). To build a template, place all FOD images in a single folder. We also recommend placing a set of corresponding mask images (with the same prefix as the FOD images) in another folder. Using masks can speed up registration significantly:

mkdir -p ../template/fod_input
mkdir ../template/mask_input

Symbolic link all FOD images (and masks) into a single input folder. If you have fewer than 40 subjects in your study, you can use the entire population to build the template:

foreach * : ln -sr IN/fod_bias_norm.mif ../template/fod_input/PRE.mif
foreach * : ln -sr IN/dwi_mask_upsampled.mif ../template/mask_input/PRE.mif

Alternatively, if you have more than 40 subjects you can randomly select a subset of the individuals. If your study has multiple groups, then ideally you want to select the same number of subjects from each group to ensure the template is un-biased. Assuming the subject directory labels can be used to identify members of each group, you could use:

foreach `ls -d *patient | sort -R | tail -20` : ln -sr IN/fod_bias_norm.mif ../template/fod_input/PRE.mif ";" ln -sr IN/dwi_mask_upsampled.mif ../template/mask_input/PRE.mif
foreach `ls -d *control | sort -R | tail -20` : ln -sr IN/fod_bias_norm.mif ../template/fod_input/PRE.mif ";" ln -sr IN/dwi_mask_upsampled.mif ../template/mask_input/PRE.mif

Run the template building script as follows:

population_template ../template/fod_input -mask_dir ../template/mask_input ../template/fod_template.mif

If you are building a template from your entire study population, run the population_template script use the -warp_dir warps option to output a directory containing all subject warps to the template. Saving the warps here will enable you to skip the next step. Note that the warps used (and therefore output) from the population_template script are 5D images containing both forward and reverse warps (see :ref:`mrregister`for more info). After population template creation is complete, to convert this warp format to a more conventional 4D deformation field format ready for the subsequent steps, run:

foreach ../template/warps/* : warpconvert -type warpfull2deformation -template ../template/fod_template.mif IN PRE/subject2template_warp.mif

9. Register all subject FOD images to the FOD template

Register the FOD image from all subjects to the FOD template image. Note you can skip this step if you built your template from your entire population and saved the warps (see previous step):

foreach * : mrregister IN/fod_wm_bias_norm.mif -mask1 IN/dwi_mask_upsampled.mif ../template/fod_template.mif -nl_warp IN/subject2template_warp.mif IN/template2subject_warp.mif

10. Compute the intersection of all subject masks in template space

Different subjects will have subtly different brain coverage. To ensure subsequent analysis is performed in voxels that contain data from all subjects, we warp all subject masks into template space and compute the mask intersection. For each subject:

foreach * : mrtransform IN/dwi_mask_upsampled.mif -warp IN/subject2template_warp.mif -interp nearest IN/dwi_mask_in_template_space.mif

Compute the intersection of all warped masks:

mrmath */dwi_mask_in_template_space.mif min ../template/mask_intersection.mif

11. Compute a white matter analysis voxel & fixel mask

Here we first identify all voxels having some white matter by thresholding the DC term (first SH coefficient) of the multi-tissue FOD image:

mrconvert ../template/fod_template.mif -coord 3 0 - | mrthreshold - ../template/voxel_mask.mif

Next we segment all fixels from each FOD in the template image (see here for more information about a analysis fixel mask). Note that the fixel image output from this step is stored using the Fixel image (directory) format, which exploits the filesystem to store all fixel data in a directory:

fod2fixel -mask ../template/voxel_mask.mif -fmls_peak_value 0.2 ../template/fod_template.mif ../template/fixel_mask

You can visualise the output fixels using the fixel plot tool from mrview, and opening either the index.mif or directions.mif found in ../template/fixel_mask. The automatic thresholding step used above should give you a mask that nicely covers all of white matter, however if not you can always try manually adjusting the threshold with the mrthreshold -abs option.


We recommend having no more than 500,000 fixels in the analysis_fixel_mask (you can check this by mrinfo -size ../template/fixel_mask/directions.mif, and looking at the size of the image along the 1st dimension), otherwise downstream statistical analysis (using fixelcfestats) will run out of RAM). A mask with 500,000 fixels will require a PC with 128GB of RAM for the statistical analysis step. To reduce the number of fixels, try either reducing the number of voxels in the voxel mask by applying a manual threshold using -abs, increasing the -fmls_peak_value, or reducing the extent of upsampling in step 4.

12. Warp FOD images to template space

Note that here we warp FOD images into template space without FOD reorientation. Reorientation will be performed in a separate subsequent step:

foreach * : mrtransform IN/fod_bias_norm.mif -warp IN/subject2template_warp.mif -noreorientation IN/fod_in_template_space.mif

13. Segment FOD images to estimate fixels and their apparent fibre density (FD)

Here we segment each FOD lobe to identify the number and orientation of fixels in each voxel. The output also contains the apparent fibre density (AFD) value per fixel estimated as the FOD lobe integral (see here for details on FOD segmentation). Note that in the following steps we will use a more generic shortened acronym - Fibre Density (FD) instead of AFD, since the following steps can also apply for other measures of fibre density (see the note below). The terminology is also consistent with our recent work:

foreach * : fod2fixel IN/fod_in_template_space.mif -mask ../template/voxel_mask.mif IN/fixel_in_template_space -afd fd.mif

14. Reorient fixel orientations

Here we reorient the direction of all fixels based on the Jacobian matrix (local affine transformation) at each voxel in the warp. Note that in-place fixel reorientation can be performed by specifing the output fixel folder to be the same as the input, and using the -force option:

foreach * : fixelreorient IN/fixel_in_template_space IN/subject2template_warp.mif IN/fixel_in_template_space --force

15. Assign subject fixels to template fixels

In step 8 & 9 we obtained spatial correspondence between subject and template. In step 14 we corrected the fixel orientations to ensure angular correspondence of the segmented peaks of subject and template. Here, for each fixel in the template fixel analysis mask, we identify the corresponding fixel in each voxel of the subject image and assign the FD value of the subject fixel to the corresponding fixel in template space. If no fixel exists in the subject that corresponds to the template fixel then it is assigned a value of zero. See this paper for more information. In the command below, you will note that the output fixel directory is the same for all subjects. This directory now stores data for all subjects at corresponding fixels, ready for input to fixelcfestats in step 20 below:

foreach * : fixelcorrespondence IN/fixel_in_template_space/fd.mif ../template/fixel_mask ../template/fd PRE.mif

16. Compute fibre cross-section (FC) metric

Apparent fibre density, and other related measures that are influenced by the quantity of restricted water, only permit the investigation of group differences in the number of axons that manifest as a change to within-voxel density. However, depending on the disease type and stage, changes to the number of axons may also manifest as macroscopic differences in brain morphology. This step computes a fixel-based metric related to morphological differences in fibre cross-section, where information is derived entirely from the warps generated during registration (see this paper for more information):

foreach * : warp2metric IN/subject2template_warp.mif -fc ../template/fixel_mask ../template/fc IN.mif

Note that the FC files will be used in the next step. However, for group statistical analysis of FC we recommend taking the log(FC) to ensure data are centred about zero and normally distributed. We could place all the log(FC) fixel data files in the same fixel directory as the FC files (as long as they are named differently. However to keep things tidy, create a separate fixel directory to store the log(FC) data and copy the fixel index and directions file across:

mkdir ../template/log_fc
cp ../template/fc/index.mif ../template/fc/directions.mif ../template/log_fc
foreach * : mrcalc ../template/fc/IN.mif -log ../template/log_fc/IN.mif

17. Compute a combined measure of fibre density and cross-section (FDC)

To account for changes to both within-voxel fibre density and macroscopic atrophy, fibre density and fibre cross-section must be combined (a measure we call fibre density & cross-section, FDC). This enables a more complete picture of group differences in white matter. Note that as discussed in this paper, group differences in FD or FC alone must be interpreted with care in crossing-fibre regions. However group differences in FDC are more directly interpretable. To generate the combined measure we ‘modulate’ the FD by FC:

mkdir ../template/fdc
cp ../template/fc/index.mif cp ../template/fc/directions.mif ../template/fdc
foreach * : mrcalc ../template/fd/IN.mif ../template/fc/IN.mif -mult ../template/fdc/IN.mif

18. Perform whole-brain fibre tractography on the FOD template

Statistical analysis using connectivity-based fixel enhancement exploits connectivity information derived from probabilistic fibre tractography. To generate a whole-brain tractogram from the FOD template. Note the remaining steps from here on are executed from the template directory:

cd ../template
tckgen -angle 22.5 -maxlen 250 -minlen 10 -power 1.0 fod_template.mif -seed_image voxel_mask.mif -mask voxel_mask.mif -select 20000000 tracks_20_million.tck

19. Reduce biases in tractogram densities

Perform SIFT to reduce tractography biases in the whole-brain tractogram:

tcksift tracks_20_million.tck fod_template.mif tracks_2_million_sift.tck -term_number 2000000

20. Perform statistical analysis of FD, FC, and FDC

Statistical analysis is performed using connectivity-based fixel enhancement, with a separate analysis for FD, FC, and FDC as follows:

fixelcfestats fd files.txt design_matrix.txt contrast_matrix.txt input_tracks_2_million_sift.tck stats_fd
fixelcfestats log_fc files.txt design_matrix.txt contrast_matrix.txt input_tracks_2_million_sift.tck stats_log_fc
fixelcfestats fdc files.txt design_matrix.txt contrast_matrix.txt input_tracks_2_million_sift.tck stats_fdc

Where the input files.txt is a text file containing the filename of each file (i.e. not the full path) to be analysed inside the input fixel directory, each filename on a separate line. The line ordering should correspond to the lines in the design_matrix.txt. Note that for correlation analysis, a column of 1’s will not be automatically included (as per FSL randomise). Note that fixelcfestats currently only accepts a single contrast. However if the opposite (negative) contrast is also required (i.e. a two-tailed test), then use the -neg option. Several output files will generated all starting with the supplied prefix.


We recommend having no more than 500,000 fixels in the analysis_fixel_mask (you can check this by mrinfo -size ../template/fixel_template/mask.mif, and looking at the size of the image along the 1st dimension), otherwise fixelcfestats will run out of RAM. A mask with 500,000 fixels will require a PC with 128GB of RAM for the statistical analysis step. To reduce RAM requirements, you could reduce the number of fixels by reducing the extent of upsampling, or apply a higher threshold when computing the fixel analysis mask (at the risk of removing WM regions from your analysis).

21. Visualise the results

To view the results load the population FOD template image in mrview, and overlay the fixel images using the vector plot tool. Note that p-value images are saved as 1-p-value. Therefore to visualise all p-values < 0.05, threshold the fixels using the fixel plot tool at 0.95.