Small script that fetch DSI data and uses Dipy recon module to recreate figure 6 and 7 from Paquette et al. "Optimal DSI reconstruction parameter recommendations: better ODFs and better connectivity".

dsi_parameters_paper_fig_6_7.py

#!/usr/bin/env python
"""
This example reproduces the ODFs in figure 6 and 7 from Paquette et al. 
"Optimal DSI reconstruction parameter recommendations: better ODFs and better connectivity".

The Optimal DSI (figure 7), is obtained using the parameters from the last
line of Tbl. 5 as the snr was estimated at 38 following an approach detailed 
in https://github.com/nipy/dipy/blob/master/doc/examples/snr_in_cc.py

Required python modules (with their respective reqs): 
numpy: http://www.numpy.org/
nibabel: http://nipy.org/nibabel/
dipy: http://nipy.org/dipy/
"""

import os
from os.path import join as pjoin

from dipy.data.fetcher import check_md5, _get_file_data

import numpy as np
import nibabel as nib

from dipy.core.gradients import gradient_table
from dipy.io.gradients import read_bvals_bvecs

from dipy.reconst.dti import TensorModel
from dipy.reconst.dsi import DiffusionSpectrumModel, DiffusionSpectrumDeconvModel
from dipy.data import get_sphere

from dipy.viz import fvtk


# Matplotlib version check
# The visualisation is tested for 1.4 and known to fail for 1.2 and 1.3
import matplotlib as plt

plt_ver = int(plt.__version__.split('.')[0])
plt_subver = int(plt.__version__.split('.')[1])

if (plt_ver < 1) or ((plt_ver == 1) and (plt_subver < 4)):
    print("The visualisation part of this example requires matplotlib >= 1.4 (currently using {})".format(plt.__version__))


def main():

    # Get the data and gradients
    fetch_dsi515_crop()
    img, bvals, bvecs = read_dsi515_crop()
    gtab = gradient_table(bvals, bvecs)
    data = img.get_data()
    aff  = img.get_affine()

    # Estimate tensor
    tenmodel = TensorModel(gtab)
    tensorfit = tenmodel.fit(data)

    print("Computing tensor")
    fa = tensorfit.fa.astype(np.float32)
    img_fa = nib.nifti1.Nifti1Image(fa, aff)
    print('Saving FA as fa_dsi515_crop.nii.gz')
    nib.save(img, 'fa_dsi515_crop.nii.gz')
    del img_fa

    # Approximate white matter mask by FA threshold
    wm_mask = fa > 0.3

    # Get sphere for ODF computation
    print("Modify tess_order (line 59) to 1 or 2 for higher ODF tessellation")
    tess_order = 0
    sphere = get_sphere('symmetric724').subdivide(tess_order)

    ## Parameters of DSI:
    # G: Grid size
    # a: lower ODF integration relative bound
    # b: upper ODF integration relative bound
    # B: beta parameter for low-pass Kaiser window


    ## The DSI models:

    # Plain DSI (PDSI)

    [G, a, b, B] = [11, 0.0, 1.00, 0.0]
    width = beta2width(B)
    half_grid_size = (G+1)//2
    pdsimodel = DiffusionSpectrumModel(gtab,
                            qgrid_size = G,
                            r_start = a * half_grid_size,
                            r_end = b * half_grid_size,
                            r_step = (b - a) * half_grid_size / 19.,
                            filter_width = width,
                            normalize_peaks = False)

    # Classic DSI (CDSI)

    [G, a, b, B] = [17, 0.25, 0.75, 2.5]
    width = beta2width(B)
    half_grid_size = (G+1)//2
    cdsimodel = DiffusionSpectrumModel(gtab,
                            qgrid_size = G,
                            r_start = a * half_grid_size,
                            r_end = b * half_grid_size,
                            r_step = (b - a) * half_grid_size / 19.,
                            filter_width = width,
                            normalize_peaks = False)

    # Optimal DSI (ODSI) estimated SNR 38

    [G, a, b, B] = [35, 0.4, 0.8, 0.0]
    width = beta2width(B)
    half_grid_size = (G+1)//2
    odsimodel = DiffusionSpectrumModel(gtab,
                            qgrid_size = G,
                            r_start = a * half_grid_size,
                            r_end = b * half_grid_size,
                            r_step = (b - a) * half_grid_size / 19.,
                            filter_width = width,
                            normalize_peaks = False)

    # DSI Deconvolution (DDSI)

    [G, a, b, B] = [35, 0.2, 0.8, 0.0]
    width = beta2width(B)
    half_grid_size = (G+1)//2
    ddsimodel = DiffusionSpectrumDeconvModel(gtab,
                            qgrid_size = G,
                            r_start = a * half_grid_size,
                            r_end = b * half_grid_size,
                            r_step = (b - a) * half_grid_size / 19.,
                            filter_width = width,
                            normalize_peaks = False)



    ## Computing and saving ODFs

    # ODFs from figure 6: PDSI
    print("Computing ODFs from figure 6: PDSI")
    dsifit = pdsimodel.fit(data = data, mask = wm_mask)
    dsiodf = dsifit.odf(sphere)
    print("Saving ODFs as odfs_pdsi.png")
    save_ss(dsiodf, sphere, fname = 'odfs_pdsi.png')

    # ODFs from figure 6: CDSI
    print("Computing ODFs from figure 6: CDSI")
    dsifit = cdsimodel.fit(data = data, mask = wm_mask)
    dsiodf = dsifit.odf(sphere)
    print("Saving ODFs as odfs_cdsi.png")
    save_ss(dsiodf, sphere, fname = 'odfs_cdsi.png')

    # ODFs from figure 7: ODSI
    print("Computing ODFs from figure 7: ODSI")
    dsifit = odsimodel.fit(data = data, mask = wm_mask)
    dsiodf = dsifit.odf(sphere)
    print("Saving ODFs as odfs_odsi.png")
    save_ss(dsiodf, sphere, fname = 'odfs_odsi.png')

    # ODFs from figure 7: DDSI
    print("Computing ODFs from figure 7: DDSI")
    dsifit = ddsimodel.fit(data = data, mask = wm_mask)
    dsiodf = dsifit.odf(sphere)
    print("Saving ODFs as odfs_ddsi.png")
    save_ss(dsiodf, sphere, fname = 'odfs_ddsi.png')


def fetch_dsi515_crop():
    """ Download the DSI515 crop dataset
    """
    dipy_home = pjoin(os.path.expanduser('~'), '.dipy')
    uraw  = 'https://www.dropbox.com/s/uot8siaw01vaa7n/crop.nii.gz?dl=1'
    ubval = 'https://www.dropbox.com/s/jlnko3373vqycwn/dsi-scheme.bval?dl=1'
    ubvec = 'https://www.dropbox.com/s/ogr4dp3w4eorb9d/dsi-scheme.bvec?dl=1'
    folder = pjoin(dipy_home, 'dsi515_crop')


    md5_list = ['be22e7c55209e66d77088d561a5391fd', # data
                '1b6053f8ffdbdc81bf7db3f10f235c3d', # bval
                '622a3e0ed53d6832cb77d1a2d9d6f426'] # bvec

    url_list = [uraw, ubval, ubvec]
    fname_list = ['dsi515_crop.nii.gz', 'dsi515.bval', 'dsi515.bvec']

    if not os.path.exists(folder):
        print('Creating new directory %s' % folder)
        os.makedirs(folder)
        print('Downloading raw DSI crop data (<1MB)...')

        for i in range(len(md5_list)):
            _get_file_data(pjoin(folder, fname_list[i]), url_list[i])
            check_md5(pjoin(folder, fname_list[i]), md5_list[i])

        print('Done.')
        print('Files copied in folder %s' % folder)
    else:
        print('Dataset is already in place. If you want to fetch it again, please first remove the folder %s ' % folder)


def read_dsi515_crop():
    """ Load DSI515 crop dataset

    Returns
    -------
    img : obj,
        Nifti1Image
    gtab : obj,
        GradientTable
    """
    dipy_home = pjoin(os.path.expanduser('~'), '.dipy')
    folder = pjoin(dipy_home, 'dsi515_crop')
    fraw = pjoin(folder, 'dsi515_crop.nii.gz')
    fbval = pjoin(folder, 'dsi515.bval')
    fbvec = pjoin(folder, 'dsi515.bvec')

    md5_dict = {'data': 'be22e7c55209e66d77088d561a5391fd',
                'bval': '1b6053f8ffdbdc81bf7db3f10f235c3d',
                'bvec': '622a3e0ed53d6832cb77d1a2d9d6f426'}

    check_md5(fraw, md5_dict['data'])
    check_md5(fbval, md5_dict['bval'])
    check_md5(fbvec, md5_dict['bvec'])

    bvals, bvecs = read_bvals_bvecs(fbval, fbvec)

    # gtab = gradient_table(bvals, bvecs)
    img = nib.load(fraw)
    return img, bvals, bvecs#, gtab

def beta2width(beta):
    """
    Convert the Kaiser window Beta parameter to its 
    Hann window equivalent following Appendix A
    """
    if beta != 0:
        return 2 * np.pi * 11 / np.arccos((2 / np.i0(beta)) - 1)
    else:
        # no filter
        return np.inf

def save_ss(odf, sphere, fname = None):
    """
    Visualise or save as PNG a screenshot of a slice of ODFs
    """
    r = fvtk.ren()
    sfu = fvtk.sphere_funcs(odf, sphere, scale=2.2, norm=True)
    sfu.RotateX(90)
    sfu.RotateY(180)
    fvtk.add(r, sfu)
    if fname == None:
        fvtk.show(r)
    else:
        fvtk.record(r, n_frames=1, out_path= fname, size=(3000, 1500), magnification = 2)


if __name__ == "__main__":
    main()