Merge pull request #548 from PennLINC/replace-fsl-split

Replace fsl split with non-FSL tools

qsiprep/interfaces/images.py

@@ -11,6 +11,8 @@
 
 import os
 from subprocess import Popen, PIPE
+import subprocess
+import glob
 from textwrap import indent
 import numpy as np
 import nibabel as nb
@@ -19,7 +21,9 @@
 from nipype import logging
 from nipype.utils.filemanip import fname_presuffix
 from nipype.interfaces.base import (isdefined, traits, TraitedSpec, BaseInterfaceInputSpec,
-                                    SimpleInterface, File, InputMultiObject, OutputMultiObject)
+                                    SimpleInterface, File, InputMultiObject, OutputMultiObject,
+                                    OutputMultiPath)
+from nipype.interfaces.afni.base import (AFNICommand, AFNICommandInputSpec, AFNICommandOutputSpec)
 from nipype.interfaces import fsl
 # from qsiprep.interfaces.images import (
 #    nii_ones_like, extract_wm, SignalExtraction, MatchHeader,
@@ -30,7 +34,45 @@
 LOGGER = logging.getLogger('nipype.interface')
 
 
-class SplitDWIsInputSpec(BaseInterfaceInputSpec):
+class SplitDWIsBvalsInputSpec(BaseInterfaceInputSpec):
+    split_files = InputMultiObject(desc='pre-split DWI images')
+    bvec_file = File(desc='the bvec file')
+    bval_file = File(desc='the bval file')
+    deoblique_bvecs = traits.Bool(False, usedefault=True,
+                                  desc='write LPS+ world coordinate bvecs')
+    b0_threshold = traits.Int(50, usedefault=True,
+                              desc='Maximum b-value that can be considered a b0')
+
+
+class SplitDWIsBvalsOutputSpec(TraitedSpec):
+    bval_files = OutputMultiObject(File(exists=True), desc='single volume bvals')
+    bvec_files = OutputMultiObject(File(exists=True), desc='single volume bvecs')
+    b0_images = OutputMultiObject(File(exists=True), desc='just the b0s')
+    b0_indices = traits.List(desc='list of original indices for each b0 image')
+
+
+class SplitDWIsBvals(SimpleInterface):
+    input_spec = SplitDWIsBvalsInputSpec
+    output_spec = SplitDWIsBvalsOutputSpec
+
+    def _run_interface(self, runtime):
+
+        split_bval_files, split_bvec_files = split_bvals_bvecs(
+            self.inputs.bval_file, self.inputs.bvec_file, 
+            self.inputs.split_files, self.inputs.deoblique_bvecs, 
+            runtime.cwd)
+
+        bvalues = np.loadtxt(self.inputs.bval_file)
+        b0_indices = np.flatnonzero(bvalues < self.inputs.b0_threshold)
+        b0_paths = [self.inputs.split_files[idx] for idx in b0_indices]
+        self._results['bval_files'] = split_bval_files
+        self._results['bvec_files'] = split_bvec_files
+        self._results['b0_images'] = b0_paths
+        self._results['b0_indices'] = b0_indices.tolist()
+
+        return runtime
+
+class SplitDWIsFSLInputSpec(BaseInterfaceInputSpec):
     dwi_file = File(desc='the dwi image')
     bvec_file = File(desc='the bvec file')
     bval_file = File(desc='the bval file')
@@ -40,17 +82,17 @@ class SplitDWIsInputSpec(BaseInterfaceInputSpec):
                               desc='Maximum b-value that can be considered a b0')
 
 
-class SplitDWIsOutputSpec(TraitedSpec):
+class SplitDWIsFSLOutputSpec(TraitedSpec):
     dwi_files = OutputMultiObject(File(exists=True), desc='single volume dwis')
     bval_files = OutputMultiObject(File(exists=True), desc='single volume bvals')
     bvec_files = OutputMultiObject(File(exists=True), desc='single volume bvecs')
     b0_images = OutputMultiObject(File(exists=True), desc='just the b0s')
     b0_indices = traits.List(desc='list of original indices for each b0 image')
 
 
