Skip to content
This repository was archived by the owner on Dec 27, 2022. It is now read-only.

Commit 85ddb1c

Browse files
arokemarokem
committed
Remove cruft, move better version of this into run.py.
1 parent b4c9605 commit 85ddb1c

File tree

2 files changed

+114
-337
lines changed

2 files changed

+114
-337
lines changed

‎dmriprep/run.py

Lines changed: 114 additions & 126 deletions
Original file line numberDiff line numberDiff line change
@@ -1,62 +1,70 @@
11
import os.path as op
2-
from glob import glob
32
from shutil import copyfile
43

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_dmriprep(dwi_file, bvec_file, bval_file,
4+
5+
def run_dmriprep(dwi_file, dwi_file_AP, dwi_file_PA, bvec_file, bval_file,
206
subjects_dir, working_dir, out_dir):
21-
wf = create_dmri_preprocessing(name='dmriprep',
22-
use_fieldmap=False,
23-
fieldmap_registration=False)
24-
wf.inputs.inputnode.ref_num = 0
25-
wf.inputs.inputnode.in_file = dwi_file
26-
wf.inputs.inputnode.in_bvec = bvec_file
7+
"""This is for HBN diffusion data
8+
9+
Assuming phase-encode direction is AP/PA for topup
10+
"""
2711

12+
import nipype.interfaces.freesurfer as fs
13+
import nipype.interfaces.fsl as fsl
14+
import nipype.interfaces.io as nio
15+
import nipype.interfaces.utility as niu
16+
import nipype.pipeline.engine as pe
17+
from nipype.interfaces.dipy import DTI
18+
from nipype.workflows.dmri.fsl.artifacts import all_fsl_pipeline
19+
20+
# some bookkeeping (getting the filename, gettings the BIDS subject name)
2821
dwi_fname = op.split(dwi_file)[1].split(".nii.gz")[0]
2922
bids_sub_name = dwi_fname.split("_")[0]
3023
assert bids_sub_name.startswith("sub-")
3124

32-
# inputnode = wf.get_node("inputnode")
33-
outputspec = wf.get_node("outputnode")
25+
# Grab the preprocessing all_fsl_pipeline
26+
# AK: watch out, other datasets might be encoded LR
27+
epi_AP = {'echospacing': 66.5e-3, 'enc_dir': 'y-'}
28+
epi_PA = {'echospacing': 66.5e-3, 'enc_dir': 'y'}
29+
prep = all_fsl_pipeline(epi_params=epi_AP, altepi_params=epi_PA)
30+
31+
# initialize an overall workflow
32+
wf = pe.Workflow(name="dmriprep")
33+
wf.base_dir = op.abspath(working_dir)
3434

35-
# QC: FLIRT translation and rotation parameters
36-
flirt = wf.get_node("motion_correct.coregistration")
37-
# flirt.inputs.save_mats = True
35+
prep.inputs.inputnode.in_file = dwi_file
36+
# prep.inputs.inputnode.alt_file = dwi_file_PA
37+
prep.inputs.inputnode.in_bvec = bvec_file
38+
prep.inputs.inputnode.in_bval = bval_file
39+
eddy = prep.get_node('fsl_eddy')
40+
eddy.inputs.repol = True
41+
eddy.inputs.niter = 1 # TODO: change back to 5 when running for real
42+
43+
merge = pe.Node(fsl.Merge(dimension='t'), name="mergeAPPA")
44+
merge.inputs.in_files = [dwi_file_AP, dwi_file_PA]
45+
wf.connect(merge, 'merged_file', prep, 'inputnode.alt_file')
46+
47+
fslroi = pe.Node(fsl.ExtractROI(t_min=0, t_size=1), name="fslroi")
48+
wf.connect(prep, "outputnode.out_file", fslroi, "in_file")
49+
50+
bbreg = pe.Node(fs.BBRegister(contrast_type="t2", init="coreg",
51+
out_fsl_file=True,
52+
subjects_dir=subjects_dir,
53+
epi_mask=True),
54+
name="bbreg")
55+
bbreg.inputs.subject_id = 'freesurfer' # bids_sub_name
56+
wf.connect(fslroi, "roi_file", bbreg, "source_file")
3857

3958
get_tensor = pe.Node(DTI(), name="dipy_tensor")
40-
wf.connect(outputspec, "dmri_corrected", get_tensor, "in_file")
41-
# wf.connect(inputspec2,"bvals", get_tensor, "in_bval")
59+
wf.connect(prep, "outputnode.out_file", get_tensor, "in_file")
4260
get_tensor.inputs.in_bval = bval_file
43-
wf.connect(outputspec, "bvec_rotated", get_tensor, "in_bvec")
61+
wf.connect(prep, "outputnode.out_bvec", get_tensor, "in_bvec")
4462

4563
scale_tensor = pe.Node(name='scale_tensor', interface=fsl.BinaryMaths())
4664
scale_tensor.inputs.operation = 'mul'
4765
scale_tensor.inputs.operand_value = 1000
4866
wf.connect(get_tensor, 'out_file', scale_tensor, 'in_file')
4967

50-
fslroi = pe.Node(fsl.ExtractROI(t_min=0, t_size=1), name="fslroi")
51-
wf.connect(outputspec, "dmri_corrected", fslroi, "in_file")
52-
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")
59-
6068
voltransform = pe.Node(fs.ApplyVolTransform(inverse=True),
6169
iterfield=['source_file', 'reg_file'],
6270
name='transform')
@@ -65,7 +73,12 @@ def run_dmriprep(dwi_file, bvec_file, bval_file,
6573
vt2 = voltransform.clone("transform_aparcaseg")
6674
vt2.inputs.interp = "nearest"
6775

76+
vt3 = voltransform.clone("transform_orig")
77+
6878
def binarize_aparc(aparc_aseg):
79+
import nibabel as nib
80+
from nipype.utils.filemanip import fname_presuffix
81+
import os.path as op
6982
img = nib.load(aparc_aseg)
7083
data, aff = img.get_data(), img.affine
7184
outfile = fname_presuffix(
@@ -75,7 +88,6 @@ def binarize_aparc(aparc_aseg):
7588
nib.Nifti1Image((data > 0).astype(float), aff).to_filename(outfile)
7689
return outfile
7790

78-
# wf.connect(inputspec2, "mask_nii", voltransform, "target_file")
7991
create_mask = pe.Node(niu.Function(input_names=["aparc_aseg"],
8092
output_names=["outfile"],
8193
function=binarize_aparc),
@@ -84,140 +96,116 @@ def binarize_aparc(aparc_aseg):
8496
def get_aparc_aseg(subjects_dir, sub):
8597
return op.join(subjects_dir, sub, "mri", "aparc+aseg.mgz")
8698

99+
def get_orig(subjects_dir, sub):
100+
return op.join(subjects_dir, sub, "mri", "orig.mgz")
101+
87102
create_mask.inputs.aparc_aseg = get_aparc_aseg(subjects_dir, 'freesurfer')
88103
wf.connect(create_mask, "outfile", voltransform, "target_file")
89104

90105
wf.connect(fslroi, "roi_file", voltransform, "source_file")
91106
wf.connect(bbreg, "out_reg_file", voltransform, "reg_file")
92107

93108
vt2.inputs.target_file = get_aparc_aseg(subjects_dir, 'freesurfer')
94-
# wf.connect(inputspec2, "aparc_aseg", vt2, "target_file")
95109
wf.connect(fslroi, "roi_file", vt2, "source_file")
96110
wf.connect(bbreg, "out_reg_file", vt2, "reg_file")
97111

112+
vt3.inputs.target_file = get_orig(subjects_dir, 'freesurfer')
113+
wf.connect(fslroi, "roi_file", vt3, "source_file")
114+
wf.connect(bbreg, "out_reg_file", vt3, "reg_file")
115+
98116
# AK (2017): THIS NODE MIGHT NOT BE NECESSARY
99117
# AK (2018) doesn't know why she said that above..
100118
threshold2 = pe.Node(fs.Binarize(min=0.5, out_type='nii.gz', dilate=1),
101119
iterfield=['in_file'],
102120
name='threshold2')
103121
wf.connect(voltransform, "transformed_file", threshold2, "in_file")
104122

105-
# wf.connect(getmotion, "motion_params", datasink, "dti.@motparams")
106-
107-
def get_flirt_motion_parameters(flirt_out_mats):
108-
def get_params(A):
109-
"""This is a copy of spm's spm_imatrix where
110-
we already know the rotations and translations matrix,
111-
shears and zooms (as outputs from fsl FLIRT/avscale)
112-
Let A = the 4x4 rotation and translation matrix
113-
R = [ c5*c6, c5*s6, s5]
114-
[-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
115-
[-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]
116-
"""
117-
def rang(b):
118-
a = min(max(b, -1), 1)
119-
return a
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))
123-
124-
if (abs(Ry)-np.pi/2)**2 < 1e-9:
125-
Rx = 0
126-
Rz = np.arctan2(-rang(A[1, 0]), rang(-A[2, 0]/A[0, 2]))
127-
else:
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))
131-
132-
rotations = [Rx, Ry, Rz]
133-
translations = [A[0, 3], A[1, 3], A[2, 3]]
134-
135-
return rotations, translations
136-
137-
motion_params = open(op.abspath('motion_parameters.par'), 'w')
138-
for mat in flirt_out_mats:
139-
res = AvScale(mat_file=mat).run()
140-
A = np.asarray(res.outputs.rotation_translation_matrix)
141-
rotations, translations = get_params(A)
142-
for i in rotations+translations:
143-
motion_params.write('%f ' % i)
144-
motion_params.write('\n')
145-
motion_params.close()
146-
motion_params = op.abspath('motion_parameters.par')
147-
return motion_params
148-
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-
)
156-
157-
wf.connect(flirt, "out_matrix_file", getmotion, "flirt_out_mats")
158-
159-
art = pe.Node(interface=ArtifactDetect(), name="art")
160-
art.inputs.use_differences = [True, True]
161-
art.inputs.save_plot = False
162-
art.inputs.use_norm = True
163-
art.inputs.norm_threshold = 3
164-
art.inputs.zintensity_threshold = 9
165-
art.inputs.mask_type = 'spm_global'
166-
art.inputs.parameter_source = 'FSL'
167-
168-
wf.connect(getmotion, "motion_params", art, "realignment_parameters")
169-
wf.connect(outputspec, "dmri_corrected", art, "realigned_files")
170-
171123
datasink = pe.Node(nio.DataSink(), name="sinker")
172124
datasink.inputs.base_directory = out_dir
173125
datasink.inputs.substitutions = [
174126
("vol0000_flirt_merged.nii.gz", dwi_fname + '.nii.gz'),
175127
("stats.vol0000_flirt_merged.txt", dwi_fname + ".art.json"),
176128
("motion_parameters.par", dwi_fname + ".motion.txt"),
177129
("_rotated.bvec", ".bvec"),
178-
("aparc+aseg_warped_out", dwi_fname.replace("_dwi", "_aparc+aseg")),
179130
("art.vol0000_flirt_merged_outliers.txt", dwi_fname + ".outliers.txt"),
180131
("vol0000_flirt_merged", dwi_fname),
181132
("_roi_bbreg_freesurfer", "_register"),
133+
("dwi/eddy_corrected", "dwi/%s" % dwi_fname),
134+
("dti/eddy_corrected", "dti/%s" % dwi_fname.replace("_dwi", "")),
135+
("reg/eddy_corrected", "reg/%s" % dwi_fname.replace("_dwi", "")),
182136
("aparc+asegbin_warped_thresh", dwi_fname.replace("_dwi", "_mask")),
137+
("aparc+aseg_warped_out", dwi_fname.replace("_dwi", "_aparc+aseg")),
138+
("art.eddy_corrected_outliers", dwi_fname.replace("dwi", "outliers")),
139+
("color_fa", "colorfa"),
140+
("orig_warped_out", dwi_fname.replace("_dwi", "_T1w")),
141+
# ("eddy_corrected_", dwi_fname.replace("dwi", "")),
142+
("stats.eddy_corrected", dwi_fname.replace("dwi", "artStats")),
143+
("eddy_corrected.eddy_parameters", dwi_fname+".eddy_parameters"),
144+
("qc/eddy_corrected", "qc/"+dwi_fname),
183145
("derivatives/dmriprep", "derivatives/{}/dmriprep".format(bids_sub_name))
184146
]
185147

