Multi-fibre/HARDI Reconstruction

1. Overview

This tutorial deals with reconstruction of diffusion MRI data using multiple-fibre algorithms and voxel classification. We start by looking at the multi-tensor model of diffusion and the use of voxel classification as a way of improving the procedure. We then go on to look at multiple fibre reconstructions (specifically QBall and PASMRI) and the Camino functions that can be used to extract useful information/statistics from them.

Multiple fibre reconstruction can be quite time consuming. Therefore, prior to attempting a multiple-fibre reconstruction, we highly recommend checking the data for artefacts and other problems and then fitting the diffusion tensor (see diffusion tensor imaging tutorial) to make sure that there are no issues with scheme file. You should only proceed when you are completely satisfied that the data and schemefile are in order.

2. Voxel Classification

The command voxelclassify classifies each voxel as isotropic, anisotropic Gaussian or non-Gaussian using the algorithm in [Alexander, Barker and Arridge, MRM, Vol 48, pp 331-340, 2002]. In isotropic and anisotropic Gaussian voxels, the diffusion tensor model works well. Non-Gaussian voxels are likely fibre crossings where the diffusion tensor fails. The algorithm requires you to select thresholds for separating the three classes. First run:

Black voxels are isotropic ("Order 0") or background, dark gray are anisotropic Gaussian ("Order 2") and light gray are non-Gaussian ("Order 4"). The light gray ones are likely fibre crossings. Scroll through the slices using the list on the right. Set the thresholds using the scroll bars or text fields at the bottom. At good threshold settings, the order 4 voxels should be fairly clustered, ie as few isolated order 4 voxels as possible with still a reasonable fraction of order 4 voxels overall. Around 5 to 10 percent order 4 voxels is typically about right.

Once you have chosen thresholds, make a note of them and generate the classified volume using voxelclassify with specified thresholds:

The output of this command has type int and every voxel contains either -1 for background or 0, 2 or 4 indicating the classification order. We can convert the results to NIfTI format and overlay them on the FA map using ITK SNAP:
voxelclassify -inputfile dwi.Bfloat -bgthresh 200 -schemefile 4Ddwi_b1000_bvector.scheme -order 4 -ftest 1.0E-09 1.0E-01 1.0E-03 > dwi_VC.Bint

Although the two-tensor model produces a reasonable output in multiple-fibre voxels, the fit in single-fibre voxels is misleading. For example, the genu of the corpus callosum contains fibre-orientation estimates that do not reflect the underlying fibre-orientation. The reconstruction can be improved by using voxelclassify to determine whether a single-tensor or two-tensor model should be fit on a voxel-by-voxel basis. To do this, we use the voxel classification generated in Voxel Classification section, above, in the multitenfit command:

sfplot can be used to generate images of the multi-tensor reconstructions, as opposed to just the principal orientations. To do this, we first need to split both multitensor.Bdouble and fa.img into slices. For this step we use the UNIX split command:

In the command above, the -projection 1 -2 flag is used to flip the y-direction for the ODFs, ensuring that they are displayed in the correct orientation. The -minmaxnorm flag normalises the ODF to exaggerate its peaks.

To view the image on a UNIX machine with Image Magick installed, use the command display -size 3584x3584 dwi_multitensor_slicear.rgb (some versions of ImageMagick may require the flag -depth 8 in addition to the -size 3584x3584 flag). The file extension directly affects how display reads the file. Alternatively, you can use the ImageMagick convert command to convert dwi_multitensor_slicear.rgb to a png or other common image formats. For those that are not using UNIX, the file can be loaded into matlab using the following command:

4. QBall

The first stage required for QBall reconstruction is the generation of the QBall reconstruction matrix. This matrix transforms the data from a sphere in q-space to the ODF. To calculate the QBall reconstruction matrix we use qballmx:

@qballmx -schemefile 4Ddwi_b1000_bvector.scheme > qballMatrix.Bdouble

This command uses sensible defaults for the various parameters of the QBall algorithm and is equivalent to

The extra options can be changed to adapt the behaviour of the algorithm, see [1,2] for details. If the program has run successfully, you will see information about the basis functions used as well as their settings on the command console. Record this information for future reference as you will need it later.

We can check that the settings (i.e. the number ODF basis functions and their widths) give reasonable results by generating some test ODFs using synthetic data from datasynth using the commands

We choose the value for -bgthresh by thresholding the b0 image such that the background is masked while brain voxels remain. Use any image viewer (for example MRIcro, ITKsnap, Matlab or FSLView) to compare foreground and background intensities in the b0 image and pick a value that lies between the two.

