Tutorial: Motion Correction for Diffusion Weighted Images

This tutorial gives a general introduction of using function mbalign in Camino package to align the diffusion-weighted images within a single acquisition.

SYNOPSIS
mbalign [options]

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>.

And also we need to make some other information of input image data clear and specify a few other options

-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 sigma is approximate noise standard deviation. A suggested value is sqrt(E(S^2)/2), where S is the signal in background and E denotes expectation over an ROI.

A camino program, datastats, can work it out for you as well.

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

where S32_M100.Bshort is a mask file for ROI (we will discuss more about generating mask in Section 3).

then the screen output would be:
=============================================================================================
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 7.868863E04 1.294450E04 1.137739E02
5 7.893223E01 7.438512E03 1.208215E03 3.475939E01
6 7.813621E01 7.298962E03 1.193695E03 3.454989E01
7 8.115909E01 7.870362E03 1.283565E03 3.582688E01
8 7.870379E01 7.465566E03 1.271279E03 3.565500E01
9 8.104558E01 7.776110E03 1.207724E03 3.475232E01
10 8.053585E01 7.688916E03 1.202894E03 3.468276E01
11 8.019136E01 7.728569E03 1.297915E03 3.602659E01
12 7.751786E01 7.154714E03 1.145696E03 3.384814E01
13 8.196583E01 8.091765E03 1.373367E03 3.705897E01
14 7.725234E01 7.106027E03 1.138104E03 3.373580E01
15 8.016002E01 7.625693E03 1.200064E03 3.464194E01
==============================================================================================

So we can choose 7.5 as the value of E(S^2). Then sigma, sqrt(E(S^2)/2), would be 1.9.

1.4.Make slice to check input volume

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 camino function shredder to extract one slice from each of the 32 3D components, and build a 3D image shown in Fig1:

Fig.1 An image volume containing 32 slices from 64 diffusion components.

To make the command easy to read, we can use variables representing data information:

COMPONENT=64

DATADIM_X=128
DATADIM_Y=128
DATADIM_Z=32

Since our dataset used in the example is float, so

TYPESIZE=4

If we would like to extract the middle slice along z direction, we can set

OFFSET=$(($DATADIM_X*$DATADIM_Y*$((DATADIM_Z/2))*$TYPESIZE))

Then, calling shredder

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 (*.hdr) to make *.img file readable by many visualisation softwares:

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

Using some visualisation tools, such like MRIcro(Fig.2) and camino, we can check the slice motion of 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 create 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, computer need to have 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, but more easily, once we have camino and FSL installed, default value of variable DIR_FSL can be set which can be easily found inside mbalign.

During the running of mbalign, there could be a lot of warning messages showing inside terminal window. Normally, just do not worry about it too much. Once the registration has got done, 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 program finishes, the temporary directory would be deleted (Removing junk files...), but we can still use -keepjunk to make it kept, which could be useful to help us to analyse the final output. Certainly, if the program is interrupted, this temporary fold will remain as junk files in computer.

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, such like

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 involve more measurement 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 a mask file can improve the quality of the reference images used in mbalign registration. And the data type of mask file should be "short".

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 the diffusion tensors.

-searchrange <angle>
Sometimes, the pitch of histogram will cause the failure of registration. Simply narrow down the angle search range can cover this problem for 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 indicate 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 UpdateGradient.m 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 most often to use.

# 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

Acknowledgement

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

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.