Co-registration of Diffusion-Weighted Images for Motion Correction

This tutorial introduces the function mbalign in Camino, which is used to align the diffusion-weighted images within a single acquisition.

HEALTH WARNING: We no longer directly support this module. It also relies on an old version of FSL. The tutorial remains because some users still find mbalign useful. The technique can be quite effective, but the code and documentation is a bit flaky, it can be hard to tune for good performance, and the original author has left the Camino team. So use at your own risk!

1.  Preparation before running mbalign

1.1  General data files and information

Before running mbalign, we need to have input and scheme files ready, and use:

-inputfile <Input voxel-order file>
-schemefile <Scheme file name>

to specify in command options. If the input file is in scanner-order, we can use:

-scanner -inputfile <Input scanner-order file>.

In addition we need to specify the size of the images, the voxel size and the standard deviation of the noise in the data:

-datadims X Y Z <Number of voxels in each dimension>
-voxeldims x y z <Voxel size in mm>
-sigma <Standard deviation of noise>

1.2  Sigma

The value of sigma is the approximate standard deviation of the noise. An estimate of sigma is sqrt(E(S^2)/2), where S is the signal in background and E denotes expectation over an ROI. The Camino program, datastats, can work it out for you:

datastats -schemefile S.scheme -bgmask S32_BG.Bshort -inputfile S32.Bfloat

where S32_BG.Bshort is a binary image containing an ROI in a background region of the image.

The output looks like this: 
Foreground voxel count: 142584 
Component E(S) E(S^2) Var(S) Std(S) 
1 2.882503E02 9.834995E04 1.526169E04 1.235382E02 
2 3.010773E02 1.070645E05 1.641700E04 1.281288E02 
3 2.876059E02 9.854322E04 1.582610E04 1.258018E02 
4 2.564062E02 9.868863E04 1.294450E04 1.137739E02 
5 2.017598E02 9.636593E04 1.028376E04 1.293741E02 
6 2.915845E02 9.591012E04 1.817574E04 1.197654E02 

A well chosen background region that genuinely contains no signal should show similar statistics in the output above in both non-diffusion-weighted (b=0) images and diffusion weighted images. The first four images in the example above have b=0 and the subsequent ones have b>0. We see little difference in the statistics, which suggests the background genuinely contains no signal.

The value of E(S^2) is around 1E5, so we might pick sigma = sqrt(E(S^2)/2), ie around 224.

However, you may need to play with the setting of sigma to get the best results out of mbalign. The value is somewhat artificial, because really we just seek a good threshold that reliably rejects measurements from corrupted images stopping them contributing to the diffusion tensor fit. Several users have reported that they have to set sigma 100 or 1000 times smaller than the estimate of sigma described above to get good results. Others report having to set sigma much higher, so this seems to be data dependent.

1.3  Make slice to check input volume

This section shows a simple way to check that the input volume is what the program expects and also to give you a feel for how much realignment is required. It creates an image volume containing the corresponding slice from each image volume in the data set.

If the image file is in voxel-order, we need to transfer it to scanner-order first:

voxel2scanner -voxels $((128*128*32)) -inputdatatype float -outputdatatype float -components 64 -inputfile S32.Bfloat > S32.scan.Bfloat

Then we can use the Camino function shredder to extract one slice from each of the 64 3D components, and build a 3D image shown in Fig1:

Fig.1 An image volume containing 1 slice from each of the 64 component images in the data set.

Here is an example script for creating the image volume above:

COMPONENT=64 
DATADIM_X=128
DATADIM_Y=128
DATADIM_Z=32

#Since our dataset used in the example is float, so
TYPESIZE=4

# Extract the middle slice
OFFSET=$(($DATADIM_X*$DATADIM_Y*$((DATADIM_Z/2))*$TYPESIZE))
shredder $OFFSET $(($DATADIM_X*$DATADIM_Y*$TYPESIZE)) $(($DATADIM_X*$DATADIM_Y*$(($DATADIM_Z-1))*$TYPESIZE)) < S32.scan.Bfloat > S32._SLICECHECK.img

# Make a header file
VOXELDIM_X=1.88
VOXELDIM_Y=1.88
VOXELDIM_Z=2.0

analyzeheader -voxeldims $VOXELDIM_X $VOXELDIM_Y $VOXELDIM_Z -datadims $DATADIM_X $DATADIM_Y $COMPONENT -datatype float > S32._SLICECHECK.hdr

Now we can use visualisation tools, like MRIcro (Fig.2), to check the alignment of the input data set.

Fig.2 Using MRIcro to see slice volume

2.  Run mbalign

2.1  Run mbalign simply

An example of using mbalign in the simplest way is:

mbalign -datatype float -schemefile S.scheme -datadims 128 128 60 -voxeldims 1.88 1.88 2.0 -sigma 1.9 -inputfile S32.Bfloat 

then the screen output would be:

:-) I have everything I need!
WARNING: No -bgmask and -bgthresh input, using zero background threshold. Performance may improve with better threshold.
WARNING: No -components input. Using 64 components from schemefile. If wrong, restart with -components option.
64 components inside inputfile.
No temp directory specified. Trying to use /tmp/S32_17.03.08_180229.
Successfully created /tmp/S32_17.03.08_180229.
Linux system detected.
No output file name specified. Output file will be: /tmp/S32.out.Bfloat 
Disk space temporarily used during calculation is about 2264924160. Make sure space is available!

The program will create a temporary directory used for calculation process. It can be either specified by -tmpdir or automatically created according to the file name and current time. We do not have to use -outputfile to specify the output file name, and the program can generate it according to the input file name, which is like the example shown above. To run mbalign, your computer needs to have the registration software flirt(http://www.fmrib.ox.ac.uk/fsl/flirt/index.html) installed, which is part of FSL library(http://www.fmrib.ox.ac.uk/fsl/). We also need to specify the flirt direction by using -fsldir. However an easier way once we have Camino and FSL installed is to set the default value of the variable DIR_FSL inside @~/camino/bin/mbalign@. During the running of mbalign, there could be a lot of warning messages showing inside terminal window. Normally, just don't worry about it too much. Once the registration is complete, we may see the message shown below.

... 
Registering ori.ScannerOrder.Bfloat.ck by using ref.ScannerOrder.Bfloat.ck
Registering ori.ScannerOrder.Bfloat.cl by using ref.ScannerOrder.Bfloat.cl
Transferring output to Big-endian format... 
Making img and hdr files for slice checking...
Transfer output file to voxel-order... 
Removing junk files... 
Scheme file not updated.
Aligned data set output to /tmp/S32.out.Bfloat.
Program finished at 
Mon Mar 24 23:36:24 GMT 2008

After the program finishes, the temporary directory would be deleted (Removing junk files...), but we can still use -keepjunk to keep it. This could be useful to help us to analyse the final output. Please note if the program is interrupted, this temporary folder will remain as junk files on the computer and need to be manually deleted. If we did not specify -outputfile, the program can generate it according to the input file name (Aligned data set output to /tmp/S32.out.Bfloat). We will discuss scheme file updating in Section 4.Update Gradient.

2.2  Run mbalign in an advanced way

Some other options mignt need to use: -flirtsearchcost <Search cost function used in flirt>
Default cost function is mutualinfo (Mutual Information). Other options are corratio,normcorr,normmi,leastsq. -flirttransform <Transformation used in flirt>
Default transformation is affine. The other option is rigid. -searchrange <angle>
Default is 90, which means search range is between -90 and 90 in all x, y and z directions. -eddy
Specifies registration for eddy-current induced distortion. -datatype <Data type for input and output files>
Default is float. -scanout <output scanner-order file>
Adds an extra output file in scanner-order. This won't stop default voxel-order output. -omat <File name>
Output transform matrix in ascii format. -slicecheck <File name>
Output a pair of <File name>.img and <File name>.hrd files. Default is no calculation.
When all the options are decided, we can run mbalign in an advanced way e.g. mbalign -datatype float -schemefile S.scheme -datadims 128 128 60 -voxeldims 1.88 1.88 2.0 -sigma 1.9 -fsldir /cs/research/medim/common0/green/common/fsl/fslRH9/ -inputfile S32.Bfloat -slicecheck S32.fmr.slice.check -outputfile S32.fmr.Bfloat -omat S32.fmr.mat.txt -scanout S32.fmr.scanout.Bfloat -keepjunk -tmpdir tmp.fmr

3.  Improve performance

There are a few options can be used to improve the performance of mbalign. -sigma
High sigma allows the program to include more measurements in the DT fitting, and low sigma leads rejections during the DT fitting. Based on this theory, we can change the value of sigma to improve the reference making. -bgmask <Mask file>'
Use of a mask file can improve the quality of the reference images used in mbalign registration. The data type of mask file should be "short". The Camino function mask can help to create a background mask from a voxel-ordered DW data file by thresholding the average b=0 measurement. mask -inputfile S32.Bfloat -inputdatatype float -schemefile S.scheme -bgthresh 100 -outputdatatype short > S32_M100.Bshort

Fig.3 View projection of FA map generated by Camino

Fig.4 View projection of mask file generated by Camino

We can also use Matlab to make a mask file. -bgthresh <Background threshold>
Decide the value of threshold, and improve fitting of the diffusion tensors. -searchrange <angle>
Sometimes, the pitch of the histogram will cause the failure of registration. Simply narrowing down the angle search range can solve this problem most of the time.

4.  Update Gradient

Updating diffusion gradients after registration is not an essential procedure, but can improve the registration result. Fig.5 explains the reason why diffusion gradients need to be updated after registration.

Fig.5 Illustration of how rotation affects the effective gradient direction. The arrows indicates gradient directions. (a) Head without rotation. Suppose the mouth is a fibre. The signal is high because the gradient is perpendicular to the fibre. (b) Head with rotation. The signal is lower in the mouth fibre. (c) The unrotated head (after registration). The effective gradient direction for the corrected image is rotated. We can use our Matlab function to update diffusion gradient for registered data set, and the usage method is explained inside the matlab file. Since the transformations used in the registration contribute to the gradient updating, we MUST save the transform matrix when running mbalign, using -omat.

5.  More hints

Quite a few default options can be change in the source code of mbalign. Make the default values to the ones most frequently used can make the everyday use of mbalign much simpler. Inside mbalign source code, we can find and change the default options from the following part.

####################################################
##### Change default variables to match system ##### 
####################################################
# Hint: To make your input arguments simple, set
# default input which you use most often
# FSL directory
DIR_FSL=/cs/research/medim/common0/green/common/fsl/fslRH9
# LIM_ROTATE is default for -searchrange
LIM_ROTATE=90
# Available cost functions are:
# mutualinfo corratio,normcorr,normmi,leastsq.
SEARCH_COST=mutualinfo
#Degree of freedom
# 12 for affine; 6 for rigid.
DOF=12

6.  Acknowledgement

We would like to thank Geoff Parker and Karl Embleton - University of Manchester, for providing the brain data.

7.  Reference

Y. Bai and D. C. Alexander, “Model-Based Registration to Correct for Motion between Acquisitions in Diffusion MR Imaging”, The Fifth IEEE International Symposium on Biomedical Imaging ( ISBI 2008), May 2008.