A new way to select b0 reference. Close #941

Author: oesteban

nipype/workflows/dmri/fsl/artifacts.py

@@ -342,17 +342,15 @@ def hmc_pipeline(name='motion_correct'):
     from nipype.workflows.data import get_flirt_schedule
 
     params = dict(dof=6, bgvalue=0, save_log=True, no_search=True,
-                  #cost='mutualinfo', cost_func='mutualinfo', bins=64,
+                  # cost='mutualinfo', cost_func='mutualinfo', bins=64,
                   schedule=get_flirt_schedule('hmc'))
 
     inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'ref_num',
                         'in_bvec', 'in_bval', 'in_mask']), name='inputnode')
-    refid = pe.Node(niu.IdentityInterface(fields=['ref_num']),
-                    name='RefVolume')
-    pick_ref = pe.Node(fsl.ExtractROI(t_size=1), name='GetRef')
-    pick_mov = pe.Node(niu.Function(input_names=['in_file', 'in_bval', 'volid'],
-                       output_names=['out_file', 'out_bval'],
-                       function=remove_comp), name='GetMovingDWs')
+    split = pe.Node(niu.Function(function=hmc_split,
+                    input_names=['in_file', 'in_bval', 'ref_num'],
+                    output_names=['out_ref', 'out_mov', 'out_bval', 'volid']),
+                    name='SplitDWI')
     flirt = dwi_flirt(flirt_param=params)
     insmat = pe.Node(niu.Function(input_names=['inlist', 'volid'],
                      output_names=['out'], function=insert_mat),
@@ -366,18 +364,15 @@ def hmc_pipeline(name='motion_correct'):
 
     wf = pe.Workflow(name=name)
     wf.connect([
-        (inputnode,  refid,      [(('ref_num', _checkrnum), 'ref_num')]),
-        (inputnode,  pick_ref,   [('in_file', 'in_file')]),
-        (inputnode,  pick_mov,   [('in_file', 'in_file'),
-                                  ('in_bval', 'in_bval')]),
+        (inputnode,     split,   [('in_file', 'in_file'),
+                                  ('in_bval', 'in_bval'),
+                                  ('ref_num', 'ref_num')]),
         (inputnode,  flirt,      [('in_mask', 'inputnode.ref_mask')]),
-        (refid,      pick_ref,   [('ref_num', 't_min')]),
-        (refid,      pick_mov,   [('ref_num', 'volid')]),
-        (pick_mov,   flirt,      [('out_file', 'inputnode.in_file'),
+        (split,      flirt,      [('out_ref', 'inputnode.reference'),
+                                  ('out_mov', 'inputnode.in_file'),
                                   ('out_bval', 'inputnode.in_bval')]),
-        (pick_ref,   flirt,      [('roi_file', 'inputnode.reference')]),
         (flirt,      insmat,     [('outputnode.out_xfms', 'inlist')]),
-        (refid,      insmat,     [('ref_num', 'volid')]),
+        (split,      insmat,     [('volid', 'volid')]),
         (inputnode,  rot_bvec,   [('in_bvec', 'in_bvec')]),
         (insmat,     rot_bvec,   [('out', 'in_matrix')]),
         (rot_bvec,   outputnode, [('out_file', 'out_bvec')]),

nipype/workflows/dmri/fsl/utils.py

@@ -229,6 +229,52 @@ def extract_bval(in_dwi, in_bval, b=0, out_file=None):
     return out_file
 
 
+def hmc_split(in_file, in_bval, ref_num=0, lowbval=5.0):
+    """
+    Selects the reference and moving volumes from a dwi dataset
+    for the purpose of HMC.
+    """
+    import numpy as np
+    import nibabel as nb
+    import os.path as op
+    from nipype.interfaces.base import isdefined
+
+    im = nb.load(in_file)
+    data = im.get_data()
+    hdr = im.get_header().copy()
+    bval = np.loadtxt(in_bval)
+
+    lowbs = np.where(bval <= lowbval)[0]
+
+    volid = lowbs[0]
+    if (isdefined(ref_num) and (ref_num < len(lowbs))):
+        volid = [ref_num]
+
+    if volid == 0:
+        data = data[..., 1:]
+        bval = bval[1:]
+    elif volid == (data.shape[-1] - 1):
+        data = data[..., :-1]
+        bval = bval[:-1]
+    else:
+        data = np.concatenate((data[..., :volid], data[..., (volid + 1):]),
+                              axis=3)
+        bval = np.hstack((bval[:volid], bval[(volid + 1):]))
+
+    out_ref = op.abspath('hmc_ref.nii.gz')
+    out_mov = op.abspath('hmc_mov.nii.gz')
+    out_bval = op.abspath('bval_split.txt')
+
+    hdr.set_data_shape(refdata.shape)
+    refdata = data[..., volid]
+    nb.Nifti1Image(refdata, im.get_affine(), hdr).to_filename(out_ref)
+
+    hdr.set_data_shape(data.shape)
+    nb.Nifti1Image(data, im.get_affine(), hdr).to_filename(out_mov)
+    np.savetxt(out_bval, bval)
+    return [out_ref, out_mov, out_bval, volid]
+
+
 def remove_comp(in_file, in_bval, volid=0, out_file=None):
     """
     Removes the volume ``volid`` from the 4D nifti file

nipype/workflows/fmri/spm/preprocess.py

@@ -254,11 +254,11 @@ def create_DARTEL_template(name='dartel_template'):
                              iterfield=['channel_files'],
                              name='segment')
     workflow.connect(inputnode, 'structural_files', segment, 'channel_files')
-    
+
     version = spm.Info.version()
     if version:
         spm_path = version['path']
-        if version['name'] == 'SPM8':          
+        if version['name'] == 'SPM8':
             tissue1 = ((os.path.join(spm_path,'toolbox/Seg/TPM.nii'), 1), 2, (True,True), (False, False))
             tissue2 = ((os.path.join(spm_path,'toolbox/Seg/TPM.nii'), 2), 2, (True,True), (False, False))
             tissue3 = ((os.path.join(spm_path,'toolbox/Seg/TPM.nii'), 3), 2, (True,False), (False, False))
@@ -275,11 +275,11 @@ def create_DARTEL_template(name='dartel_template'):
             tissue6 = ((os.path.join(spm_path,'tpm/TPM.nii'), 6), 2, (False,False), (False, False))
         else:
             logger.critical('Unsupported version of SPM')
-            
+
         segment.inputs.tissues = [tissue1, tissue2, tissue3, tissue4, tissue5, tissue6]
     else:
         logger.critical('SPM not found')
-        
+
     dartel = pe.Node(spm.DARTEL(), name='dartel')
 
     """Get the gray and white segmentation classes generated by NewSegment