-
Notifications
You must be signed in to change notification settings - Fork 1
/
fiber_tracking.py
144 lines (111 loc) · 7.95 KB
/
fiber_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
import os
import shutil
import subprocess
import csv
import paths
from utils import ReadImage
import numpy as np
def MoveLesionsMNIMask(isles_train_path, dst_dir):
'''Move lesion in MNI251 space from isles_train_path to dst_dir'''
lesion_mni_path = [os.path.join(root, name) for root, dirs, files in os.walk(isles_train_path) for name in files if 'MNI152_T1_1mm' in name and 'prob' not in name and name.endswith('nii.gz')]
for src in lesion_mni_path:
lesion = os.path.split(src)[1]
print('Moving ', lesion)
dst = os.path.join(dst_dir, lesion)
shutil.copy(src, dst)
# Change the working directory to dsistudio
work_dir = paths.dsi_studio_path
# create a folder called gt_stroke to store all stroke lesion in MNI152 space
dst_dir = os.path.join(work_dir, 'gt_stroke')
if not os.path.exists(dst_dir):
os.mkdir(dst_dir)
# move the strokes in MNI space to gt_stroke directory
if not os.listdir(dst_dir):
MoveLesionsMNIMask(paths.isles2017_training_dir, dst_dir)
# create a folder to save network measures
network_measures_dir = os.path.join(work_dir, 'network_measures')
if not os.path.exists(network_measures_dir):
os.mkdir(network_measures_dir)
if not os.path.exists(os.path.join(network_measures_dir, 'gt_stroke')):
os.mkdir(os.path.join(network_measures_dir, 'gt_stroke'))
# create a folder to save connectogram
connectogram_dir = os.path.join(work_dir, 'connectogram')
if not os.path.exists(connectogram_dir):
os.mkdir(connectogram_dir)
if not os.path.exists(os.path.join(connectogram_dir, 'gt_stroke')):
os.mkdir(os.path.join(connectogram_dir, 'gt_stroke'))
# create a folder to save connectivity matrix
connectivity_dir = os.path.join(work_dir, 'connectivity')
if not os.path.exists(connectivity_dir):
os.mkdir(connectivity_dir)
if not os.path.exists(os.path.join(connectivity_dir, 'gt_stroke')):
os.mkdir(os.path.join(connectivity_dir, 'gt_stroke'))
# The average fiber tracts from 1021 HCP subjects
source ='HCP1021.1mm.fib.gz'
assert (os.path.exists(os.path.join(work_dir, source))), 'HCP1021 template is not in the dsi studio directory'
# tracking parameter
#parameter_id = '--parameter_id=3D69233E9A99193F32318D24ba3Fba3Fb404b0FA4340420Fca01cb01d'
# change the working directory
os.chdir(work_dir)
# find all stroke files and sort them
stroke_dir = dst_dir
stroke_files_dir = os.listdir(dst_dir)
stroke_files_dir.sort()
assert(len(stroke_files_dir)==43)
# region property
# seed in the whole brain region and find the tracts passing through stroke lesion
region_prop ='--roi='
# seed in the stroke lesion regions and find the tracts passing through stroke lesion
#region_prop = '--seed='
#region_prop ='--roa='
# pass type of connectivity matrices
for idx, stroke_file in enumerate(stroke_files_dir):
pat_name = stroke_file[:stroke_file.find('_MNI152_T1_1mm.nii.gz')]
print(idx, 'Working on creating pass-type of connectivity matrix of ', pat_name)
#stroke_nda = ReadImage(os.path.join(stroke_dir, stroke_file))
#number_of_seed = np.count_nonzero(stroke_nda)
#number_of_seed = 964748
number_of_seed = 2235858
# connectivity matrix's parameters
connectivity_type = '--connectivity_type=pass'
connectivity_value = '--connectivity_value=count'
connectivity_threshold = '--connectivity_threshold=0'
#subprocess.call(['./dsi_studio', '--action=trk', '--source='+source, region_prop+os.path.join(stroke_dir, stroke_file), parameter_id,'--seed_count='+str(number_of_seed), '--output=no_file', '--connectivity=aal', connectivity_type, connectivity_value, connectivity_threshold])
# fiber tracking and setting the tracking parameters
subprocess.call(['./dsi_studio', '--action=trk', '--source='+source, region_prop+os.path.join(stroke_dir, stroke_file), '--seed_count='+str(number_of_seed), '--fa_threshold=0.15958', '--seed_plan=1', '--initial_dir=2', '--interpolation=0', '--turning_angle=90.0', '--step_size=.5', '--smoothing=.5', '--min_length=3', '--max_length=500', '--tip_iteration=1', '--thread_count=8', '--output=no_file', '--connectivity=aal', connectivity_type, connectivity_value, connectivity_threshold])
# move the files
network_measure_files = [os.path.join(root, name) for root, dirs, files in os.walk(work_dir) for name in files if 'network_measures' in name and name.endswith('.txt')]
network_measure_file_dst = os.path.join(os.path.split(network_measure_files[0])[0], 'network_measures', 'gt_stroke', os.path.split(network_measure_files[0].replace(source, pat_name))[1])
shutil.move(network_measure_files[0], network_measure_file_dst)
connectogram_files = [os.path.join(root, name) for root, dirs, files in os.walk(work_dir) for name in files if 'connectogram' in name and name.endswith('.txt')]
connectogram_file_dst = os.path.join(os.path.split(connectogram_files[0])[0], 'connectogram', 'gt_stroke', os.path.split(connectogram_files[0].replace(source, pat_name))[1])
shutil.move(connectogram_files[0], connectogram_file_dst)
connectivity_files = [os.path.join(root, name) for root, dirs, files in os.walk(work_dir) for name in files if 'connectivity' in name and name.endswith('.mat')]
connectivity_file_dst = os.path.join(os.path.split(connectivity_files[0])[0], 'connectivity', 'gt_stroke', os.path.split(connectivity_files[0].replace(source, pat_name))[1])
shutil.move(connectivity_files[0], connectivity_file_dst)
print('---'*30)
# end type of connectivity matrices
for idx, stroke_file in enumerate(stroke_files_dir):
pat_name = stroke_file[:stroke_file.find('_MNI152_T1_1mm.nii.gz')]
print(idx, 'Working on creating end-type of connectivity matrix of ', pat_name)
#stroke_nda = ReadImage(os.path.join(stroke_dir, stroke_file))
#number_of_seed = np.count_nonzero(stroke_nda)
#number_of_seed = 964748
number_of_seed = 2235858
# connectivity matrix's parameters
connectivity_type = '--connectivity_type=end'
connectivity_value = '--connectivity_value=count'
connectivity_threshold = '--connectivity_threshold=0'
# fiber tracking and setting the tracking parameters
#subprocess.call(['./dsi_studio', '--action=trk', '--source='+source, region_prop+os.path.join(stroke_dir, stroke_file), parameter_id, '--output=no_file', '--connectivity=aal', connectivity_type, connectivity_value, connectivity_threshold])
subprocess.call(['./dsi_studio', '--action=trk', '--source='+source, region_prop+os.path.join(stroke_dir, stroke_file), '--seed_count='+str(number_of_seed), '--fa_threshold=0.15958', '--seed_plan=1', '--initial_dir=2', '--interpolation=0', '--turning_angle=90.0', '--step_size=.5', '--smoothing=.5', '--min_length=3', '--max_length=500', '--tip_iteration=1' '--thread_count=8', '--output=no_file', '--connectivity=aal', connectivity_type, connectivity_value, connectivity_threshold])
network_measure_files = [os.path.join(root, name) for root, dirs, files in os.walk(work_dir) for name in files if 'network_measures' in name and name.endswith('.txt')]
network_measure_file_dst = os.path.join(os.path.split(network_measure_files[0])[0], 'network_measures', 'gt_stroke', os.path.split(network_measure_files[0].replace(source, pat_name))[1])
shutil.move(network_measure_files[0], network_measure_file_dst)
connectogram_files = [os.path.join(root, name) for root, dirs, files in os.walk(work_dir) for name in files if 'connectogram' in name and name.endswith('.txt')]
connectogram_file_dst = os.path.join(os.path.split(connectogram_files[0])[0], 'connectogram', 'gt_stroke', os.path.split(connectogram_files[0].replace(source, pat_name))[1])
shutil.move(connectogram_files[0], connectogram_file_dst)
connectivity_files = [os.path.join(root, name) for root, dirs, files in os.walk(work_dir) for name in files if 'connectivity' in name and name.endswith('.mat')]
connectivity_file_dst = os.path.join(os.path.split(connectivity_files[0])[0], 'connectivity', 'gt_stroke', os.path.split(connectivity_files[0].replace(source, pat_name))[1])
shutil.move(connectivity_files[0], connectivity_file_dst)
print('---'*30)