186-
wf.connect(art, "statistic_files", datasink, "dmriprep.art.@artstat")
187-
wf.connect(art, "outlier_files", datasink, "dmriprep.art.@artoutlier")
188-
wf.connect(outputspec, "dmri_corrected", datasink, "dmriprep.dwi.@corrected")
189-
wf.connect(outputspec, "bvec_rotated", datasink, "dmriprep.dwi.@rotated")
190-
wf.connect(getmotion, "motion_params", datasink, "dmriprep.art.@motion")
148+
wf.connect(prep, "outputnode.out_file", datasink, "dmriprep.dwi.@corrected")
149+
wf.connect(prep, "outputnode.out_bvec", datasink, "dmriprep.dwi.@rotated")
150+
wf.connect(prep, "fsl_eddy.out_parameter",
151+
datasink, "dmriprep.qc.@eddyparams")
152+
153+
wf.connect(prep, "fsl_eddy.out_movement_rms",
154+
datasink, "dmriprep.qc.@eddyparamsrms")
155+
wf.connect(prep, "fsl_eddy.out_outlier_report",
156+
datasink, "dmriprep.qc.@eddyparamsreport")
157+
wf.connect(prep, "fsl_eddy.out_restricted_movement_rms",
158+
datasink, "dmriprep.qc.@eddyparamsrestrictrms")
159+
wf.connect(prep, "fsl_eddy.out_shell_alignment_parameters",
160+
datasink, "dmriprep.qc.@eddyparamsshellalign")
191161

