-
Notifications
You must be signed in to change notification settings - Fork 3
/
init_tracking.py
226 lines (184 loc) · 7.56 KB
/
init_tracking.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
# The code below is an example of how to track using Dipy on data on Karst
#
# The script assumes that your environment is set up correctly.
#
# This script will require the following set of commands:
#
# (1) module unload python
# (2) module load anaconda2
# (3) git clone git@github.com:nipy/dipy.git
# (4) git clone git@github.com:nipy/nibabel.git
# (5) Make sure that Dipy is on the PYTHONPATH:
# (5.1) add Dipy to the shell environment
# Pythonpath is an environmental
# variable that keeps tracks of
# what Python can access to.
# export PYTHONPATH =
# "/N/dc2/projects/lifebid/code/franpest/dipy/:
# /N/dc2/projects/lifebid/code/franpest/nibabel:
# $PYTHONPATH"
# (5.2) cd dipy/
# make clean # makes a
# make ext # compiles all the files in dipy that
# - needs compilations, compiles them and
# - saves the compiled copies loacally under
# - dipy (not on the system level).
# - This is the command being called:
# - python setup.py built_ext -inplace)
#
# (5.3) open IPython; in the shell type ipython
#
#
# Franco Pestilli
# Import tools to environment
import numpy as np
import nibabel as nib
import dipy as dipy
# Load methods from Dipy for tracking and signal voxel reconstruction
from dipy.reconst.dti import TensorModel, fractional_anisotropy
from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,
auto_response, response_from_mask)
from dipy.direction import peaks_from_model
from dipy.tracking.eudx import EuDX
from dipy.data import get_sphere
from dipy.io.gradients import read_bvals_bvecs
from dipy.core.gradients import gradient_table
from dipy.align.reslice import reslice
# A few requirements for visualizing the outputs (VTK)
from dipy.viz import fvtk
from dipy.viz.colormap import line_colors
from dipy.viz import actor, window
# We initialize pointers to file paths
data_path = '/N/dc2/projects/lifebid/HCP7/108323/'
data_file = data_path + 'diffusion_data/'+'data.nii.gz'
data_bvec = data_path + 'diffusion_data/'+'data.bvec'
data_bval = data_path + 'diffusion_data/'+'data.bval'
data_brainmask = data_path + 'diffusion_data/'+'nodif_brain_mask.nii.gz'
fs_path = '/N/dc2/projects/lifebid/HCP/Brent/7t_rerun/freesurfer/7t_108323'
data_fs_seg = fs_path + '/mri/aparc+aseg+LAS.nii.gz'
# Load the data
dmri_image = nib.load(data_file)
dmri = dmri_image.get_data()
affine = dmri_image.affine
brainmask_im = nib.load(data_brainmask)
brainmask = brainmask_im.get_data()
bm_affine = brainmask_im.affine
aparc_im = nib.load(data_fs_seg)
aparc = aparc_im.get_data()
aparc_affine = brainmask_im.affine
# Freesurfer parecellation ROIs
wm_regions = [2, 41, 16, 17, 28, 60, 51, 53, 12, 52, 12, 52, 13, 18,
54, 50, 11, 251, 252, 253, 254, 255, 10, 49, 46, 7]
wm_mask = np.zeros(aparc.shape)
for l in wm_regions:
wm_mask[aparc==l] = 1
callosal_regions = [255]
callosum = np.zeros(aparc.shape)
for c in callosal_regions:
callosum[aparc==c] = 1
# reslice the masks to dmri space
#current_zooms = aparc_im.header.get_zooms()[:3]
#new_zooms = dmri_image.header.get_zooms()[:3]
#callosum_r, call_r_affine = reslice(callosum, aparc_affine,
# current_zooms,
# new_zooms)
#wm_mask_r, wm_mask_r_affine = reslice(wm_mask, aparc_affine,
# current_zooms,
# new_zooms)
show_two_slices(volume1=dmri[...,2],
affine1=affine,
volume2=wm_mask,
affine2=aparc_affine,
show_axes=True,
shift=200, opacity=[0.8, 0.7])
# Load the bvals and bvecs from disk
# ideally we should use the following:
# bvals, bvecs = read_bvals_bvecs(data_bval, data_bvec)
# in practice this data set does not conform to a standard (FSL)
bvals = np.loadtxt(data_bval,delimiter=',')
bvecs = np.loadtxt(data_bvec).T
gtab = gradient_table(bvals, bvecs, b0_threshold=100)
# After loading all files we can start voxel reconstruction
# Estimate impulse response for deconvolution (check that the estimation is good)
# http://nipy.org/dipy/examples_built/reconst_csd.html
# response, ratio = response_from_mask(gtab, dmri, callosum)
res, ratio = response_from_mask (gtab, dmri, callosum_r)
# Build a CSD model kernek
csd_model = ConstrainedSphericalDeconvModel(gtab, res);
# Build the Diffusion Tensor Model.
dt_model = TensorModel(gtab)
# Initialize a sphere
sphere = get_sphere('repulsion724')
# Find fiber-peaks from the CSD fit
csd_peaks = peaks_from_model(model=csd_model,
data=dmri,
sphere=sphere,
mask=brainmask,
relative_peak_threshold=.5,
min_separation_angle=25,
parallel=True,
nbr_processes=10)
print('done estimating CSD_peaks')
# Fit the diffusion tensor model to the data
dt_fit = dt_model.fit(data=dmri, mask=brainmask)
print('done estimating DTI')
# Next we visualize the peaks in a single brain slice
ren = fvtk.ren()
slice = 70
fodf_peaks = fvtk.peaks(csd_peaks.peak_dirs[:,:,slice], csd_peaks.peak_values[:,:,slice], scale=1.3)
fvtk.add(ren, fodf_peaks)
fvtk.show(ren)
ren.azimuth(90)
fvtk.show(ren)
# We now classify the tissue to decide where to stop tracking.
from dipy.tracking.local import ThresholdTissueClassifier
classifier = ThresholdTissueClassifier(dt_fit.fa, .1)
# We will need some seeds for tracking. There are many different functions for seeds.
# Hereafter we will use seeds from a precomputed mask, for us that mask will be the WM mask
from dipy.tracking import utils
seeds = utils.seeds_from_mask(wm_mask_r, density=1)
seeds_im = fvtk.dots(seeds, (1, 0.5, 0))
#fvtk.add(ren,seeds_im)
# Now we have prepared all we need and we can track
from dipy.tracking.local import LocalTracking
# Creating the tracking model
streamlines = LocalTracking(csd_peaks, classifier, seeds, affine=np.eye(4), step_size=.5)
# Tracking
streamlines = list(streamlines)
# visualzie the streamlines
streamlines_actor = fvtk.line(streamlines[:100000])
#fvtk.rm(ren,fodf_peaks)
#fvtk.rm(ren, seeds_im)
fvtk.add(ren,streamlines_actor)
fvtk.show(ren)
from nibabel.streamlines import Tractogram, save
tractogram = Tractogram(streamlines, affine_to_rasmm=affine)
save(tractogram, 'test.trk')
# For the next time we will save to disk.
# After that we will do Anatomically Constrained Tracking.
def show_slice(volume, affine=None, show_axes=False, k=None):
ren = window.Renderer()
slicer_actor = actor.slicer(volume, affine)
slicer_actor.display(None, None, k)
ren.add(slicer_actor)
if show_axes:
ren.add(actor.axes((100, 100, 100)))
window.show(ren)
def show_two_slices(volume1, affine1, volume2, affine2=None,
show_axes=False, k=None, shift=None, opacity=[0.8, 0.4]):
ren = window.Renderer()
slicer_actor = actor.slicer(volume1, affine1)
slicer_actor.display(None, None, k)
ren.add(slicer_actor)
slicer_actor2 = actor.slicer(volume2, affine2)
if shift is not None:
slicer_actor2.SetPosition(shift, 0, 0)
slicer_actor2.display(None, None, k)
ren.add(slicer_actor2)
if opacity is not None:
slicer_actor.opacity(opacity[0])
slicer_actor2.opacity(opacity[1])
if show_axes:
ren.add(actor.axes((100, 100, 100)))
window.show(ren)
# END