dipy logo

Site Navigation

NIPY Community

Previous topic

Reconstruct with Constant Solid Angle (QBall)

Next topic

Reconstruct with Generalized Q-Sampling Imaging

Reconstruct with Diffusion Spectrum ImagingΒΆ

We show how to apply Diffusion Spectrum Imaging (Wedeen et al. Science 2012) to diffusion MRI datasets of Cartesian keyhole diffusion gradients.

First import the necessary modules:

import nibabel as nib
from dipy.data import fetch_taiwan_ntu_dsi, read_taiwan_ntu_dsi, get_sphere
from dipy.align.aniso2iso import resample
from dipy.reconst.dsi import DiffusionSpectrumModel

Download and read the data for this tutorial.

fetch_taiwan_ntu_dsi()
img, gtab = read_taiwan_ntu_dsi()

img contains a nibabel Nifti1Image object (data) and gtab contains a GradientTable object (gradient information e.g. b-values). For example to read the b-values it is possible to write print(gtab.bvals).

Load the raw diffusion data and the affine.

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

data.shape (96, 96, 60, 203)

This dataset has anisotropic voxel sizes, therefore reslicing is necessary.

affine = img.get_affine()

Read the voxel size from the image header.

voxel_size = img.get_header().get_zooms()[:3]

Instantiate the Model and apply it to the data.

dsmodel = DiffusionSpectrumModel(gtab)

Lets just use one slice only from the data.

dataslice = data[:, :, data.shape[2] / 2]

dsfit = dsmodel.fit(dataslice)

Load an odf reconstruction sphere

sphere = get_sphere('symmetric724')

Calculate the ODFs with this specific sphere

ODF = dsfit.odf(sphere)

print('ODF.shape (%d, %d, %d)' % ODF.shape)

ODF.shape (96, 96, 724)

In a similar fashion it is possible to calculate the PDFs of all voxels in one call with the following way

PDF = dsfit.pdf()

print('PDF.shape (%d, %d, %d, %d, %d)' % PDF.shape)

PDF.shape (96, 96, 17, 17, 17)

We see that even for a single slice this PDF array is close to 345 MBytes so we really have to be careful with memory usage when use this function with a full dataset.

The simple solution is to generate/analyze the ODFs/PDFs by iterating through each voxel and not store them in memory if that is not necessary.

from dipy.core.ndindex import ndindex

for index in ndindex(dataslice.shape[:2]):
    pdf = dsmodel.fit(dataslice[index]).pdf()

If you really want to save the PDFs of a full dataset on the disc we recommend using memory maps (numpy.memmap) but still have in mind that if you do that for example for a dataset of volume size (96, 96, 60) you will need about 20 which can take less space when reasonable spheres (with < 1000 vertices) are.

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.