In [1]:
import numpy as np
import nibabel as nib
import dsiadapt as dsi
from dipy.viz import fvtk
import dipy.core.gradients as grad
from dipy.data import get_sphere
from dipy.tracking.local import LocalTracking
from dipy.viz import fvtk
from dipy.viz.colormap import line_colors
from dipy.tracking import utils
from dipy.tracking.local import ThresholdTissueClassifier
from dipy.direction import peaks_from_model
import time
from dipy.io.trackvis import save_trk
from dipy.io.pickles import save_pickle
from dipy.io.pickles import load_pickle
import glob
from dipy.io.trackvis import save_trk
%pylab inline
np.set_printoptions(threshold=numpy.nan)
In [2]:
# ODF sphere directions
sphere = get_sphere('symmetric724')
#sphere = sphere.subdivide(2)
In [3]:
datapath = '/data/qytian/QTResearchResult/DSIQspace/DSIData'; # Whole brain data is only locally available
wintypeArr = np.array(['none', 'hanning', 'hamming', 'blackman', 'none']);
qgridsz = 17;
for ii in np.arange(fnArr.shape[0]):
fn = fnArr[ii];
data = nib.load(datapath + '/' + fn + '/' + fn + '.nii.gz').get_data();
affine = np.identity(4);
wmmask = nib.load(datapath + '/' + fn + '/' + fn + '_wmmask.nii.gz').get_data();
brainmask = nib.load(datapath + '/' + fn + '/' + fn + '_masked_mask.nii.gz').get_data();
classifier = ThresholdTissueClassifier(brainmask, 0.5)
gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs_dipy.txt');
seedmask = nib.load('data/' + fn + '_roi_mask.nii.gz').get_data();
mdd = np.loadtxt('data/' + fn + '_stats.txt')[0];
fov = np.loadtxt('data/' + fn + '_stats.txt')[1];
intdistArr = np.array([fov/2, fov/2, fov/2, fov/2, mdd]);
for jj in np.arange(wintypeArr.shape[0]):
wintype = wintypeArr[jj];
intdist = intdistArr[jj];
rend = intdist / fov * qgridsz;
winwidth = 2. * dsi.create_qtable(gtab).max();
rend = intdist / fov * qgridsz;
dsimodel = dsi.DiffusionSpectrumModel(gtab, qgrid_size=qgridsz,
filter_width=winwidth, filter_type=wintype, r_start=0, r_end=rend, r_step=0.2);
peakname = 'figure6/' + fn + '_peak_qgridsz' + str(qgridsz) + '_odfdirnum' + str(sphere.vertices.shape[0]) + '_' + wintype + '_rend' + str(intdist) + '.pkl';
if len(glob.glob(peakname)) == 0:
# Compute ODF peaks
tic = time.time()
dsipeaks = peaks_from_model(dsimodel, data, sphere,
relative_peak_threshold=0.2,
min_separation_angle=25, normalize_peaks=True, mask=brainmask);
toc = time.time() - tic
disp(toc)
save_pickle(peakname, dsipeaks);
else:
dsipeaks = load_pickle(peakname);
# Track
dennum = 2;
seeds = utils.seeds_from_mask(seedmask, density=[dennum, dennum, dennum], affine=affine)
# Initialization of LocalTracking. The computation happens in the next step.
streamlines = LocalTracking(dsipeaks, classifier, seeds, affine, step_size=.5)
# Compute streamlines and store as a list.
streamlines = list(streamlines)
# Prepare the display objects.
color = line_colors(streamlines)
streamlines_actor = fvtk.line(streamlines, line_colors(streamlines))
trackname = 'figure6/' + fn + '_track_qgridsz' + str(qgridsz) + '_odfdirnum' + str(sphere.vertices.shape[0]) + '_' + wintype + '_rend' + str(intdist) + '_density' + str(dennum) + '.trk';
save_trk(trackname, streamlines, affine, wmmask.shape)
# Create the 3d display.
# r = fvtk.ren()
# fvtk.add(r, streamlines_actor)
# fvtk.camera(r, pos=(0, 1, 0), viewup=(0, 0, 1))
# fvtk.show(r)