192162
wf.connect(get_tensor, "out_file", datasink, "dmriprep.dti.@tensor")
193163
wf.connect(get_tensor, "fa_file", datasink, "dmriprep.dti.@fa")
194164
wf.connect(get_tensor, "md_file", datasink, "dmriprep.dti.@md")
195165
wf.connect(get_tensor, "ad_file", datasink, "dmriprep.dti.@ad")
196166
wf.connect(get_tensor, "rd_file", datasink, "dmriprep.dti.@rd")
167+
wf.connect(get_tensor, "color_fa_file", datasink, "dmriprep.dti.@colorfa")
197168
wf.connect(get_tensor, "out_file", datasink, "dmriprep.dti.@scaled_tensor")
198169

199170
wf.connect(bbreg, "min_cost_file", datasink, "dmriprep.reg.@mincost")
200171
wf.connect(bbreg, "out_fsl_file", datasink, "dmriprep.reg.@fslfile")
201172
wf.connect(bbreg, "out_reg_file", datasink, "dmriprep.reg.@reg")
202173
wf.connect(threshold2, "binary_file", datasink, "dmriprep.anat.@mask")
203-
# wf.connect(vt2, "transformed_file", datasink, "dwi.@aparc_aseg")
204174

205175
convert = pe.Node(fs.MRIConvert(out_type="niigz"), name="convert2nii")
206176
wf.connect(vt2, "transformed_file", convert, "in_file")
207177
wf.connect(convert, "out_file", datasink, "dmriprep.anat.@aparc_aseg")
208178

