In [1]:

    

import numpy as np
import nibabel as nib
import dipy.reconst.dti as dti
from dipy.data import fetch_stanford_hardi
fetch_stanford_hardi()
from dipy.data import read_stanford_hardi


    

    




Dataset is already in place. If you want to fetch it again please first remove the folder C:\Users\lk_zh\.dipy\stanford_hardi 

In [2]:

    

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


    

    




Dataset is already in place. If you want to fetch it again please first remove the folder C:\Users\lk_zh\.dipy\stanford_hardi 
data.shape (81, 106, 76, 160)

In [3]:

    

from dipy.segment.mask import median_otsu
maskdata, mask = median_otsu(data, 3, 1, True,
                             vol_idx=range(10, 50), dilate=2)
print('maskdata.shape (%d, %d, %d, %d)' % maskdata.shape)


    

    




C:\Users\lk_zh\Anaconda2\lib\site-packages\skimage\filter\__init__.py:6: skimage_deprecation: The `skimage.filter` module has been renamed to `skimage.filters`.  This placeholder module will be removed in v0.13.
  warn(skimage_deprecation('The `skimage.filter` module has been renamed '






    




maskdata.shape (71, 87, 62, 160)

In [4]:

    

tenmodel = dti.TensorModel(gtab)
tenfit = tenmodel.fit(maskdata)


    

In [5]:

    

print('Computing anisotropy measures (FA, MD, RGB)')
from dipy.reconst.dti import fractional_anisotropy, color_fa, lower_triangular

FA = fractional_anisotropy(tenfit.evals)


    

    




Computing anisotropy measures (FA, MD, RGB)

In [6]:

    

FA[np.isnan(FA)] = 0


    

In [7]:

    

FA.shape


    

    
Out[7]:





(71L, 87L, 62L)

In [8]:

    

fa_img = nib.Nifti1Image(FA.astype(np.float32), img.affine)
nib.save(fa_img, 'tensor_fa.nii.gz')


    

In [9]:

    

fa_img.get_data().shape


    

    
Out[9]:





(71L, 87L, 62L)

In [10]:

    

evecs_img = nib.Nifti1Image(tenfit.evecs.astype(np.float32), img.affine)
nib.save(evecs_img, 'tensor_evecs.nii.gz')


    

In [11]:

    

evecs_img.get_data().shape


    

    
Out[11]:





(71L, 87L, 62L, 3L, 3L)

In [12]:

    

MD1 = dti.mean_diffusivity(tenfit.evals)
nib.save(nib.Nifti1Image(MD1.astype(np.float32), img.affine), 'tensors_md.nii.gz')


    

In [13]:

    

MD1.shape


    

    
Out[13]:





(71L, 87L, 62L)

In [14]:

    

MD2 = tenfit.md


    

In [15]:

    

MD2.shape


    

    
Out[15]:





(71L, 87L, 62L)

In [16]:

    

FA = np.clip(FA, 0, 1)
RGB = color_fa(FA, tenfit.evecs)
nib.save(nib.Nifti1Image(np.array(255 * RGB, 'uint8'), img.affine), 'tensor_rgb.nii.gz')


    

In [17]:

    

print FA.shape, RGB.shape


    

    




(71L, 87L, 62L) (71L, 87L, 62L, 3L)

In [18]:

    

print('Computing tensor ellipsoids in a part of the splenium of the CC')

from dipy.data import get_sphere
sphere = get_sphere('symmetric724')

from dipy.viz import fvtk
ren = fvtk.ren()

evals = tenfit.evals[13:43, 44:74, 28:29]
evecs = tenfit.evecs[13:43, 44:74, 28:29]


    

    




Computing tensor ellipsoids in a part of the splenium of the CC

In [19]:

    

print evals.shape, evecs.shape


    

    




(30L, 30L, 1L, 3L) (30L, 30L, 1L, 3L, 3L)

In [20]:

    

cfa = RGB[13:43, 44:74, 28:29]
cfa /= cfa.max()

fvtk.add(ren, fvtk.tensor(evals, evecs, cfa, sphere))

print('Saving illustration as tensor_ellipsoids.png')
fvtk.record(ren, n_frames=1, out_path='tensor_ellipsoids.png', size=(600, 600))


    

    




Saving illustration as tensor_ellipsoids.png

In [21]:

    

print cfa.shape


    

    




(30L, 30L, 1L, 3L)

In [22]:

    

fvtk.clear(ren)

tensor_odfs = tenmodel.fit(data[20:50, 55:85, 38:39]).odf(sphere)

fvtk.add(ren, fvtk.sphere_funcs(tensor_odfs, sphere, colormap=None))
fvtk.show(ren)
print('Saving illustration as tensor_odfs.png')
fvtk.record(ren, n_frames=1, out_path='tensor_odfs.png', size=(600, 600))


    

    




Saving illustration as tensor_odfs.png

In [23]:

    

fa_img = nib.load('tensor_fa.nii.gz')
FA = fa_img.get_data()
evecs_img = nib.load('tensor_evecs.nii.gz')
evecs = evecs_img.get_data()


    

In [24]:

    

print FA.shape, evecs.shape


    

    




(71L, 87L, 62L) (71L, 87L, 62L, 3L, 3L)

In [25]:

    

print FA.astype('f8').shape


    

    




(71L, 87L, 62L)

In [26]:

    

FA[np.isnan(FA)] = 0
from dipy.data import get_sphere

sphere = get_sphere('symmetric724')
from dipy.reconst.dti import quantize_evecs

peak_indices = quantize_evecs(evecs, sphere.vertices)

from dipy.tracking.eudx import EuDX

eu = EuDX(FA.astype('f8'), peak_indices, seeds=50000, odf_vertices = sphere.vertices, a_low=0.2)

tensor_streamlines = [streamline for streamline in eu]


    

In [27]:

    

hdr = nib.trackvis.empty_header()
hdr['voxel_size'] = fa_img.get_header().get_zooms()[:3]
hdr['voxel_order'] = 'LAS'
hdr['dim'] = FA.shape


    

    




C:\Users\lk_zh\Anaconda2\lib\site-packages\ipykernel\__main__.py:2: DeprecationWarning: get_header method is deprecated.
Please use the ``img.header`` property instead.

* deprecated from version: 2.1
* Will raise <class 'nibabel.deprecator.ExpiredDeprecationError'> as of version: 4.0
  from ipykernel import kernelapp as app

In [28]:

    

tensor_streamlines_trk = ((sl, None, None) for sl in tensor_streamlines)

ten_sl_fname = 'tensor_streamlines.trk'


    

In [29]:

    

nib.trackvis.write(ten_sl_fname, tensor_streamlines_trk, hdr, points_space='voxel')


    

In [30]:

    

ren = fvtk.ren()
from dipy.viz.colormap import line_colors
fvtk.add(ren, fvtk.streamtube(tensor_streamlines, line_colors(tensor_streamlines)))

print('Saving illustration as tensor_tracks.png')

ren.SetBackground(1, 1, 1)
fvtk.record(ren, n_frames=1, out_path='tensor_tracks.png', size=(600, 600))


    

    




Saving illustration as tensor_tracks.png

In [31]:

    

fvtk.show(ren)


    

In [ ]: