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: Download and get the data filenames for this tutorial. img contains a nibabel Nifti1Image object (data) and gtab contains a
Load the raw diffusion data and the affine. data.shape This dataset has anisotropic voxel sizes, therefore reslicing is necessary. Instantiate the model and apply it to the data. The parameter Lets just use one slice only from the data. Load an ODF reconstruction sphere Calculate the ODFs with this specific sphere ODF.shape Using It is also possible to calculate GFA. With parameter This ODF will be of course identical to the ODF calculated above as long as
the same data and mask are used. True The advantage of using Yeh, F-C et al., Generalized Q-sampling imaging, IEEE Transactions
on Medical Imaging, vol 29, no 9, 2010. Example source code You can download Reconstruct with Generalized Q-Sampling Imaging
import numpy as np
from dipy.core.gradients import gradient_table
from dipy.data import get_fnames, get_sphere
from dipy.io.gradients import read_bvals_bvecs
from dipy.io.image import load_nifti
from dipy.reconst.gqi import GeneralizedQSamplingModel
from dipy.direction import peaks_from_model
fraw, fbval, fbvec = get_fnames('taiwan_ntu_dsi')
GradientTable
object (gradient information e.g. b-values). For example to
read the b-values it is possible to write:print(gtab.bvals)
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)
(96, 96, 60, 203)
gqmodel = GeneralizedQSamplingModel(gtab, sampling_length=3)
sampling_length
is used here todataslice = data[:, :, data.shape[2] // 2]
mask = dataslice[..., 0] > 50
gqfit = gqmodel.fit(dataslice, mask=mask)
sphere = get_sphere('repulsion724')
ODF = gqfit.odf(sphere)
print('ODF.shape (%d, %d, %d)' % ODF.shape)
(96, 96, 724)
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
GFA = gqpeaks.gfa
print('GFA.shape (%d, %d)' % GFA.shape)
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)
print(np.sum(gqpeaks.odf != ODF) == 0)
peaks_from_model
is that it calculates the ODF only
once and saves it or deletes if it is not necessary to keep.the full source code of this example
. This same script is also included in the dipy source distribution under the doc/examples/
directory.