Multi Fiber Tractography

1.  Overview

Note: Parts of this tutorial are out of date. OOGL output is no longer supported, use vtkstreamlines instead and render in a VTK viewer such as Paraview.

This tutorial runs through deterministic and probabilistic tractography using multiple-fiber reconstructions like PASMRI and q-ball. I'll start with a data file of diffusion weighted measurements in Camino voxel-order format, called B4069.Bfloat, and corresponding schemefile, called B4069.scheme. The data set has 105 gradient directions with b value 4069s/mm^2 (it actually comes from a fixed brain) and 25 b=0 images. Here are the stages we'll go through: Deterministic and probabilistic DTI tractography as a baseline PAS reconstruction using reduced encoding PASMRI (see Sweet and Alexander ISMRM 2010). Deterministic tractography using PASMRI Probabilistic tractography using PASMRI q-ball reconstruction and tractography for comparison Define some variables used throughout, which are specific to my data set:

export DSID = 'B4069'
export SCHEMEFILE = 'B4069.scheme1'
export XSIZE = 128
export YSIZE = 141
export YSIZE = 30
export COMPONENTS=130
export SEEDFILEROOT = 'MidSagCC'

2.  DTI Tractography

This is covered in detail in other tutorials, so I'll keep it brief. Fit the diffusion tensor and check that the tensor orientations correspond with the image orientation:

dtfit ${DSID}.Bfloat $SCHEMEFILE > /tmp/${DSID}_DT.Bdouble
dteig < /tmp/${DSID}_DT.Bdouble | pdview -datadims $XSIZE $YSIZE $ZSIZE

Run deterministic tractography from a seed region (mine is the midsagittal corpus callosum):

track -inputmodel dt -seedfile $SEEDFILEROOT -outputroot ${DSID}_${SEEDFILEROOT}_DT -outputtracts oogl < /tmp/${DSID}_DT.Bdouble

Display the tract:

geomview ${DSID}_${SEEDFILEROOT}_DT1.oogl &

Here's what we get

Now run some probabilistic tractography. Start by generating the look up table and create the pdfs in each voxel

dtlutgen -schemefile $SCHEMEFILE -trace 8E-10 -snr 16.0 -inversion 1 > ${DSID}_DT_LUT.Bdouble
cat /tmp/${DSID}_DT.Bdouble | picopdfs -inputmodel dt -luts ${DSID}_DT_LUT.Bdouble > ${DSID}_DT_PDFs.Bdouble

Now do the tracking (you may want to increase the number of iterations in a real application):

track -inputmodel pico -seedfile $SEEDFILEROOT -iterations 10 -outputroot ${DSID}_${SEEDFILEROOT}_DT10 -outputtracts oogl < ${DSID}_DT_PDFs.Bdouble

And display the results for comparison:

geomview ${DSID}_${SEEDFILEROOT}_DT101.oogl &

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

#export RE=30; export REID=30;
#export RE=16; export REID=16;
export RE=5; export REID=05;

I'm 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 mkdir PAS${REID}_${DSID}@] If you only have one processor, use the following for loop (or just run the whole file in one go):

for ((i=1; i<=${ZSIZE}; i=i+1)); do

I happen to have 8, so I ran a separate process on each using the following set of 8 different for loops over the slices:

