dipy logo

Site Navigation

NIPY Community

Table Of Contents

Previous topic

A note about the examples

Next topic

Find correspondence between tractographies

From raw data to tractographies

Overview

This example gives a tour of some simple features of dipy.

First import the necessary modules

numpy is for numerical computation

import numpy as np

nibabel is for data formats

import nibabel as nib

dipy.reconst is for the reconstruction algorithms which we use to create directionality models for a voxel from the raw data.

import dipy.reconst.gqi as gqi
import dipy.reconst.dti as dti
dipy.tracking is for tractography algorithms which create sets of tracks by integrating
directionality models across voxels.
from dipy.tracking.propagation import EuDX

dipy.data is for small datasets we use in tests and examples.

from dipy.data import get_data

Isotropic voxel sizes required

dipy requires its datasets to have isotropic voxel size. If you have datasets with anisotropic voxel size then you need to resample with isotropic voxel size. We have provided an algorithm for this. You can have a look at the example resample_aniso_2_iso.py

Accessing the necessary datasets

get_data provides data for a small region of interest from a real diffusion weighted MR dataset acquired with 102 gradients (including one for b=0).

In order to make this work with your data you should comment out the line below and add the paths for your nifti file (*.nii or *.nii.gz) and your *.bvec and *.bval files.

If you are not using nifti files or you don’t know how to create the *.bvec and *.bval files from your raw dicom (*.dcm) data then you can either try recent module nibabel.nicom

try:
    from nibabel.nicom.dicomreaders import read_mosaic_dir
except:
    print('nicom for dicom is not installed')

or to convert the dicom files to nii, bvec and bval files using dcm2nii.

fimg,fbvals,fbvecs=get_data('small_101D')

Load the nifti file found at path fimg as an Nifti1Image.

img=nib.load(fimg)

Read the datasets from the Nifti1Image.

data=img.get_data()
print('data.shape (%d,%d,%d,%d)' % data.shape)

This produces the output:

data.shape (6,10,10,102)

As you would expect, the raw diffusion weighted MR data is 4-dimensional as we have one 3-d volume (6 by 10 by 10) for each gradient direction.

Read the affine matrix
which gives the mapping between volume indices (voxel coordinates) and world coordinates.
affine=img.get_affine()

Read the b-values which are a function of the strength, duration, temporal spacing and timing parameters of the specific paradigm used in the scanner, one per gradient direction.

bvals=np.loadtxt(fbvals)

Read the b-vectors, the unit gradient directions.

gradients=np.loadtxt(fbvecs).T

Calculating models and parameters of directionality

We are now set up with all the data and parameters to start calculating directional models for voxels and their associated parameters, e.g. anisotropy.

Calculate the Single Tensor Model (STM).

ten=dti.Tensor(data,bvals,gradients,thresh=50)

Calculate Fractional Anisotropy (FA) from STM

FA=ten.fa()
print('FA.shape (%d,%d,%d)' % FA.shape)

As expected the FA is a 3-d array with one value per voxel:

FA.shape (6,10,10)

Generate a tractography

Here we use the Euler Delta Crossings (EuDX) algorithm. The main input parameters of EuDX are

  • an anisotropic scalar metric e.g. FA
  • the indices for the peaks on the sampling sphere.

Other important options are

  • the number of random seeds where the track propagation is initiated,
  • a stopping criterion, for example a low threshold for anisotropy. For instance if we are using Fractional Anisotropy (FA) a typical threshold value might be a_low=.2
eu=EuDX(a=FA,ind=ten.ind(),seeds=10000,a_low=.2)

EuDX returns a generator class which yields a further track each time this class is called. In this way we can generate millions of tracks without using a substantial amount of memory. For an example of what to do when you want to generate millions of tracks with minimum memory usage have a look at save_dpy.py in the examples directory. However, in the current example that we only have 10000 seeds, and we can load all tracks in a list using list comprehension([]) without having to worry about memory.

ten_tracks=[track for track in eu]

In dipy we usually represent tractography as a list of tracks. Every track is a numpy array of shape (N,3) where N is the number of points in the track.

print ('The number of FA tracks is %d' % len(ten_tracks))
print ('The number of points in ten_tracks[130] is %d' % len(ten_tracks[130]))
print ('The points in ten_tracks[130] are:')
print ten_tracks[130]

As we use random seeding for the tractography the results will differ when repeated, however one run gave us the following information:

The number of FA tracks is 8280
The number of points in ten_track[130] is 7
The points in ten_tracks[130] are:
[[ 1.73680878  5.08249903  4.48492956]
 [ 1.45797026  4.76981783  4.21201992]
 [ 1.14244306  4.46308756  3.97461915]
 [ 0.84001541  4.14648438  3.73316503]
 [ 0.53758776  3.82988143  3.49171114]
 [ 0.22055824  3.52935386  3.24845099]
 [ 0.22055824  3.52935386  3.24845099]]

Another way to represent tractography is as a numpy array of numpy objects. This way has an additional advantage that it can be saved very easily using numpy utilities. In theory, in a list it is faster to append an element, and in an array is faster to access. In other words both representations have different pros and cons. Other representations are possible too e.g. graphtheoretic etc.

ten_tracks_asobj=np.array(ten_tracks,dtype=np.object)
np.save('ten_tracks.npy',ten_tracks_asobj)
print('FA tracks saved in ten_tracks.npy')

Crossings and Generalized Q-Sampling

You probably have heard about the problem of crossings in diffusion MRI. The single tensor model cannot detect a simple crossing of two fibres. However with Generalized Q-Sampling (GQS) this is possible even up to a quadruple crossing or higher depending on the resolution of your datasets. Resolution will typically depend on signal-to-noise ratio and voxel-size.

gqs=gqi.GeneralizedQSampling(data,bvals,gradients)

A useful metric derived from GQS is Quantitative Anisotropy (QA).

QA=gqs.qa()
print('QA.shape (%d,%d,%d,%d)' % QA.shape)

QA is a 4-d array with up to 5 peak QA values for each voxel:

QA.shape (6,10,10,5)

QA array is significantly different in shape from the FA array, however it too can be directly input to the EuDX class:

eu2=EuDX(a=QA,ind=gqs.ind(),seeds=10000,a_low=.0239)
This shows one of the advantages of our EuDX algorithm: it can be used with a wide range of model-based methods, such as
  • Single Tensor
  • Multiple Tensor
  • Stick & Ball
  • Higher Order Tensor
and model-free methods such as
  • DSI
  • QBall
  • GQI etc.

We designed the algorithm this way so we that we can compare directly tractographies generated from the same dataset with very different models and/or choices of threshold.

Now we look at the QA tractography:

gqs_tracks=[track for track in eu2]
print('The number of QA tracks is %d' % len(gqs_tracks))

with output:

The number of QA tracks is 14022

Note the difference between the number of gqs_tracks and ten_tracks. There are more with QA than with FA. This is because of the presence of crossings which GQI can detect but STM cannot. When the underlying directionality model supports crossings then distinct tracks will be propagated from a seed towards the different directions in equal abundance.

In dipy it is very easy to count the number of crossings in a voxel, volume or region of interest

gqs_tracks_asobj=np.array(gqs_tracks,dtype=np.object)
np.save('gqs_tracks.npy',gqs_tracks_asobj)
print('QA tracks saved in gqs_tracks.npy')

This is the end of this very simple example You can reload the saved tracks using np.load from your current directory. You can optionaly install python-vtk and visualize the tracks using fvtk:

from dipy.viz import fvtk
r=fvtk.ren()
fvtk.add(r,fvtk.line(ten_tracks,fvtk.red,opacity=0.05))
gqs_tracks2=[t+np.array([10,0,0]) for t in gqs_tracks]
fvtk.add(r,fvtk.line(gqs_tracks2,fvtk.green,opacity=0.05))

Press ‘s’ to save this screenshot when you have displayed it with fvtk.show. Or you can even record a video using fvtk.record.

You would show the figure with something like:

fvtk.show(r,png_magnify=1,size=(600,600))

To record a video of 50 frames of png, something like:

fvtk.record(r,cam_pos=(0,40,-40),cam_focal=(5,0,0),n_frames=50,magnification=1,out_path='nii_2_tracks',size=(600,600),bgr_color=(0,0,0))
../_images/nii_2_tracks1000000.png

Same region of interest with different underlying voxel representations generates different tractographies.

# Here's how we make the figure.
print('Saving illustration as nii_2_tracks1000000.png')
fvtk.record(r,n_frames=1,out_path='nii_2_tracks',size=(600,600))

Example source code

You can download the full source code of this example. This same script is also included in the dipy source distribution under the doc/examples/ directory.