1
- def run_preAFQ (dwi_file , bvec_file , bval_file , subjects_dir , working_dir , sink_dir ):
2
- import nipype .interfaces .freesurfer as fs
3
- import nipype .interfaces .fsl as fsl
4
- from nipype .interfaces .fsl .utils import AvScale
5
- import nipype .pipeline .engine as pe
6
- import nipype .interfaces .utility as niu
7
- import nipype .interfaces .io as nio
8
- from nipype .workflows .dmri .fsl .epi import create_dmri_preprocessing
9
- from nipype .interfaces .dipy import DTI
10
- from nipype .utils .filemanip import fname_presuffix
11
- import os
12
- from nipype .interfaces .dipy import DTI
13
- from nipype .interfaces .c3 import C3dAffineTool
14
-
15
- wf = create_dmri_preprocessing (name = 'preAFQ' ,
16
- use_fieldmap = False ,
17
- fieldmap_registration = False )
1
+ import os .path as op
2
+ from glob import glob
3
+ from shutil import copyfile
4
+
5
+ import nibabel as nib
6
+ import nipype .interfaces .freesurfer as fs
7
+ import nipype .interfaces .fsl as fsl
8
+ import nipype .interfaces .io as nio
9
+ import nipype .interfaces .utility as niu
10
+ import nipype .pipeline .engine as pe
11
+ import numpy as np
12
+ from nipype .algorithms .rapidart import ArtifactDetect
13
+ from nipype .interfaces .dipy import DTI
14
+ from nipype .interfaces .fsl .utils import AvScale
15
+ from nipype .utils .filemanip import fname_presuffix
16
+ from nipype .workflows .dmri .fsl .epi import create_dmri_preprocessing
17
+
18
+
19
+ def run_preAFQ (dwi_file , bvec_file , bval_file ,
20
+ subjects_dir , working_dir , out_dir ):
21
+ wf = create_dmri_preprocessing (name = 'preAFQ' ,
22
+ use_fieldmap = False ,
23
+ fieldmap_registration = False )
18
24
wf .inputs .inputnode .ref_num = 0
19
25
wf .inputs .inputnode .in_file = dwi_file
20
26
wf .inputs .inputnode .in_bvec = bvec_file
21
27
22
- dwi_filename = os . path .split (dwi_file )[1 ].split (".nii.gz" )[0 ]
23
- bids_subject_name = dwi_filename .split ("_" )[0 ]
24
- assert bids_subject_name .startswith ("sub-" )
28
+ dwi_fname = op .split (dwi_file )[1 ].split (".nii.gz" )[0 ]
29
+ bids_sub_name = dwi_fname .split ("_" )[0 ]
30
+ assert bids_sub_name .startswith ("sub-" )
25
31
26
- inputnode = wf .get_node ("inputnode" )
32
+ # inputnode = wf.get_node("inputnode")
27
33
outputspec = wf .get_node ("outputnode" )
28
34
29
- #QC: FLIRT translation and rotation parameters
35
+ # QC: FLIRT translation and rotation parameters
30
36
flirt = wf .get_node ("motion_correct.coregistration" )
31
- #flirt.inputs.save_mats = True
37
+ # flirt.inputs.save_mats = True
32
38
33
39
get_tensor = pe .Node (DTI (), name = "dipy_tensor" )
34
- wf .connect (outputspec , "dmri_corrected" ,get_tensor ,"in_file" )
35
- #wf.connect(inputspec2,"bvals", get_tensor, "in_bval")
40
+ wf .connect (outputspec , "dmri_corrected" , get_tensor , "in_file" )
41
+ # wf.connect(inputspec2,"bvals", get_tensor, "in_bval")
36
42
get_tensor .inputs .in_bval = bval_file
37
43
wf .connect (outputspec , "bvec_rotated" , get_tensor , "in_bvec" )
38
44
@@ -44,40 +50,39 @@ def run_preAFQ(dwi_file, bvec_file, bval_file, subjects_dir, working_dir, sink_d
44
50
fslroi = pe .Node (fsl .ExtractROI (t_min = 0 , t_size = 1 ), name = "fslroi" )
45
51
wf .connect (outputspec , "dmri_corrected" , fslroi , "in_file" )
46
52
47
- bbreg = pe .Node (fs .BBRegister (contrast_type = "t2" , init = "fsl" , out_fsl_file = True ,
48
- subjects_dir = subjects_dir ,
49
- epi_mask = True ),
50
- name = "bbreg" )
51
- #wf.connect(inputspec2,"fsid", bbreg,"subject_id")
52
- bbreg .inputs .subject_id = 'freesurfer' #bids_subject_name
53
- wf .connect (fslroi ,"roi_file" , bbreg , "source_file" )
53
+ bbreg = pe .Node (fs .BBRegister (contrast_type = "t2" , init = "fsl" ,
54
+ out_fsl_file = True , subjects_dir = subjects_dir ,
55
+ epi_mask = True ), name = "bbreg" )
56
+ # wf.connect(inputspec2,"fsid", bbreg,"subject_id")
57
+ bbreg .inputs .subject_id = 'freesurfer' # bids_sub_name
58
+ wf .connect (fslroi , "roi_file" , bbreg , "source_file" )
54
59
55
60
voltransform = pe .Node (fs .ApplyVolTransform (inverse = True ),
56
- iterfield = ['source_file' , 'reg_file' ],
57
- name = 'transform' )
61
+ iterfield = ['source_file' , 'reg_file' ],
62
+ name = 'transform' )
58
63
voltransform .inputs .subjects_dir = subjects_dir
59
64
60
65
vt2 = voltransform .clone ("transform_aparcaseg" )
61
66
vt2 .inputs .interp = "nearest"
62
67
63
68
def binarize_aparc (aparc_aseg ):
64
- import nibabel as nib
65
- from nipype .utils .filemanip import fname_presuffix
66
- import os
67
69
img = nib .load (aparc_aseg )
68
70
data , aff = img .get_data (), img .affine
69
- outfile = fname_presuffix (aparc_aseg , suffix = "bin.nii.gz" , newpath = os .path .abspath ("." ), use_ext = False )
71
+ outfile = fname_presuffix (
72
+ aparc_aseg , suffix = "bin.nii.gz" ,
73
+ newpath = op .abspath ("." ), use_ext = False
74
+ )
70
75
nib .Nifti1Image ((data > 0 ).astype (float ), aff ).to_filename (outfile )
71
76
return outfile
72
77
73
- #wf.connect(inputspec2, "mask_nii", voltransform, "target_file")
74
- create_mask = pe .Node (niu .Function (input_names = ["aparc_aseg" ],
75
- output_names = ["outfile" ],
78
+ # wf.connect(inputspec2, "mask_nii", voltransform, "target_file")
79
+ create_mask = pe .Node (niu .Function (input_names = ["aparc_aseg" ],
80
+ output_names = ["outfile" ],
76
81
function = binarize_aparc ),
77
- name = "bin_aparc" )
82
+ name = "bin_aparc" )
78
83
79
84
def get_aparc_aseg (subjects_dir , sub ):
80
- return os . path .join (subjects_dir , sub , "mri" , "aparc+aseg.mgz" )
85
+ return op .join (subjects_dir , sub , "mri" , "aparc+aseg.mgz" )
81
86
82
87
create_mask .inputs .aparc_aseg = get_aparc_aseg (subjects_dir , 'freesurfer' )
83
88
wf .connect (create_mask , "outfile" , voltransform , "target_file" )
@@ -86,25 +91,20 @@ def get_aparc_aseg(subjects_dir, sub):
86
91
wf .connect (bbreg , "out_reg_file" , voltransform , "reg_file" )
87
92
88
93
vt2 .inputs .target_file = get_aparc_aseg (subjects_dir , 'freesurfer' )
89
- #wf.connect(inputspec2, "aparc_aseg", vt2, "target_file")
94
+ # wf.connect(inputspec2, "aparc_aseg", vt2, "target_file")
90
95
wf .connect (fslroi , "roi_file" , vt2 , "source_file" )
91
96
wf .connect (bbreg , "out_reg_file" , vt2 , "reg_file" )
92
97
93
- # AK (2017): THIS NODE MIGHT NOT BE NECESSARY
98
+ # AK (2017): THIS NODE MIGHT NOT BE NECESSARY
94
99
# AK (2018) doesn't know why she said that above..
95
100
threshold2 = pe .Node (fs .Binarize (min = 0.5 , out_type = 'nii.gz' , dilate = 1 ),
96
- iterfield = ['in_file' ],
97
- name = 'threshold2' )
101
+ iterfield = ['in_file' ],
102
+ name = 'threshold2' )
98
103
wf .connect (voltransform , "transformed_file" , threshold2 , "in_file" )
99
104
100
-
101
- #wf.connect(getmotion, "motion_params", datasink, "dti.@motparams")
105
+ # wf.connect(getmotion, "motion_params", datasink, "dti.@motparams")
102
106
103
107
def get_flirt_motion_parameters (flirt_out_mats ):
104
- import nipype .interfaces .fsl .utils as fsl
105
- import os
106
- import numpy as np
107
-
108
108
def get_params (A ):
109
109
"""This is a copy of spm's spm_imatrix where
110
110
we already know the rotations and translations matrix,
@@ -114,54 +114,51 @@ def get_params(A):
114
114
[-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
115
115
[-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]
116
116
"""
117
- import numpy as np
118
-
119
117
def rang (b ):
120
118
a = min (max (b , - 1 ), 1 )
121
119
return a
122
- Ry = np .arcsin (A [0 ,2 ])
123
- #Rx = np.arcsin(A[1,2]/ np.cos(Ry))
124
- #Rz = np.arccos(A[0,1]/ np.sin(Ry))
120
+ Ry = np .arcsin (A [0 , 2 ])
121
+ # Rx = np.arcsin(A[1, 2] / np.cos(Ry))
122
+ # Rz = np.arccos(A[0, 1] / np.sin(Ry))
125
123
126
124
if (abs (Ry )- np .pi / 2 )** 2 < 1e-9 :
127
125
Rx = 0
128
- Rz = np .arctan2 (- rang (A [1 ,0 ]), rang (- A [2 ,0 ]/ A [0 ,2 ]))
126
+ Rz = np .arctan2 (- rang (A [1 , 0 ]), rang (- A [2 , 0 ]/ A [0 , 2 ]))
129
127
else :
130
- c = np .cos (Ry )
131
- Rx = np .arctan2 (rang (A [1 ,2 ]/ c ), rang (A [2 ,2 ]/ c ))
132
- Rz = np .arctan2 (rang (A [0 ,1 ]/ c ), rang (A [0 ,0 ]/ c ))
128
+ c = np .cos (Ry )
129
+ Rx = np .arctan2 (rang (A [1 , 2 ]/ c ), rang (A [2 , 2 ]/ c ))
130
+ Rz = np .arctan2 (rang (A [0 , 1 ]/ c ), rang (A [0 , 0 ]/ c ))
133
131
134
132
rotations = [Rx , Ry , Rz ]
135
- translations = [A [0 ,3 ], A [1 ,3 ], A [2 ,3 ]]
133
+ translations = [A [0 , 3 ], A [1 , 3 ], A [2 , 3 ]]
136
134
137
135
return rotations , translations
138
136
139
-
140
- motion_params = open (os .path .abspath ('motion_parameters.par' ),'w' )
137
+ motion_params = open (op .abspath ('motion_parameters.par' ), 'w' )
141
138
for mat in flirt_out_mats :
142
- res = fsl . AvScale (mat_file = mat ).run ()
139
+ res = AvScale (mat_file = mat ).run ()
143
140
A = np .asarray (res .outputs .rotation_translation_matrix )
144
141
rotations , translations = get_params (A )
145
142
for i in rotations + translations :
146
- motion_params .write ('%f ' % i )
143
+ motion_params .write ('%f ' % i )
147
144
motion_params .write ('\n ' )
148
145
motion_params .close ()
149
- motion_params = os . path .abspath ('motion_parameters.par' )
146
+ motion_params = op .abspath ('motion_parameters.par' )
150
147
return motion_params
151
148
152
-
153
- getmotion = pe .Node (niu .Function (input_names = ["flirt_out_mats" ],
154
- output_names = ["motion_params" ],
155
- function = get_flirt_motion_parameters ),
156
- name = "get_motion_parameters" ,iterfield = "flirt_out_mats" )
149
+ getmotion = pe .Node (
150
+ niu .Function (input_names = ["flirt_out_mats" ],
151
+ output_names = ["motion_params" ],
152
+ function = get_flirt_motion_parameters ),
153
+ name = "get_motion_parameters" ,
154
+ iterfield = "flirt_out_mats"
155
+ )
157
156
158
157
wf .connect (flirt , "out_matrix_file" , getmotion , "flirt_out_mats" )
159
158
160
- from nipype .algorithms .rapidart import ArtifactDetect
161
-
162
159
art = pe .Node (interface = ArtifactDetect (), name = "art" )
163
160
art .inputs .use_differences = [True , True ]
164
- art .inputs .save_plot = False
161
+ art .inputs .save_plot = False
165
162
art .inputs .use_norm = True
166
163
art .inputs .norm_threshold = 3
167
164
art .inputs .zintensity_threshold = 9
@@ -171,25 +168,27 @@ def rang(b):
171
168
wf .connect (getmotion , "motion_params" , art , "realignment_parameters" )
172
169
wf .connect (outputspec , "dmri_corrected" , art , "realigned_files" )
173
170
174
- datasink = pe .Node (nio .DataSink (),name = "sinker" )
175
- datasink .inputs .base_directory = sink_dir
176
- datasink .inputs .substitutions = [("vol0000_flirt_merged.nii.gz" , dwi_filename + '.nii.gz' ),
177
- ("stats.vol0000_flirt_merged.txt" , dwi_filename + ".art.json" ),
178
- ("motion_parameters.par" , dwi_filename + ".motion.txt" ),
179
- ("_rotated.bvec" , ".bvec" ),
180
- ("aparc+aseg_warped_out" , dwi_filename .replace ("_dwi" , "_aparc+aseg" )),
181
- ("art.vol0000_flirt_merged_outliers.txt" , dwi_filename + ".outliers.txt" ),
182
- ("vol0000_flirt_merged" , dwi_filename ),
183
- ("_roi_bbreg_freesurfer" , "_register" ),
184
- ("aparc+asegbin_warped_thresh" , dwi_filename .replace ("_dwi" , "_mask" )),
185
- ("derivatives/preafq" , "derivatives/{}/preafq" .format (bids_subject_name ))]
171
+ datasink = pe .Node (nio .DataSink (), name = "sinker" )
172
+ datasink .inputs .base_directory = out_dir
173
+ datasink .inputs .substitutions = [
174
+ ("vol0000_flirt_merged.nii.gz" , dwi_fname + '.nii.gz' ),
175
+ ("stats.vol0000_flirt_merged.txt" , dwi_fname + ".art.json" ),
176
+ ("motion_parameters.par" , dwi_fname + ".motion.txt" ),
177
+ ("_rotated.bvec" , ".bvec" ),
178
+ ("aparc+aseg_warped_out" , dwi_fname .replace ("_dwi" , "_aparc+aseg" )),
179
+ ("art.vol0000_flirt_merged_outliers.txt" , dwi_fname + ".outliers.txt" ),
180
+ ("vol0000_flirt_merged" , dwi_fname ),
181
+ ("_roi_bbreg_freesurfer" , "_register" ),
182
+ ("aparc+asegbin_warped_thresh" , dwi_fname .replace ("_dwi" , "_mask" )),
183
+ ("derivatives/preafq" , "derivatives/{}/preafq" .format (bids_sub_name ))
184
+ ]
186
185
187
186
wf .connect (art , "statistic_files" , datasink , "preafq.art.@artstat" )
188
187
wf .connect (art , "outlier_files" , datasink , "preafq.art.@artoutlier" )
189
188
wf .connect (outputspec , "dmri_corrected" , datasink , "preafq.dwi.@corrected" )
190
- wf .connect (outputspec ,"bvec_rotated" , datasink , "preafq.dwi.@rotated" )
189
+ wf .connect (outputspec , "bvec_rotated" , datasink , "preafq.dwi.@rotated" )
191
190
wf .connect (getmotion , "motion_params" , datasink , "preafq.art.@motion" )
192
-
191
+
193
192
wf .connect (get_tensor , "out_file" , datasink , "preafq.dti.@tensor" )
194
193
wf .connect (get_tensor , "fa_file" , datasink , "preafq.dti.@fa" )
195
194
wf .connect (get_tensor , "md_file" , datasink , "preafq.dti.@md" )
@@ -198,10 +197,10 @@ def rang(b):
198
197
wf .connect (get_tensor , "out_file" , datasink , "preafq.dti.@scaled_tensor" )
199
198
200
199
wf .connect (bbreg , "min_cost_file" , datasink , "preafq.reg.@mincost" )
201
- wf .connect (bbreg ,"out_fsl_file" , datasink , "preafq.reg.@fslfile" )
200
+ wf .connect (bbreg , "out_fsl_file" , datasink , "preafq.reg.@fslfile" )
202
201
wf .connect (bbreg , "out_reg_file" , datasink , "preafq.reg.@reg" )
203
202
wf .connect (threshold2 , "binary_file" , datasink , "preafq.anat.@mask" )
204
- #wf.connect(vt2, "transformed_file", datasink, "dwi.@aparc_aseg")
203
+ # wf.connect(vt2, "transformed_file", datasink, "dwi.@aparc_aseg")
205
204
206
205
convert = pe .Node (fs .MRIConvert (out_type = "niigz" ), name = "convert2nii" )
207
206
wf .connect (vt2 , "transformed_file" , convert , "in_file" )
@@ -210,16 +209,15 @@ def rang(b):
210
209
wf .base_dir = working_dir
211
210
wf .run ()
212
211
213
- from glob import glob
214
- import os
215
- from shutil import copyfile
216
-
217
- copyfile (bval_file , os .path .join (sink_dir , bids_subject_name , "preafq" , "dwi" , os .path .split (bval_file )[1 ]))
218
-
219
- dmri_corrected = glob (os .path .join (sink_dir , '*/preafq/dwi' , '*.nii.gz' ))[0 ]
220
- bvec_rotated = glob (os .path .join (sink_dir , '*/preafq/dwi' , '*.bvec' ))[0 ]
221
- bval_file = glob (os .path .join (sink_dir , '*/preafq/dwi' , '*.bval' ))[0 ]
222
- art_file = glob (os .path .join (sink_dir , '*/preafq/art' , '*.art.json' ))[0 ]
223
- motion_file = glob (os .path .join (sink_dir , '*/preafq/art' , '*.motion.txt' ))[0 ]
224
- outlier_file = glob (os .path .join (sink_dir , '*/preafq/art' , '*.outliers.txt' ))[0 ]
225
- return dmri_corrected , bvec_rotated , art_file , motion_file , outlier_file
212
+ copyfile (bval_file , op .join (
213
+ out_dir , bids_sub_name , "preafq" , "dwi" ,
214
+ op .split (bval_file )[1 ]
215
+ ))
216
+
217
+ dmri_corrected = glob (op .join (out_dir , '*/preafq/dwi' , '*.nii.gz' ))[0 ]
218
+ bvec_rotated = glob (op .join (out_dir , '*/preafq/dwi' , '*.bvec' ))[0 ]
219
+ bval_file = glob (op .join (out_dir , '*/preafq/dwi' , '*.bval' ))[0 ]
220
+ art_file = glob (op .join (out_dir , '*/preafq/art' , '*.art.json' ))[0 ]
221
+ motion_file = glob (op .join (out_dir , '*/preafq/art' , '*.motion.txt' ))[0 ]
222
+ outlier_file = glob (op .join (out_dir , '*/preafq/art' , '*.outliers.txt' ))[0 ]
223
+ return dmri_corrected , bvec_rotated , art_file , motion_file , outlier_file
0 commit comments