Merge pull request #60 from josephmje/develop

enh/fix: add eddy_cuda and minor updates to pybids and eddy quad

Author: richford

README.rst

@@ -31,11 +31,11 @@ You should have raw data organized in the BIDS format. Also, you should have run
 |   |-- sub-01
 |       |-- dwi
 |           |-- sub-01_ses-01_dwi.nii.gz
-|           |-- sub-01_ses-01_dwi.bvals
-|           |-- sub-01_ses-01_dwi.bvecs
+|           |-- sub-01_ses-01_dwi.bval
+|           |-- sub-01_ses-01_dwi.bvec
 |       |-- fmap
-|           |-- sub-01_ses-01_dir-AP_dwi.nii.gz
-|           |-- sub-01_ses-01_dir-PA_dwi.nii.gz
+|           |-- sub-01_ses-01_acq-dwi_dir-AP_epi.nii.gz
+|           |-- sub-01_ses-01_acq-dwi_dir-PA_epi.nii.gz
 |   |-- derivatives
 |       |-- sub-01
 |           |-- freesurfer

dmriprep/cli.py

@@ -1,13 +1,15 @@
 # -*- coding: utf-8 -*-
 
 """Console script for dmriprep."""
+import os
 import sys
+import warnings
+
 import click
-from . import run
+
 from . import io
+from . import run
 from .data import get_dataset
-import os
-import warnings
 
 # Filter warnings that are visible whenever you import another package that
 # was compiled against an older numpy than is installed.

dmriprep/io.py

@@ -2,9 +2,10 @@
 BIDS-functions to return inputs for the run.py functions.
 
 """
-from glob import glob
 import os
 import os.path as op
+from glob import glob
+
 import bids
 
 
@@ -20,18 +21,26 @@ def get_bids_subject_input_files(subject_id, bids_input_directory):
     subjects = layout.get_subjects()
     assert subject_id in subjects, "subject {} is not in the bids folder".format(subject_id)
 
-    ap_file = layout.get(subject=subject_id, fmap="epi", modality="fmap", dir="AP")
+    ap_file = layout.get(subject=subject_id,
+                         datatype='fmap',
+                         suffix='epi',
+                         dir='AP',
+                         extensions=['.nii', '.nii.gz'])
     assert len(ap_file) == 1, 'found {} ap fieldmap files and we need just 1'.format(len(ap_file))
 
-    pa_file = layout.get(subject=subject_id, fmap="epi", modality="fmap", dir="PA")
+    pa_file = layout.get(subject=subject_id,
+                         datatype='fmap',
+                         suffix='epi',
+                         dir='PA',
+                         extensions=['.nii', '.nii.gz'])
     assert len(pa_file) == 1, 'found {} pa fieldmap files and we need just 1'.format(len(pa_file))
 
-    dwi_files = layout.get(subject=subject_id, modality="dwi")
+    dwi_files = layout.get(subject=subject_id, datatype='dwi', suffix='dwi')
     valid_dwi_files = []
 
     for d in dwi_files:
-        if d.filename.startswith(op.abspath(op.join(bids_input_directory, 'sub-' + subject_id))):
-            valid_dwi_files.append(d.filename)
+        if d.path.startswith(op.abspath(op.join(bids_input_directory, 'sub-' + subject_id))):
+            valid_dwi_files.append(d.path)
 
     dwi_file = [d for d in valid_dwi_files if d.endswith('.nii.gz') and not "TRACE" in d]
     assert len(dwi_file) == 1, 'found {} dwi files and we need just 1'.format(len(dwi_file))
@@ -51,8 +60,8 @@ def get_bids_subject_input_files(subject_id, bids_input_directory):
 
     outputs = dict(subject_id="sub-"+subject_id,
                    dwi_file=dwi_file[0],
-                   dwi_file_AP=ap_file[0].filename,
-                   dwi_file_PA=pa_file[0].filename,
+                   dwi_file_AP=ap_file[0].path,
+                   dwi_file_PA=pa_file[0].path,
                    bvec_file=bvec_file[0],
                    bval_file=bval_file[0],
                    subjects_dir=op.abspath(subjects_dir))

dmriprep/qc.py

@@ -1,13 +1,15 @@
-import nibabel as nib
-import numpy as np
+import base64
 import os.path as op
-from dipy.segment.mask import median_otsu
 from io import BytesIO
-from nipype.utils.filemanip import save_json, load_json
-import base64
+
 import matplotlib
 matplotlib.use('agg')
 import matplotlib.pyplot as plt
+import nibabel as nib
+import numpy as np
+
+from dipy.segment.mask import median_otsu
+from nipype.utils.filemanip import save_json, load_json
 
 
 def reorient_array(data, aff):

dmriprep/run.py

@@ -1,5 +1,5 @@
-import os.path as op
 import os
+import os.path as op
 from shutil import copyfile
 
 
@@ -363,6 +363,12 @@ def get_dmriprep_pe_workflow():
     eddy.inputs.residuals = True
     import multiprocessing
     eddy.inputs.num_threads = multiprocessing.cpu_count()
+    import numba.cuda
+    try:
+        if numba.cuda.gpus:
+            eddy.inputs.use_cuda = True
+    except:
+        eddy.inputs.use_cuda = False
 
     topup = prep.get_node('peb_correction.topup')
     topup.inputs.checksize = True
@@ -374,12 +380,14 @@ def get_dmriprep_pe_workflow():
 
     eddy_quad = pe.Node(fsl.EddyQuad(verbose=True, checksize=True), name="eddy_quad")
     get_path = lambda x: x.split('.nii.gz')[0].split('_fix')[0]
-    wf.connect(prep, ('fsl_eddy.out_corrected', get_path), eddy_quad, "base_name")
+    get_qc_path = lambda x: x.split('.nii.gz')[0] + '.qc'
+    wf.connect(prep, ('fsl_eddy.out_corrected', get_path), eddy_quad, 'base_name')
     wf.connect(inputspec, 'bval_file', eddy_quad, 'bval_file')
     wf.connect(prep, 'Rotate_Bvec.out_file', eddy_quad, 'bvec_file')
     wf.connect(prep, 'peb_correction.topup.out_field', eddy_quad, 'field')
     wf.connect(prep, 'gen_index.out_file', eddy_quad, 'idx_file')
     wf.connect(prep, 'peb_correction.topup.out_enc_file', eddy_quad, 'param_file')
+    wf.connect(prep, ('fsl_eddy.out_corrected', get_qc_path), eddy_quad, 'output_dir')
 
     # need a mask file for eddy_quad. lets get it from the B0.
     def get_b0_mask_fn(b0_file):
@@ -396,7 +404,6 @@ def get_b0_mask_fn(b0_file):
         return mask_file
 
 
-
     def id_outliers_fn(outlier_report, threshold, dwi_file):
         """Get list of scans that exceed threshold for number of outliers
 