To create a colour-coded image using sfplot, you first need to split the data into slices. As before, we will use the UNIX split command:

Camino also implements the spherical harmonic (SH) QBall reconstructions of Descoteaux [3] and Mukherjee [4]. We will now create an SH representation of the QBall ODFs and use sfpeaks to find the peaks of the ODF in each voxel, which we will view in pdview. sfpeaks is a program that finds a wide range of features for any multiple-fibre reconstruction (such as PASMRI, QBall, Spherical Deconvolution). In particular, for any given voxel it outputs the number of peaks found, the mean of the function, the orientations and strengths of the peaks as well as the Hessian (or matrix of second partial-derivatives, which describes the curvature of the peak in two orthogonal directions) of each peak. This information is used in several Camino programs, including the tractography and statistics programs. To create a spherical harmonic representation of the ODF, you will need to calculate a new QBall matrix with the flag -basistype sh and then run linrecon. The full commands are:

In addition to displaying the principal eigenvectors of the diffusion tensor, pdview also allows you to view the output of sfpeaks. To do this, use the flag -inputmodel pds. The background map uses the trace of the Hessian (i.e. the sharpness) of the dominant peak of each ODF. Alternatively, you can specify another greyscale background map (for example, a fractional anisotropy image) using the flag -scalarfile [filename]. The color-coding is added automatically by pdview. The command is:

5. Reduced encoding PASMRI

Now lets run reduced encoding (RE) PASMRI to recover PAS functions in each voxel. The following variable, RE, specifies the number of basis functions in the PAS representation. By default, the number of basis functions is the same as the number of gradient directions. The idea of RE-PASMRI is that we can use fewer basis functions, which speeds up the estimation of the PAS. The smaller the number, the quicker the reconstruction, but the more crude the representation. Here we'll try 3 levels of reduced encoding. In practice, you'll need to experiment a bit to pick a setting that recovers crossing fibres reliably and has manageable computation time. RE=5 is probably too low, but is a good setting to use as a first run to get the whole pipeline working, as PASMRI reconstruction is fast with this low setting. RE=16 we often find is a good setting with a good trade-off between accuracy and computation time. Worth comparing with a higher setting such as RE=30 though to see whether any significant improvement occurs. To process the data, we first select a reduced encoding and run mesd:

In the above command, the -fastmesd flag is used to speed up the reconstruction by skipping some of the integrity checks. Once the data have been reconstructed, we use sfpeaks to find the peak orientations, which we use as our fibre-orientation estimates:

The images below show a pdview screenshot for each of the reduced-encoding reconstructions (from left to right: RE=5, RE=16, RE=30). If there are voxels in white-matter that appear to have no fibre-orientation estimates (these usually show up as black squares with a white dot in the centre), you probably have the thresholds set too high. Reducing the threshold and/or removing the -fastmesd flag can improve results, although removing the -fastmesd flag can significantly increase processing times.

In addition to pdview, you can also use sfplot to generate an image for each slice of the reconstruction. We show an example sfplot command in the section on exploiting multi-core processors, below.

5.1 Exploiting multi-core processors for PASMRI reconstruction

If you have a multi-core processor, splitting the dataset into chunks and processing the data on several cores can speed up the reconstruction. In this example I am going to run the reconstruction separately for each slice, which is useful for generating nice visualizations of the results using sfplot. This directory will contain the output for each slice

You'll need to adapt some of the commands above for your data (specifically -voxeldims, -bgthresh and -size). You can find the script StringZeroPad in camino/SGE. convert is an ImageMagick command, which you may need to install.

Once that's all done, you can combine the results into a single file ready for tractography

6. Multi-fibre stats

In this section, we use mfrstats to generate several useful of statistics from the output of sfpeaks. First, we generate some synthetic data, noting the test function and rotation angle used. Here we create data for 100 crossing-fibre voxels with SNR=20 using the command

The output of mfrstats can be used to calculate measures of accuracy and precision, such as direction concentration and angle bias (see [1,2] for details).

Another useful tool is consfrac, which calculates the consistency fraction. This program requires a ground-truth, which in this example is given by the settings used to generate the synthetic data (namely the test function and the rotation angle of the second tensor). For a set of reconstructions to be consistent, they should have the same number of fibre-orientation estimates as the test function and the peak orientations should be within a given tolerance of the true fibre-orientation. To run consfrac, we use the command