-
Notifications
You must be signed in to change notification settings - Fork 59
/
scil_compute_fodf.py
executable file
·134 lines (104 loc) · 4.41 KB
/
scil_compute_fodf.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Script to compute Constrained Spherical Deconvolution (CSD) fiber ODFs.
By default, will output all possible files, using default names. Specific names
can be specified using the file flags specified in the "File flags" section.
If --not_all is set, only the files specified explicitly by the flags
will be output.
See [Tournier et al. NeuroImage 2007] and [Cote et al Tractometer MedIA 2013]
for quantitative comparisons with Sharpening Deconvolution Transform (SDT)
"""
import argparse
import logging
from dipy.io.gradients import read_bvals_bvecs
from dipy.direction.peaks import reshape_peaks_for_visualization
import nibabel as nib
import numpy as np
from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
assert_outputs_exist, add_force_b0_arg,
add_sh_basis_args)
from scilpy.reconst.fodf import compute_fodf
def _build_arg_parser():
p = argparse.ArgumentParser(
description=__doc__, formatter_class=argparse.RawTextHelpFormatter)
p.add_argument('input',
help='Path of the input diffusion volume.')
p.add_argument('bvals',
help='Path of the bvals file, in FSL format.')
p.add_argument('bvecs',
help='Path of the bvecs file, in FSL format.')
p.add_argument('frf_file',
help='Path of the FRF file')
p.add_argument(
'--sh_order', metavar='int', default=8, type=int,
help='SH order used for the CSD. (Default: 8)')
p.add_argument(
'--mask', metavar='',
help='Path to a binary mask. Only the data inside the mask will be '
'used for computations and reconstruction.')
p.add_argument(
'--processes', dest='nbr_processes', metavar='NBR', type=int,
help='Number of sub processes to start. Default : cpu count')
p.add_argument(
'--not_all', action='store_true',
help='If set, only saves the files specified using the file flags. '
'(Default: False)')
add_force_b0_arg(p)
add_sh_basis_args(p)
g = p.add_argument_group(title='File flags')
g.add_argument(
'--fodf', metavar='file', default='',
help='Output filename for the fiber ODF coefficients.')
g.add_argument(
'--peaks', metavar='file', default='',
help='Output filename for the extracted peaks.')
g.add_argument(
'--peak_indices', metavar='file', default='',
help='Output filename for the generated peaks indices on the sphere.')
add_overwrite_arg(p)
return p
def main():
parser = _build_arg_parser()
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
# Checking args
if not args.not_all:
args.fodf = args.fodf or 'fodf.nii.gz'
args.peaks = args.peaks or 'peaks.nii.gz'
args.peak_indices = args.peak_indices or 'peak_indices.nii.gz'
arglist = [args.fodf, args.peaks, args.peak_indices]
if args.not_all and not any(arglist):
parser.error('When using --not_all, you need to specify at least '
'one file to output.')
assert_inputs_exist(parser, [args.input, args.bvals, args.bvecs,
args.frf_file])
assert_outputs_exist(parser, args, arglist)
full_frf = np.loadtxt(args.frf_file)
vol = nib.load(args.input)
data = vol.dataobj
bvals, bvecs = read_bvals_bvecs(args.bvals, args.bvecs)
if args.mask is None:
mask = None
else:
mask = np.asanyarray(nib.load(args.mask).dataobj).astype(np.bool)
# Computing fODF
peaks_csd = compute_fodf(data, bvals, bvecs, full_frf,
sh_order=args.sh_order,
nbr_processes=args.nbr_processes,
mask=mask, sh_basis=args.sh_basis,
return_sh=True,
force_b0_threshold=args.force_b0_threshold)
# Saving results
if args.fodf:
nib.save(nib.Nifti1Image(peaks_csd.shm_coeff.astype(np.float32),
vol.affine), args.fodf)
if args.peaks:
nib.save(nib.Nifti1Image(
reshape_peaks_for_visualization(peaks_csd), vol.affine),
args.peaks)
if args.peak_indices:
nib.save(nib.Nifti1Image(peaks_csd.peak_indices, vol.affine),
args.peak_indices)
if __name__ == "__main__":
main()