-class SplitDWIs(SimpleInterface):
-    input_spec = SplitDWIsInputSpec
-    output_spec = SplitDWIsOutputSpec
+class SplitDWIsFSL(SimpleInterface):
+    input_spec = SplitDWIsFSLInputSpec
+    output_spec = SplitDWIsFSLOutputSpec
 
     def _run_interface(self, runtime):
         split = fsl.Split(dimension='t', in_file=self.inputs.dwi_file)
@@ -71,7 +113,6 @@ def _run_interface(self, runtime):
 
         return runtime
 
-
 def _flatten(in_list):
     out_list = []
     for item in in_list:
@@ -623,3 +664,52 @@ def to_lps(input_img, new_axcodes=("L", "P", "S")):
         return reoriented_img
     else:
         return input_img
+
+class TSplitInputSpec(AFNICommandInputSpec):
+    in_file = File(
+        desc="input file to 3dTsplit4D",
+        argstr=" %s",
+        position=-1,
+        mandatory=True,
+        copyfile=False,
+    )
+    out_name = File(
+        mandatory=True,
+        desc="output image file name",
+        argstr="-prefix %s.nii",
+    )
+    digits = traits.Int(
+        argstr="-digits %d", desc="Number of digits to include in split file names"
+    )
+
+class TSplitOutputSpec(TraitedSpec):
+    out_files = OutputMultiPath(File(exists=True))
+
+class TSplit(AFNICommand):
+    """Converts a 3D + time dataset into multiple 3D volumes (one volume per file).
+    For complete details, see the `3dTsplit4D Documentation.
+    <https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dTsplit4D.html>`_
+    """
+
+    _cmd = "3dTsplit4D"
+    input_spec = TSplitInputSpec
+    output_spec = TSplitOutputSpec
+
+    def _list_outputs(self):
+        """Create a Bunch which contains all possible files generated
+        by running the interface.  Some files are always generated, others
+        depending on which ``inputs`` options are set.
+        Returns
+        -------
+        outputs : Bunch object
+            Bunch object containing all possible files generated by
+            interface object.
+            If None, file was not generated
+            Else, contains path, filename of generated outputfile
+        """
+        outputs = self._outputs().get()
+        ext = '.nii'
+        outputs["out_files"] = sorted(glob.glob(
+            os.path.join(os.getcwd(),'{outname}.**.nii'.format(
+            outname=self.inputs.out_name))))
+        return outputs

qsiprep/workflows/dwi/fsl.py

@@ -17,7 +17,7 @@
 from .registration import init_b0_to_anat_registration_wf, init_direct_b0_acpc_wf
 from ...interfaces.eddy import (GatherEddyInputs, ExtendedEddy, Eddy2SPMMotion,
                                 boilerplate_from_eddy_config)
-from ...interfaces.images import SplitDWIs, ConformDwi, IntraModalMerge
+from ...interfaces.images import SplitDWIsFSL, ConformDwi, IntraModalMerge
 from ...interfaces.reports import TopupSummary
 from ...interfaces.nilearn import EnhanceB0
 from ...interfaces import DerivativesDataSink
@@ -145,7 +145,7 @@ def init_fsl_hmc_wf(scan_groups,
     # Convert eddy outputs back to LPS+, split them
     back_to_lps = pe.Node(ConformDwi(orientation="LPS"), name='back_to_lps')
     cnr_lps = pe.Node(ConformDwi(orientation="LPS"), name='cnr_lps')
-    split_eddy_lps = pe.Node(SplitDWIs(b0_threshold=b0_threshold, deoblique_bvecs=True),
+    split_eddy_lps = pe.Node(SplitDWIsFSL(b0_threshold=b0_threshold, deoblique_bvecs=True),
                              name="split_eddy_lps")
 
     # Convert the b=0 template from pre_eddy_b0_ref to LPS+

qsiprep/workflows/dwi/hmc_sdc.py

@@ -13,7 +13,7 @@
 
 from .registration import init_b0_to_anat_registration_wf, init_direct_b0_acpc_wf
 from ...interfaces.gradients import SliceQC, CombineMotions, GradientRotation
-from ...interfaces.images import SplitDWIs
+from ...interfaces.images import SplitDWIsBvals, TSplit
 from ..fieldmap.base import init_sdc_wf
 from ..fieldmap.drbuddi import init_drbuddi_wf
 from ...engine import Workflow
@@ -97,8 +97,8 @@ def init_qsiprep_hmcsdc_wf(scan_groups,
     workflow = Workflow(name=name)
 
     # Split the input data into single volumes, put bvecs in LPS+ world reference frame
-    split_dwis = pe.Node(SplitDWIs(b0_threshold=b0_threshold, deoblique_bvecs=True),
-                         name='split_dwis')
+    split_dwis = pe.Node(TSplit(digits=4, out_name='vol'),name='split_dwis')
+    split_bvals = pe.Node(SplitDWIsBvals(b0_threshold=b0_threshold, deoblique_bvecs=True),name='split_bvals')
 
     # Motion correct the data
     dwi_hmc_wf = init_dwi_hmc_wf(hmc_transform, hmc_model, hmc_align_to,
@@ -115,21 +115,25 @@ def init_qsiprep_hmcsdc_wf(scan_groups,
 
     workflow.connect([
         (inputnode, split_dwis, [
-            ('dwi_file', 'dwi_file'),
+            ('dwi_file', 'in_file')]),
+        (inputnode, split_bvals, [
             ('bval_file', 'bval_file'),
             ('bvec_file', 'bvec_file')]),
-        (split_dwis, dwi_hmc_wf, [
-            ('dwi_files', 'inputnode.dwi_files'),
+        (split_dwis, split_bvals, [
+            ('out_files','split_files')]),
+        (split_bvals, dwi_hmc_wf, [
             ('bval_files', 'inputnode.bvals'),
             ('bvec_files', 'inputnode.bvecs'),
             ('b0_images', 'inputnode.b0_images'),
             ('b0_indices', 'inputnode.b0_indices')]),
+        (split_dwis, dwi_hmc_wf, [
+            ('out_files', 'inputnode.dwi_files')]),
         (inputnode, dwi_hmc_wf, [
             ('t1_brain', 'inputnode.t1_brain'),
             ('t1_mask', 'inputnode.t1_mask'),
             ('t1_seg', 'inputnode.t1_seg'),
             ('original_files', 'inputnode.original_files')]),
-        (split_dwis, slice_qc, [('dwi_files', 'uncorrected_dwi_files')]),
+        (split_dwis, slice_qc, [('out_files', 'uncorrected_dwi_files')]),
         (dwi_hmc_wf, outputnode, [
             ('outputnode.final_template', 'pre_sdc_template'),
             (('outputnode.forward_transforms', _list_squeeze),
@@ -145,7 +149,7 @@ def init_qsiprep_hmcsdc_wf(scan_groups,
             ('outputnode.noise_free_dwis', 'ideal_image_files'),
             ('outputnode.final_template_mask', 'mask_image')]),
         (summarize_motion, outputnode, [('spm_motion_file', 'motion_params')]),
-        (split_dwis, outputnode, [
+        (split_bvals, outputnode, [
             ('bvec_files', 'bvec_files_to_transform'),
             ('bval_files', 'bval_files'),
             ('b0_indices', 'b0_indices')]),
@@ -180,9 +184,9 @@ def init_qsiprep_hmcsdc_wf(scan_groups,
                 (('outputnode.forward_transforms', _list_squeeze), 
                  'affine_transforms')]),
             (split_dwis, apply_hmc_transforms, [
-                ('dwi_files', 'input_image'),
-                ('dwi_files', 'reference_image')]),
-            (split_dwis, rotate_gradients, [
+                ('out_files', 'input_image'),
+                ('out_files', 'reference_image')]),
+            (split_bvals, rotate_gradients, [
                 ('bvec_files', 'bvec_files'),
                 ('bval_files', 'bval_files')]),
             (apply_hmc_transforms, drbuddi_wf, [