for ((i=1; i<=${ZSIZE}; i=i+8)); do
#for ((i=2; i<=${ZSIZE}; i=i+8)); do
#for ((i=3; i<=${ZSIZE}; i=i+8)); do
#for ((i=4; i<=${ZSIZE}; i=i+8)); do
#for ((i=5; i<=${ZSIZE}; i=i+8)); do
#for ((i=6; i<=${ZSIZE}; i=i+8)); do
#for ((i=7; i<=${ZSIZE}; i=i+8)); do
#for ((i=8; i<=${ZSIZE}; i=i+8)); do
echo $i
export SNUM=`StringZeroPad $i 2`
shredder $((${COMPONENTS}*${XSIZE}*${YSIZE}*(i-1)*4)) $((${COMPONENTS}*${XSIZE}*${YSIZE}*4)) $((${COMPONENTS}*${XSIZE}*${YSIZE}*${ZSIZE}*4)) < ${DSID}.Bfloat > /tmp/${DSID}_Slice${SNUM}.Bfloat
dtfit /tmp/${DSID}_Slice${SNUM}.Bfloat $SCHEMEFILE | fa -outputdatatype float > /tmp/${DSID}_Slice${SNUM}_FA.img
analyzeheader -datadims ${XSIZE} ${YSIZE} 1 -voxeldims 0.5 0.5 0.5 -datatype float > /tmp/${DSID}_Slice${SNUM}_FA.hdr
mesd -schemefile $SCHEMEFILE -filter PAS 1.4 -fastmesd -mepointset $RE -bgthresh 200 < /tmp/${DSID}_Slice${SNUM}.Bfloat > PAS${REID}_${DSID}/${DSID}_Slice${SNUM}_PAS${REID}.Bdouble
sfpeaks -schemefile $SCHEMEFILE -inputmodel maxent -filter PAS 1.4 -mepointset $RE -inputdatatype double < PAS${REID}_${DSID}/${DSID}_Slice${SNUM}_PAS${REID}.Bdouble > PAS${REID}_${DSID}/${DSID}_Slice${SNUM}_PAS${REID}_PDs.Bdouble
sfplot -xsize ${YSIZE} -ysize ${XSIZE} -inputmodel maxent -filter PAS 1.4 -mepointset $RE -pointset 2 -minifigsize 30 30 -minifigseparation 2 2 -inputdatatype double -dircolcode -projection 1 2 < PAS${REID}_${DSID}/${DSID}_Slice${SNUM}_PAS${REID}.Bdouble > /tmp/${DSID}_Slice${SNUM}_PAS${REID}.rgb
convert -size 4096x4512 -depth 8 /tmp/${DSID}_Slice${SNUM}_PAS${REID}.rgb PAS${REID}_${DSID}/${DSID}_Slice${SNUM}_PAS${REID}.png
rm /tmp/${DSID}_Slice${SNUM}_PAS${REID}.rgb /tmp/${DSID}_Slice${SNUM}.Bfloat /tmp/${DSID}_Slice${SNUM}_FA.hdr /tmp/${DSID}_Slice${SNUM}_FA.img

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

