.. AUTO-GENERATED FILE -- DO NOT EDIT!

.. _example_reconst_gqi:


===============================================
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.data import fetch_taiwan_ntu_dsi, read_taiwan_ntu_dsi, get_sphere
  from dipy.reconst.gqi import GeneralizedQSamplingModel
  from dipy.direction import peaks_from_model
  

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

Read the voxel size from the image header.

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

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 = gqmodel.fit(dataslice, mask=mask)
  

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,
                             data=dataslice,
                             sphere=sphere,
                             relative_peak_threshold=.5,
                             min_separation_angle=25,
                             mask=mask,
                             return_odf=False,
                             normalize_peaks=True)
  
  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)
  

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

::
  
  gqpeaks = peaks_from_model(model=gqmodel,
                             data=dataslice,
                             sphere=sphere,
                             relative_peak_threshold=.5,
                             min_separation_angle=25,
                             mask=mask,
                             return_odf=True,
                             normalize_peaks=True)
  

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

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

True

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.

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

.. include:: ../links_names.inc



.. admonition:: Example source code

   You can download :download:`the full source code of this example <./reconst_gqi.py>`. This same script is also included in the dipy source distribution under the :file:`doc/examples/` directory.