209-
wf.base_dir = working_dir
179+
convert1 = convert.clone("convertorig2nii")
180+
wf.connect(vt3, "transformed_file", convert1, "in_file")
181+
wf.connect(convert1, "out_file", datasink, "dmriprep.anat.@anat")
182+
183+
def reportNodeFunc(dwi_corrected_file, eddy_rms, eddy_report,
184+
color_fa_file, anat_mask_file):
185+
from dmriprep.qc import create_report_json
186+
187+
report = create_report_json(dwi_corrected_file, eddy_rms, eddy_report,
188+
color_fa_file, anat_mask_file)
189+
return report
190+
191+
reportNode = pe.Node(niu.Function(
192+
input_names=['dwi_corrected_file', 'eddy_rms',
193+
'eddy_report', 'color_fa_file',
194+
'anat_mask_file'],
195+
output_names=['report'],
196+
function=reportNodeFunc
197+
), name="reportJSON")
198+
199+
wf.connect(prep, "outputnode.out_file", reportNode, 'dwi_corrected_file')
200+
wf.connect(prep, "fsl_eddy.out_movement_rms", reportNode, 'eddy_rms')
201+
wf.connect(prep, "fsl_eddy.out_outlier_report", reportNode, 'eddy_report')
202+
wf.connect(threshold2, "binary_file", reportNode, 'anat_mask_file')
203+
wf.connect(get_tensor, "color_fa_file", reportNode, 'color_fa_file')
204+
205+
wf.connect(reportNode, 'report', datasink, 'dmriprep.report.@report')
206+
210207
wf.run()
211208

212209
copyfile(bval_file, op.join(
213-
out_dir, bids_sub_name, "dmriprep", "dwi",
214-
op.split(bval_file)[1]
210+
out_dir, "dmriprep", "dwi", op.split(bval_file)[1]
215211
))
216-
217-
dmri_corrected = glob(op.join(out_dir, '*/dmriprep/dwi', '*.nii.gz'))[0]
218-
bvec_rotated = glob(op.join(out_dir, '*/dmriprep/dwi', '*.bvec'))[0]
219-
bval_file = glob(op.join(out_dir, '*/dmriprep/dwi', '*.bval'))[0]
220-
art_file = glob(op.join(out_dir, '*/dmriprep/art', '*.art.json'))[0]
221-
motion_file = glob(op.join(out_dir, '*/dmriprep/art', '*.motion.txt'))[0]
222-
outlier_file = glob(op.join(out_dir, '*/dmriprep/art', '*.outliers.txt'))[0]
223-
return dmri_corrected, bvec_rotated, art_file, motion_file, outlier_file

0 commit comments

Comments
 (0)