-
Notifications
You must be signed in to change notification settings - Fork 520
/
odf.py
399 lines (363 loc) · 13.5 KB
/
odf.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
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
# -*- coding: utf-8 -*-
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""Provides interfaces to various commands provided by diffusion toolkit
"""
import os
import re
from ...utils.filemanip import fname_presuffix, split_filename, copyfile
from ..base import (
TraitedSpec,
File,
traits,
CommandLine,
CommandLineInputSpec,
isdefined,
)
__docformat__ = "restructuredtext"
class HARDIMatInputSpec(CommandLineInputSpec):
bvecs = File(
exists=True, desc="b vectors file", argstr="%s", position=1, mandatory=True
)
bvals = File(exists=True, desc="b values file", mandatory=True)
out_file = File(
"recon_mat.dat",
desc="output matrix file",
argstr="%s",
usedefault=True,
position=2,
)
order = traits.Int(
argstr="-order %s",
desc="""maximum order of spherical harmonics. must be even number. default is 4""",
)
odf_file = File(
exists=True,
argstr="-odf %s",
desc="""\
Filename that contains the reconstruction points on a HEMI-sphere.
Use the pre-set 181 points by default""",
)
reference_file = File(
exists=True,
argstr="-ref %s",
desc="""\
Provide a dicom or nifti image as the reference for the program to
figure out the image orientation information. if no such info was
found in the given image header, the next 5 options -info, etc.,
will be used if provided. if image orientation info can be found
in the given reference, all other 5 image orientation options will
be IGNORED""",
)
image_info = File(
exists=True,
argstr="-info %s",
desc="""\
specify image information file. the image info file is generated
from original dicom image by diff_unpack program and contains image
orientation and other information needed for reconstruction and
tracking. by default will look into the image folder for .info file""",
)
image_orientation_vectors = traits.List(
traits.Float(),
minlen=6,
maxlen=6,
desc="""\
specify image orientation vectors. if just one argument given,
will treat it as filename and read the orientation vectors from
the file. if 6 arguments are given, will treat them as 6 float
numbers and construct the 1st and 2nd vector and calculate the 3rd
one automatically.
this information will be used to determine image orientation,
as well as to adjust gradient vectors with oblique angle when""",
argstr="-iop %f",
)
oblique_correction = traits.Bool(
desc="""\
when oblique angle(s) applied, some SIEMENS dti protocols do not
adjust gradient accordingly, thus it requires adjustment for correct
diffusion tensor calculation""",
argstr="-oc",
)
class HARDIMatOutputSpec(TraitedSpec):
out_file = File(exists=True, desc="output matrix file")
class HARDIMat(CommandLine):
"""Use hardi_mat to calculate a reconstruction matrix from a gradient table
"""
input_spec = HARDIMatInputSpec
output_spec = HARDIMatOutputSpec
_cmd = "hardi_mat"
def _create_gradient_matrix(self, bvecs_file, bvals_file):
_gradient_matrix_file = "gradient_matrix.txt"
bvals = [val for val in re.split(r"\s+", open(bvals_file).readline().strip())]
bvecs_f = open(bvecs_file)
bvecs_x = [val for val in re.split(r"\s+", bvecs_f.readline().strip())]
bvecs_y = [val for val in re.split(r"\s+", bvecs_f.readline().strip())]
bvecs_z = [val for val in re.split(r"\s+", bvecs_f.readline().strip())]
bvecs_f.close()
gradient_matrix_f = open(_gradient_matrix_file, "w")
for i in range(len(bvals)):
if int(bvals[i]) == 0:
continue
gradient_matrix_f.write("%s %s %s\n" % (bvecs_x[i], bvecs_y[i], bvecs_z[i]))
gradient_matrix_f.close()
return _gradient_matrix_file
def _format_arg(self, name, spec, value):
if name == "bvecs":
new_val = self._create_gradient_matrix(self.inputs.bvecs, self.inputs.bvals)
return super(HARDIMat, self)._format_arg("bvecs", spec, new_val)
return super(HARDIMat, self)._format_arg(name, spec, value)
def _list_outputs(self):
outputs = self.output_spec().get()
outputs["out_file"] = os.path.abspath(self.inputs.out_file)
return outputs
class ODFReconInputSpec(CommandLineInputSpec):
DWI = File(
desc="Input raw data", argstr="%s", exists=True, mandatory=True, position=1
)
n_directions = traits.Int(
desc="Number of directions", argstr="%s", mandatory=True, position=2
)
n_output_directions = traits.Int(
desc="Number of output directions", argstr="%s", mandatory=True, position=3
)
out_prefix = traits.Str(
"odf", desc="Output file prefix", argstr="%s", usedefault=True, position=4
)
matrix = File(
argstr="-mat %s",
exists=True,
desc="""use given file as reconstruction matrix.""",
mandatory=True,
)
n_b0 = traits.Int(
argstr="-b0 %s",
desc="""\
number of b0 scans. by default the program gets this information
from the number of directions and number of volumes in
the raw data. useful when dealing with incomplete raw
data set or only using part of raw data set to reconstruct""",
mandatory=True,
)
output_type = traits.Enum(
"nii",
"analyze",
"ni1",
"nii.gz",
argstr="-ot %s",
desc="output file type",
usedefault=True,
)
sharpness = traits.Float(
desc="""\
smooth or sharpen the raw data. factor > 0 is smoothing.
factor < 0 is sharpening. default value is 0
NOTE: this option applies to DSI study only""",
argstr="-s %f",
)
filter = traits.Bool(
desc="""apply a filter (e.g. high pass) to the raw image""", argstr="-f"
)
subtract_background = traits.Bool(
desc="""subtract the background value before reconstruction""", argstr="-bg"
)
dsi = traits.Bool(desc="""indicates that the data is dsi""", argstr="-dsi")
output_entropy = traits.Bool(desc="""output entropy map""", argstr="-oe")
image_orientation_vectors = traits.List(
traits.Float(),
minlen=6,
maxlen=6,
desc="""\
specify image orientation vectors. if just one argument given,
will treat it as filename and read the orientation vectors from
the file. if 6 arguments are given, will treat them as 6 float
numbers and construct the 1st and 2nd vector and calculate the 3rd
one automatically.
this information will be used to determine image orientation,
as well as to adjust gradient vectors with oblique angle when""",
argstr="-iop %f",
)
oblique_correction = traits.Bool(
desc="""\
when oblique angle(s) applied, some SIEMENS dti protocols do not
adjust gradient accordingly, thus it requires adjustment for correct
diffusion tensor calculation""",
argstr="-oc",
)
class ODFReconOutputSpec(TraitedSpec):
B0 = File(exists=True)
DWI = File(exists=True)
max = File(exists=True)
ODF = File(exists=True)
entropy = File()
class ODFRecon(CommandLine):
"""Use odf_recon to generate tensors and other maps
"""
input_spec = ODFReconInputSpec
output_spec = ODFReconOutputSpec
_cmd = "odf_recon"
def _list_outputs(self):
out_prefix = self.inputs.out_prefix
output_type = self.inputs.output_type
outputs = self.output_spec().get()
outputs["B0"] = os.path.abspath(
fname_presuffix("", prefix=out_prefix, suffix="_b0." + output_type)
)
outputs["DWI"] = os.path.abspath(
fname_presuffix("", prefix=out_prefix, suffix="_dwi." + output_type)
)
outputs["max"] = os.path.abspath(
fname_presuffix("", prefix=out_prefix, suffix="_max." + output_type)
)
outputs["ODF"] = os.path.abspath(
fname_presuffix("", prefix=out_prefix, suffix="_odf." + output_type)
)
if isdefined(self.inputs.output_entropy):
outputs["entropy"] = os.path.abspath(
fname_presuffix("", prefix=out_prefix, suffix="_entropy." + output_type)
)
return outputs
class ODFTrackerInputSpec(CommandLineInputSpec):
max = File(exists=True, mandatory=True)
ODF = File(exists=True, mandatory=True)
input_data_prefix = traits.Str(
"odf", desc="recon data prefix", argstr="%s", usedefault=True, position=0
)
out_file = File(
"tracks.trk", desc="output track file", argstr="%s", usedefault=True, position=1
)
input_output_type = traits.Enum(
"nii",
"analyze",
"ni1",
"nii.gz",
argstr="-it %s",
desc="input and output file type",
usedefault=True,
)
runge_kutta2 = traits.Bool(
argstr="-rk2",
desc="""\
use 2nd order Runge-Kutta method for tracking.
default tracking method is non-interpolate streamline""",
)
step_length = traits.Float(
argstr="-l %f",
desc="""\
set step length, in the unit of minimum voxel size.
default value is 0.1.""",
)
angle_threshold = traits.Float(
argstr="-at %f",
desc="""\
set angle threshold. default value is 35 degree for
default tracking method and 25 for rk2""",
)
random_seed = traits.Int(
argstr="-rseed %s",
desc="""\
use random location in a voxel instead of the center of the voxel
to seed. can also define number of seed per voxel. default is 1""",
)
invert_x = traits.Bool(argstr="-ix", desc="invert x component of the vector")
invert_y = traits.Bool(argstr="-iy", desc="invert y component of the vector")
invert_z = traits.Bool(argstr="-iz", desc="invert z component of the vector")
swap_xy = traits.Bool(argstr="-sxy", desc="swap x and y vectors while tracking")
swap_yz = traits.Bool(argstr="-syz", desc="swap y and z vectors while tracking")
swap_zx = traits.Bool(argstr="-szx", desc="swap x and z vectors while tracking")
disc = traits.Bool(argstr="-disc", desc="use disc tracking")
mask1_file = File(
desc="first mask image", mandatory=True, argstr="-m %s", position=2
)
mask1_threshold = traits.Float(
desc="threshold value for the first mask image, if not given, the program will "
"try automatically find the threshold",
position=3,
)
mask2_file = File(desc="second mask image", argstr="-m2 %s", position=4)
mask2_threshold = traits.Float(
desc="threshold value for the second mask image, if not given, the program will "
"try automatically find the threshold",
position=5,
)
limit = traits.Int(
argstr="-limit %d",
desc="""\
in some special case, such as heart data, some track may go into
infinite circle and take long time to stop. this option allows
setting a limit for the longest tracking steps (voxels)""",
)
dsi = traits.Bool(
argstr="-dsi",
desc="""\
specify the input odf data is dsi. because dsi recon uses fixed
pre-calculated matrix, some special orientation patch needs to
be applied to keep dti/dsi/q-ball consistent.""",
)
image_orientation_vectors = traits.List(
traits.Float(),
minlen=6,
maxlen=6,
desc="""\
specify image orientation vectors. if just one argument given,
will treat it as filename and read the orientation vectors from
the file. if 6 arguments are given, will treat them as 6 float
numbers and construct the 1st and 2nd vector and calculate the 3rd
one automatically.
this information will be used to determine image orientation,
as well as to adjust gradient vectors with oblique angle when""",
argstr="-iop %f",
)
slice_order = traits.Int(
argstr="-sorder %d",
desc="set the slice order. 1 means normal, -1 means reversed. default value is 1",
)
voxel_order = traits.Enum(
"RAS",
"RPS",
"RAI",
"RPI",
"LAI",
"LAS",
"LPS",
"LPI",
argstr="-vorder %s",
desc="""\
specify the voxel order in RL/AP/IS (human brain) reference. must be
3 letters with no space in between.
for example, RAS means the voxel row is from L->R, the column
is from P->A and the slice order is from I->S.
by default voxel order is determined by the image orientation
(but NOT guaranteed to be correct because of various standards).
for example, siemens axial image is LPS, coronal image is LIP and
sagittal image is PIL.
this information also is NOT needed for tracking but will be saved
in the track file and is essential for track display to map onto
the right coordinates""",
)
class ODFTrackerOutputSpec(TraitedSpec):
track_file = File(exists=True, desc="output track file")
class ODFTracker(CommandLine):
"""Use odf_tracker to generate track file
"""
input_spec = ODFTrackerInputSpec
output_spec = ODFTrackerOutputSpec
_cmd = "odf_tracker"
def _run_interface(self, runtime):
_, _, ext = split_filename(self.inputs.max)
copyfile(
self.inputs.max,
os.path.abspath(self.inputs.input_data_prefix + "_max" + ext),
copy=False,
)
_, _, ext = split_filename(self.inputs.ODF)
copyfile(
self.inputs.ODF,
os.path.abspath(self.inputs.input_data_prefix + "_odf" + ext),
copy=False,
)
return super(ODFTracker, self)._run_interface(runtime)
def _list_outputs(self):
outputs = self.output_spec().get()
outputs["track_file"] = os.path.abspath(self.inputs.out_file)
return outputs