cat PAS${REID}_${DSID}/*PAS${REID}_PDs.Bdouble > PAS${REID}_${DSID}/PDs.Bdouble

Now you can check that principal directions of the PAS functions:

pdview -inputmodel pds -datadims ${XSIZE} ${YSIZE} ${ZSIZE} -scalarfile ${DSID}_FA.img < PAS${REID}_${DSID}/PDs.Bdouble

4.  PAS Tractography

Deterministic tractography:

track -inputmodel pds -seedfile $SEEDFILEROOT -outputroot ${DSID}_${SEEDFILEROOT}_PAS${REID} -outputtracts oogl < PAS${REID}_${DSID}/PDs.Bdouble

For comparison with DT results:

geomview ${DSID}_${SEEDFILEROOT}_PAS${REID}1.oogl &

Here are results for RE=5, 16 and 30:

Now try some probabilistic tractography. Start by calibrating the pdf generator. This command creates a full size calibration data set (see sfpicocalibdata man page):

sfpicocalibdata -schemefile $SCHEMEFILE -snr 12 -infooutputfile ${DSID} > ${DSID}_CalibData.Bfloat

For my postmortem data set, I had to add -trace 8E-10, to specify a lower diffusivity for fixed tissue at room temperature than the default, which is meant for in-vivo tissue. I've removed it from the above so in-vivo users can cut and paste. I've used -snr 12, as the signal to noise of my unweighted images is around 12. Performance is reasonably robust to this setting so don't worry if the SNR varies over the image. Just pick a ball park correct figure. Now do the actual calibration. This involves running RE-PAS on the calibration data set and constructing a look-up table. Over the full calibration data set, even RE-PAS can take a while. To get a quick set of results for a first pass, have a look at the sfpicocalibdata man page to see how you can reduce the size of the data set. You may find that you need to add arguments like -binincsize 2 -minvectsperbin 10 to sflutgen below to get it to work on smaller calibration data sets; check the sflutgen man page to see what they do. If you find you need such arguments for the full calibration data set, the calibration data set is not large enough. Anyway, here goes...

mesd -schemefile $SCHEMEFILE -filter PAS 1.4 -fastmesd -mepointset $RE < ${DSID}_CalibData.Bfloat > ${DSID}_CalibDataPAS${REID}.Bdouble
sfpeaks -schemefile $SCHEMEFILE -inputmodel maxent -filter PAS 1.4 -mepointset $RE -inputdatatype double < ${DSID}_CalibDataPAS${REID}.Bdouble > ${DSID}_CalibDataPAS${REID}_PDs.Bdouble
sflutgen -infofile ${DSID} -outputstem ${DSID}_PAS${REID}_LUT -pdf watson < ${DSID}_CalibDataPAS${REID}_PDs.Bdouble

Note that the above uses a Watson model for the pdfs. You can also use Bingham; see sflutgen man page.

Now we can get the pdfs and do the tracking:

picopdfs -inputmodel pds -numpds 3 -pdf watson -luts ${DSID}_PAS${REID}_LUT_oneFibreLineCoeffs.Bdouble ${DSID}_PAS${REID}_LUT_twoFibreLineCoeffs.Bdouble ${DSID}_PAS${REID}_LUT_twoFibreLineCoeffs.Bdouble < PAS${REID}_${DSID}/PDs.Bdouble > PAS${REID}_${DSID}/PDsWatsonPDFs.Bdouble
track -inputmodel pico -pdf watson -numpds 3 -iterations 10 -seedfile $SEEDFILEROOT -outputroot ${DSID}_${SEEDFILEROOT}_PAS${REID}10 -outputtracts oogl < PAS${REID}_${DSID}/PDsWatsonPDFs.Bdouble
geomview ${DSID}_${SEEDFILEROOT}_PAS${REID}101.oogl &

Results for R=5, 16 and 30 as before:

5.  Q-ball reconstruction and tractography

Do the q-ball reconstruction using spherical harmonic representation up to order 6:

qballmx -schemefile $SCHEMEFILE -basistype sh -order 6 > ${DSID}_QBALLMX.Bdouble
linrecon $DSID.Bfloat $SCHEMEFILE ${DSID}_QBALLMX.Bdouble -normalize -bgthresh 200 > ${DSID}_QBall.Bdouble
sfpeaks -inputmodel sh -order 6 -density 100 -searchradius 1.0 < ${DSID}_QBall.Bdouble > ${DSID}_QBallPDs.Bdouble

Note the sfpeaks command can take quite a while.

Check the directions are OK:

pdview -inputmodel pds -datadims ${XSIZE} ${YSIZE} ${ZSIZE} -scalarfile ${DSID}_FA.img < ${DSID}_QBallPDs.Bdouble

Do some deterministic tractography:

track -inputmodel pds -seedfile $SEEDFILEROOT -outputroot ${DSID}_${SEEDFILEROOT}_QBall -outputtracts oogl < ${DSID}_QBallPDs.Bdouble
geomview ${DSID}_${SEEDFILEROOT}_QBall1.oogl &

Using the same calibration data set constructed above for PAS-PICo, calibrate for q-ball PICo

linrecon ${DSID}_CalibData.Bfloat $SCHEMEFILE ${DSID}_QBALLMX.Bdouble -normalize > ${DSID}_CalibDataQBall6.Bdouble
sfpeaks -inputmodel sh -order 6 -inputdatatype double -density 100 -searchradius 1.0 < ${DSID}_CalibDataQBall6.Bdouble > ${DSID}_CalibDataQBall6_PDs.Bdouble
sflutgen -infofile ${DSID} -outputstem ${DSID}_QBall6_LUT -pdf watson < ${DSID}_CalibDataQBall6_PDs.Bdouble

Get the pdfs and run the tractography:

picopdfs -inputmodel pds -numpds 3 -pdf watson -luts ${DSID}_PAS${REID}_LUT_oneFibreLineCoeffs.Bdouble ${DSID}_PAS${REID}_LUT_twoFibreLineCoeffs.Bdouble ${DSID}_PAS${REID}_LUT_twoFibreLineCoeffs.Bdouble < ${DSID}_QBallPDs.Bdouble > ${DSID}_QBallPDsWatsonPDFs.Bdouble
track -inputmodel pico -pdf watson -numpds 3 -iterations 10 -seedfile $SEEDFILEROOT -outputroot ${DSID}_${SEEDFILEROOT}_QBall10 -outputtracts oogl < ${DSID}_QBallPDsWatsonPDFs.Bdouble
geomview ${DSID}_QBallPDsWatsonPDFs.Bdouble &