@@ -596,6 +603,7 @@ def apply_transform_to_bvecs_fn(bvec_file, reg_mat_file):
         import numpy as np
         import nipype.utils.filemanip as fm
         import os
+
         aff = np.loadtxt(reg_mat_file)
         bvecs = np.loadtxt(bvec_file)
         bvec_trans = []
@@ -654,6 +662,7 @@ def binarize_aparc(aparc_aseg):
         import nibabel as nib
         from nipype.utils.filemanip import fname_presuffix
         import os.path as op
+
         img = nib.load(aparc_aseg)
         data, aff = img.get_data(), img.affine
         outfile = fname_presuffix(
@@ -746,13 +755,13 @@ def binarize_aparc(aparc_aseg):
     wf.connect(scale_tensor_eddy, "out_file", datasink, "dmriprep.dti_eddy.@scaled_tensor")
 
     # all the eddy_quad stuff
-    wf.connect(eddy_quad, 'out_qc_json', datasink, "dmriprep.qc.@eddyquad_json")
-    wf.connect(eddy_quad, 'out_qc_pdf', datasink, "dmriprep.qc.@eddyquad_pdf")
-    wf.connect(eddy_quad, 'out_avg_b_png', datasink, "dmriprep.qc.@eddyquad_bpng")
-    wf.connect(eddy_quad, 'out_avg_b0_pe_png', datasink, "dmriprep.qc.@eddyquad_b0png")
-    wf.connect(eddy_quad, 'out_cnr_png', datasink, "dmriprep.qc.@eddyquad_cnr")
-    wf.connect(eddy_quad, 'out_vdm_png', datasink, "dmriprep.qc.@eddyquad_vdm")
-    wf.connect(eddy_quad, 'out_residuals', datasink, 'dmriprep.qc.@eddyquad_resid')
+    wf.connect(eddy_quad, 'qc_json', datasink, "dmriprep.qc.@eddyquad_json")
+    wf.connect(eddy_quad, 'qc_pdf', datasink, "dmriprep.qc.@eddyquad_pdf")
+    wf.connect(eddy_quad, 'avg_b_png', datasink, "dmriprep.qc.@eddyquad_bpng")
+    wf.connect(eddy_quad, 'avg_b0_pe_png', datasink, "dmriprep.qc.@eddyquad_b0png")
+    wf.connect(eddy_quad, 'cnr_png', datasink, "dmriprep.qc.@eddyquad_cnr")
+    wf.connect(eddy_quad, 'vdm_png', datasink, "dmriprep.qc.@eddyquad_vdm")
+    wf.connect(eddy_quad, 'residuals', datasink, 'dmriprep.qc.@eddyquad_resid')
 
     # anatomical registration stuff
     wf.connect(bbreg, "min_cost_file", datasink, "dmriprep.reg.@mincost")
@@ -784,7 +793,7 @@ def report_fn(dwi_corrected_file, eddy_rms, eddy_report,
 
     # for the report, lets show the eddy corrected (full volume) image
     wf.connect(voltransform, "transformed_file", report_node, 'dwi_corrected_file')
-    wf.connect(eddy_quad, 'out_qc_json', report_node, 'eddy_qc_file')
+    wf.connect(eddy_quad, 'qc_json', report_node, 'eddy_qc_file')
 
     # add the rms movement output from eddy
     wf.connect(prep, "fsl_eddy.out_movement_rms", report_node, 'eddy_rms')
@@ -802,6 +811,7 @@ def report_fn(dwi_corrected_file, eddy_rms, eddy_report,
     # its super annoying.
     def name_files_nicely(dwi_file, subject_id):
         import os.path as op
+
         dwi_fname = op.split(dwi_file)[1].split(".nii.gz")[0]
         substitutions = [
             ("vol0000_flirt_merged.nii.gz", dwi_fname + '.nii.gz'),

setup.py

@@ -19,6 +19,10 @@
     'pandas',
     'parse',
     'tqdm',
+    'pybids',
+    'matplotlib',
+    'cytoolz',
+    'numba'
 ]
 
 setup_requirements = ['pytest-runner', ]