Reconstruct with Generalized Q-Sampling Imaging#

We show how to apply Generalized Q-Sampling Imaging [Yeh2010] to diffusion MRI datasets. You can think of GQI as an analytical version of DSI orientation distribution function (ODF) (Garyfallidis, PhD thesis, 2012).

First import the necessary modules:

import numpy as np
from dipy.core.gradients import gradient_table
from import get_fnames, get_sphere
from import read_bvals_bvecs
from import load_nifti
from dipy.reconst.gqi import GeneralizedQSamplingModel
from dipy.direction import peaks_from_model

Download and get the data filenames for this tutorial.

fraw, fbval, fbvec = get_fnames('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:


Load the raw diffusion data and the affine.

data, affine, voxel_size = load_nifti(fraw, return_voxsize=True)
bvals, bvecs = read_bvals_bvecs(fbval, fbvec)
bvecs[1:] = (bvecs[1:] /
                 np.sqrt(np.sum(bvecs[1:] * bvecs[1:], axis=1))[:, None])
gtab = gradient_table(bvals, bvecs)
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.

Instantiate the model and apply it to the data.

gqmodel = GeneralizedQSamplingModel(gtab, sampling_length=3)

The parameter sampling_length is used here to

Lets just use one slice only from the data.

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

mask = dataslice[..., 0] > 50

gqfit =, mask=mask)
  0%|                                                                                                                                          | 0/4733 [00:00<?, ?it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4733/4733 [00:00<00:00, 424379.85it/s]

Load an ODF reconstruction sphere

sphere = get_sphere('repulsion724')

Calculate the ODFs with this specific sphere

ODF = gqfit.odf(sphere)

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

Using peaks_from_model we can find the main peaks of the ODFs and other properties.

gqpeaks = peaks_from_model(model=gqmodel,

gqpeak_values = gqpeaks.peak_values

gqpeak_indices show which sphere points have the maximum values.

gqpeak_indices = gqpeaks.peak_indices

It is also possible to calculate GFA.

GFA = gqpeaks.gfa

print('GFA.shape (%d, %d)' % GFA.shape)
GFA.shape (96, 96)

With parameter return_odf=True we can obtain the ODF using gqpeaks.ODF

gqpeaks = peaks_from_model(model=gqmodel,

This ODF will be of course identical to the ODF calculated above as long as the same data and mask are used.

print(np.sum(gqpeaks.odf != ODF) == 0)

The advantage of using peaks_from_model is that it calculates the ODF only once and saves it or deletes if it is not necessary to keep.



Yeh, F-C et al., Generalized Q-sampling imaging, IEEE Transactions on Medical Imaging, vol 29, no 9, 2010.

Total running time of the script: (0 minutes 5.341 seconds)

Gallery generated by Sphinx-Gallery