1
+ def run_preAFQ (dwi_file , bvec_file , working_dir , sink_dir ):
2
+ import nipype .interfaces .fsl as fsl
3
+ from nipype .interfaces .fsl .utils import AvScale
4
+ import nipype .pipeline .engine as pe
5
+ import nipype .interfaces .utility as niu
6
+ import nipype .interfaces .io as nio
7
+ from nipype .workflows .dmri .fsl .epi import create_dmri_preprocessing
8
+ from nipype .interfaces .dipy import DTI
9
+ from nipype .utils .filemanip import fname_presuffix
10
+ import os
11
+
12
+ wf = create_dmri_preprocessing (name = 'preAFQ' ,
13
+ use_fieldmap = False ,
14
+ fieldmap_registration = False )
15
+ wf .inputs .inputnode .ref_num = 0
16
+ wf .inputs .inputnode .in_file = dwi_file
17
+ wf .inputs .inputnode .in_bvec = bvec_file
18
+
19
+ dwi_filename = os .path .split (dwi_file )[1 ].split (".nii.gz" )[0 ]
20
+ bids_subject_name = dwi_filename .split ("_" )[0 ]
21
+ assert bids_subject_name .startswith ("sub-" )
22
+
23
+ inputnode = wf .get_node ("inputnode" )
24
+ outputspec = wf .get_node ("outputnode" )
25
+
26
+ #QC: FLIRT translation and rotation parameters
27
+ flirt = wf .get_node ("motion_correct.coregistration" )
28
+ #flirt.inputs.save_mats = True
29
+
30
+
31
+ #wf.connect(getmotion, "motion_params", datasink, "dti.@motparams")
32
+
33
+ def get_flirt_motion_parameters (flirt_out_mats ):
34
+ import nipype .interfaces .fsl .utils as fsl
35
+ import os
36
+ import numpy as np
37
+
38
+ def get_params (A ):
39
+ """This is a copy of spm's spm_imatrix where
40
+ we already know the rotations and translations matrix,
41
+ shears and zooms (as outputs from fsl FLIRT/avscale)
42
+ Let A = the 4x4 rotation and translation matrix
43
+ R = [ c5*c6, c5*s6, s5]
44
+ [-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
45
+ [-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]
46
+ """
47
+ import numpy as np
48
+
49
+ def rang (b ):
50
+ a = min (max (b , - 1 ), 1 )
51
+ return a
52
+ Ry = np .arcsin (A [0 ,2 ])
53
+ #Rx = np.arcsin(A[1,2]/np.cos(Ry))
54
+ #Rz = np.arccos(A[0,1]/np.sin(Ry))
55
+
56
+ if (abs (Ry )- np .pi / 2 )** 2 < 1e-9 :
57
+ Rx = 0
58
+ Rz = np .arctan2 (- rang (A [1 ,0 ]), rang (- A [2 ,0 ]/ A [0 ,2 ]))
59
+ else :
60
+ c = np .cos (Ry )
61
+ Rx = np .arctan2 (rang (A [1 ,2 ]/ c ), rang (A [2 ,2 ]/ c ))
62
+ Rz = np .arctan2 (rang (A [0 ,1 ]/ c ), rang (A [0 ,0 ]/ c ))
63
+
64
+ rotations = [Rx , Ry , Rz ]
65
+ translations = [A [0 ,3 ], A [1 ,3 ], A [2 ,3 ]]
66
+
67
+ return rotations , translations
68
+
69
+
70
+ motion_params = open (os .path .abspath ('motion_parameters.par' ),'w' )
71
+ for mat in flirt_out_mats :
72
+ res = fsl .AvScale (mat_file = mat ).run ()
73
+ A = np .asarray (res .outputs .rotation_translation_matrix )
74
+ rotations , translations = get_params (A )
75
+ for i in rotations + translations :
76
+ motion_params .write ('%f ' % i )
77
+ motion_params .write ('\n ' )
78
+ motion_params .close ()
79
+ motion_params = os .path .abspath ('motion_parameters.par' )
80
+ return motion_params
81
+
82
+
83
+ getmotion = pe .Node (niu .Function (input_names = ["flirt_out_mats" ],
84
+ output_names = ["motion_params" ],
85
+ function = get_flirt_motion_parameters ),
86
+ name = "get_motion_parameters" ,iterfield = "flirt_out_mats" )
87
+
88
+ wf .connect (flirt , "out_matrix_file" , getmotion , "flirt_out_mats" )
89
+
90
+ from nipype .algorithms .rapidart import ArtifactDetect
91
+
92
+ art = pe .Node (interface = ArtifactDetect (), name = "art" )
93
+ art .inputs .use_differences = [True , True ]
94
+ art .inputs .save_plot = False
95
+ art .inputs .use_norm = True
96
+ art .inputs .norm_threshold = 3
97
+ art .inputs .zintensity_threshold = 9
98
+ art .inputs .mask_type = 'spm_global'
99
+ art .inputs .parameter_source = 'FSL'
100
+
101
+ wf .connect (getmotion , "motion_params" , art , "realignment_parameters" )
102
+ wf .connect (outputspec , "dmri_corrected" , art , "realigned_files" )
103
+
104
+ datasink = pe .Node (nio .DataSink (),name = "sinker" )
105
+ datasink .inputs .base_directory = sink_dir
106
+ datasink .inputs .substitutions = [("vol0000_flirt_merged.nii.gz" , dwi_filename + '_corrected.nii.gz' ),
107
+ ("stats.vol0000_flirt_merged.txt" , dwi_filename + ".motion.txt" ),
108
+ ("derivatives/dti" , "derivatives/{}/dti" .format (bids_subject_name ))]
109
+
110
+ wf .connect (art , "statistic_files" , datasink , "dti.@artstat" )
111
+ wf .connect (outputspec , "dmri_corrected" , datasink , "dti.@corrected" )
112
+ wf .connect (outputspec ,"bvec_rotated" , datasink , "dti.@rotated" )
113
+
114
+ wf .base_dir = working_dir
115
+ wf .run ()
116
+
117
+ from glob import glob
118
+ import os
119
+
120
+ dmri_corrected = glob (os .path .join (sink_dir , 'dti' , '*.nii.gz' ))[0 ]
121
+ bvec_rotated = glob (os .path .join (sink_dir , 'dti' , '*.bvec' ))[0 ]
122
+ motion_outliers = glob (os .path .join (sink_dir , 'dti' , '*.txt' ))[0 ]
123
+ return dmri_corrected , bvec_rotated , motion_outliers
0 commit comments