tested new XFibres. new parallel workflow

Author: oesteban

nipype/interfaces/fsl/base.py

@@ -255,7 +255,7 @@ class FSLXCommandInputSpec(FSLCommandInputSpec):
     bvals = File(exists=True, argstr='--bvals=%s', mandatory=True,
                  desc='b values file')
 
-    logdir = Directory('logdir', argstr='--logdir=%s', usedefault=True)
+    logdir = Directory('.', argstr='--logdir=%s', usedefault=True)
     n_fibres = traits.Range(low=1, argstr='--nfibres=%d', desc=('Maximum '
                             'number of fibres to fit in each voxel'))
     model = traits.Enum(1, 2, argstr='--model=%d',
@@ -311,34 +311,22 @@ class FSLXCommandInputSpec(FSLCommandInputSpec):
 
 
 class FSLXCommandOutputSpec(TraitedSpec):
-    merged_thsamples = OutputMultiPath(File(), desc=('Samples from '
-                                       'the distribution on theta'))
-    merged_phsamples = OutputMultiPath(File(), desc=('Samples from '
-                                       'the distribution on phi'))
-    merged_fsamples = OutputMultiPath(File(),
-                                      desc=('Samples from the distribution on '
-                                            'anisotropic volume fraction.'))
-
-    mean_thsamples = OutputMultiPath(File(), desc=('Mean of '
-                                     'distribution on theta'))
-    mean_phsamples = OutputMultiPath(File(), desc=('Mean of '
-                                     'distribution on phi'))
+    dsamples = File(desc=('Samples from the distribution on diffusivity d'))
+    d_stdsamples = File(desc=('Std of samples from the distribution d'))
+    dyads = OutputMultiPath(File(), desc=('Mean of PDD distribution'
+                            ' in vector form.'))
+    fsamples = OutputMultiPath(File(), desc=('Samples from the '
+                               'distribution on f anisotropy'))
+    mean_dsamples = File(desc='Mean of distribution on diffusivity d')
+    mean_d_stdsamples = File(desc='Mean of distribution on diffusivity d')
     mean_fsamples = OutputMultiPath(File(), desc=('Mean of '
                                     'distribution on f anisotropy'))
-
-    mean_dsamples = File(desc='Mean of distribution on '
-                         'diffusivity d')
     mean_S0samples = File(desc='Mean of distribution on T2w'
                           'baseline signal intensity S0')
     mean_tausamples = File(desc='Mean of distribution on '
                            'tau samples (only with rician noise)')
-
-    dyads = OutputMultiPath(File(), desc=('Mean of PDD distribution'
-                            ' in vector form.'))
-    dyads_disp = OutputMultiPath(File(), desc=('Uncertainty on the '
-                                 ' estimated fiber orientation'))
-    fsamples = OutputMultiPath(File(), desc=('Samples from the '
-                               'distribution on anisotropic volume fraction'))
+    phsamples = OutputMultiPath(File(), desc=('phi samples, per fiber'))
+    thsamples = OutputMultiPath(File(), desc=('theta samples, per fiber'))
 
 
 class FSLXCommand(FSLCommand):
@@ -359,45 +347,38 @@ def _list_outputs(self):
         outputs = self.output_spec().get()
         out_dir = self._out_dir
 
-        for k in outputs.keys():
-            if k not in ('outputtype', 'environ', 'args', 'bpx_out_directory',
-                         'xfms_directory', 'mean_dsamples', 'mean_S0samples',
-                         'mean_tausamples'):
-                outputs[k] = []
+        if isdefined(self.inputs.logdir):
+            out_dir = os.path.abspath(self.inputs.logdir)
+        else:
+            out_dir = os.path.abspath('logdir')
+
+        multi_out = ['dyads', 'fsamples', 'mean_fsamples',
+                     'phsamples', 'thsamples']
+        single_out = ['dsamples', 'd_stdsamples', 'mean_dsamples',
+                      'mean_S0samples', 'mean_d_stdsamples']
 
-        outputs['mean_dsamples'] = self._gen_fname('mean_dsamples',
-                                                   cwd=out_dir)
-        outputs['mean_S0samples'] = self._gen_fname('mean_S0samples',
-                                                    cwd=out_dir)
+        for k in single_out:
+            outputs[k] = self._gen_fname(k, cwd=out_dir)
 
         if isdefined(self.inputs.rician) and self.inputs.rician:
             outputs['mean_tausamples'] = self._gen_fname('mean_tausamples',
                                                          cwd=out_dir)
 
+        for k in multi_out:
+            outputs[k] = []
+
         for i in xrange(1, self.inputs.n_fibres + 1):
-            outputs['merged_thsamples'].append(self._gen_fname(('merged_th%d'
-                                               'samples') % i,
-                                               cwd=out_dir))
-            outputs['merged_phsamples'].append(self._gen_fname(('merged_ph%d'
-                                               'samples') % i,
-                                               cwd=out_dir))
-            outputs['merged_fsamples'].append(self._gen_fname(('merged_f%d'
-                                              'samples') % i,
-                                              cwd=out_dir))
-
-            outputs['mean_thsamples'].append(self._gen_fname(('mean_th%d'
-                                             'samples') % i,
-                                             cwd=out_dir))
-            outputs['mean_phsamples'].append(self._gen_fname(('mean_ph%d'
-                                             'samples') % i,
-                                             cwd=out_dir))
-            outputs['mean_fsamples'].append(self._gen_fname(('mean_f%d'
-                                            'samples') % i,
-                                            cwd=out_dir))
             outputs['dyads'].append(self._gen_fname('dyads%d' % i,
                                     cwd=out_dir))
-            outputs['dyads_disp'].append(self._gen_fname(('dyads%d'
-                                         '_dispersion') % i, cwd=out_dir))
+            outputs['fsamples'].append(self._gen_fname('f%dsamples' % i,
+                                       cwd=out_dir))
+            outputs['mean_fsamples'].append(self._gen_fname(('mean_f%d'
+                                            'samples') % i, cwd=out_dir))
+            outputs['phsamples'].append(self._gen_fname('ph%dsamples' % i,
+                                        cwd=out_dir))
+            outputs['thsamples'].append(self._gen_fname('th%dsamples' % i,
+                                        cwd=out_dir))
+
         return outputs
 
 

nipype/interfaces/fsl/dti.py

@@ -107,6 +107,23 @@ class BEDPOSTXInputSpec(FSLXCommandInputSpec):
                              'nonlinearities, default off'))
 
 
+class BEDPOSTXOutputSpec(FSLXCommandOutputSpec):
+    mean_thsamples = OutputMultiPath(File(), desc=('Mean of '
+                                     'distribution on theta'))
+    mean_phsamples = OutputMultiPath(File(), desc=('Mean of '
+                                     'distribution on phi'))
+
+    merged_thsamples = OutputMultiPath(File(), desc=('Samples from '
+                                       'the distribution on theta'))
+    merged_phsamples = OutputMultiPath(File(), desc=('Samples from '
+                                       'the distribution on phi'))
+    merged_fsamples = OutputMultiPath(File(),
+                                      desc=('Samples from the distribution on '
+                                            'anisotropic volume fraction.'))
+    dyads_disp = OutputMultiPath(File(), desc=('Uncertainty on the '
+                                 ' estimated fiber orientation'))
+
+
 class BEDPOSTX(FSLXCommand):
     """
     BEDPOSTX stands for Bayesian Estimation of Diffusion Parameters Obtained
@@ -130,13 +147,13 @@ class BEDPOSTX(FSLXCommand):
     ...                     mask='mask.nii', n_fibres=1)
     >>> bedp.cmdline
     'bedpostx . --bvals=bvals --bvecs=bvecs --data=diffusion.nii \
---forcedir --logdir=logdir --mask=mask.nii --nfibres=1'
+--forcedir --logdir=. --mask=mask.nii --nfibres=1'
 
     """
 
     _cmd = 'bedpostx'
     input_spec = BEDPOSTXInputSpec
-    output_spec = FSLXCommandOutputSpec
+    output_spec = BEDPOSTXOutputSpec
     _can_resume = True
 
     def _run_interface(self, runtime):
@@ -164,6 +181,41 @@ def _run_interface(self, runtime):
 
         return super(BEDPOSTX, self)._run_interface(runtime)
 
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        out_dir = self._out_dir
+
+        post_fields = ['merged_thsamples', 'merged_phsamples',
+                       'merged_fsamples', 'dyads_disp']
+        for k in post_fields:
+            outputs[k] = []
+
+        for i in xrange(1, self.inputs.n_fibres + 1):
+            outputs['merged_thsamples'].append(self._gen_fname(('merged_th%d'
+                                               'samples') % i,
+                                               cwd=out_dir))
+            outputs['merged_phsamples'].append(self._gen_fname(('merged_ph%d'
+                                               'samples') % i,
+                                               cwd=out_dir))
+            outputs['merged_fsamples'].append(self._gen_fname(('merged_f%d'
+                                              'samples') % i,
+                                              cwd=out_dir))
+            outputs['mean_thsamples'].append(self._gen_fname(('mean_th%d'
+                                             'samples') % i,
+                                             cwd=out_dir))
+            outputs['mean_phsamples'].append(self._gen_fname(('mean_ph%d'
+                                             'samples') % i,
+                                             cwd=out_dir))
+            outputs['dyads_disp'].append(self._gen_fname(('dyads%d'
+                                         '_dispersion') % i, cwd=out_dir))
+
+        super_out = super(BEDPOSTX, self)._list_outputs()
+
+        for k, v in super_out.iteritems():
+            outputs[k] = v
+
+        return outputs
+
 
 class XFibresInputSpec(FSLXCommandInputSpec):
     gradnonlin = File(exists=True, argstr='--gradnonlin=%s',

nipype/interfaces/fsl/tests/test_auto_BEDPOSTX.py

@@ -89,10 +89,13 @@ def test_BEDPOSTX_inputs():
             yield assert_equal, getattr(inputs.traits()[key], metakey), value
 
 def test_BEDPOSTX_outputs():
-    output_map = dict(dyads=dict(),
+    output_map = dict(d_stdsamples=dict(),
+    dsamples=dict(),
+    dyads=dict(),
     dyads_disp=dict(),
     fsamples=dict(),
     mean_S0samples=dict(),
+    mean_d_stdsamples=dict(),
     mean_dsamples=dict(),
     mean_fsamples=dict(),
     mean_phsamples=dict(),
@@ -101,6 +104,8 @@ def test_BEDPOSTX_outputs():
     merged_fsamples=dict(),
     merged_phsamples=dict(),
     merged_thsamples=dict(),
+    phsamples=dict(),
+    thsamples=dict(),
     )
     outputs = BEDPOSTX.output_spec()
 

nipype/interfaces/fsl/tests/test_auto_FSLXCommand.py

@@ -82,18 +82,17 @@ def test_FSLXCommand_inputs():
             yield assert_equal, getattr(inputs.traits()[key], metakey), value
 
 def test_FSLXCommand_outputs():
-    output_map = dict(dyads=dict(),
-    dyads_disp=dict(),
+    output_map = dict(d_stdsamples=dict(),
+    dsamples=dict(),
+    dyads=dict(),
     fsamples=dict(),
     mean_S0samples=dict(),
+    mean_d_stdsamples=dict(),
     mean_dsamples=dict(),
     mean_fsamples=dict(),
-    mean_phsamples=dict(),
     mean_tausamples=dict(),
-    mean_thsamples=dict(),
-    merged_fsamples=dict(),
-    merged_phsamples=dict(),
-    merged_thsamples=dict(),
+    phsamples=dict(),
+    thsamples=dict(),
     )
     outputs = FSLXCommand.output_spec()
 

nipype/interfaces/fsl/tests/test_auto_XFibres.py

@@ -84,18 +84,17 @@ def test_XFibres_inputs():
             yield assert_equal, getattr(inputs.traits()[key], metakey), value
 
 def test_XFibres_outputs():
-    output_map = dict(dyads=dict(),
-    dyads_disp=dict(),
+    output_map = dict(d_stdsamples=dict(),
+    dsamples=dict(),
+    dyads=dict(),
     fsamples=dict(),
     mean_S0samples=dict(),
+    mean_d_stdsamples=dict(),
     mean_dsamples=dict(),
     mean_fsamples=dict(),
-    mean_phsamples=dict(),
     mean_tausamples=dict(),
-    mean_thsamples=dict(),
-    merged_fsamples=dict(),
-    merged_phsamples=dict(),
-    merged_thsamples=dict(),
+    phsamples=dict(),
+    thsamples=dict(),
     )
     outputs = XFibres.output_spec()
 

nipype/workflows/dmri/fsl/dti.py

@@ -59,22 +59,104 @@ def create_bedpostx_pipeline(name='bedpostx', params={}):
     xfibres = pe.MapNode(xfib_if, name='xfibres',
                          iterfield=['dwi', 'mask'])
 
-    out_fields = xfib_if.output_spec().get().keys()
+    make_dyads = pe.MapNode(fsl.MakeDyadicVectors(), name="make_dyads",
+                            iterfield=['theta_vol', 'phi_vol'])
+    out_fields = ['dyads', 'dyads_disp',
+                  'thsamples', 'phsamples', 'fsamples',
+                  'mean_thsamples', 'mean_phsamples', 'mean_fsamples']
 
     outputnode = pe.Node(niu.IdentityInterface(fields=out_fields),
                          name='outputnode')
 
     wf = pe.Workflow(name=name)
     wf.connect([
-        (inputnode, slice_dwi, [('dwi', 'in_file')]),
-        (inputnode, slice_msk, [('mask', 'in_file')]),
-        (slice_dwi, mask_dwi,  [('out_files', 'in_file')]),
-        (slice_msk, mask_dwi,  [('out_files', 'in_file2')]),
-        (slice_dwi, xfibres,   [('out_files', 'dwi')]),
-        (mask_dwi, xfibres,    [('out_file', 'mask')]),
-        (inputnode, xfibres,   [('bvecs', 'bvecs'),
-                                ('bvals', 'bvals')]),
-        (xfibres, outputnode,  [((f, transpose), f) for f in out_fields])
+        (inputnode, slice_dwi,  [('dwi', 'in_file')]),
+        (inputnode, slice_msk,  [('mask', 'in_file')]),
+        (slice_dwi, mask_dwi,   [('out_files', 'in_file')]),
+        (slice_msk, mask_dwi,   [('out_files', 'in_file2')]),
+        (slice_dwi, xfibres,    [('out_files', 'dwi')]),
+        (mask_dwi, xfibres,     [('out_file', 'mask')]),
+        (inputnode, xfibres,    [('bvecs', 'bvecs'),
+                                 ('bvals', 'bvals')]),
+        (inputnode, make_dyads, [('mask', 'mask')])
     ])
 
+    mms = {}
+    for k in ['thsamples', 'phsamples', 'fsamples']:
+        mms[k] = merge_and_mean(k)
+        wf.connect([
+            (xfibres, mms[k], [(k, 'inputnode.in_files')]),
+            (mms[k], outputnode, [('outputnode.merged', k),
+                                  ('outputnode.mean', 'mean_%s' % k)])
+
+        ])
+
+    # m_mdsamples = pe.Node(fsl.Merge(dimension="z"),
+    #                       name="merge_mean_dsamples")
+    wf.connect([
+        (mms['thsamples'], make_dyads, [('outputnode.merged', 'theta_vol')]),
+        (mms['phsamples'], make_dyads, [('outputnode.merged', 'phi_vol')]),
+        #(xfibres, m_mdsamples,  [('mean_dsamples', 'in_files')]),
+        (make_dyads, outputnode, [('dyads', 'dyads'),
+                                  ('dispersion', 'dyads_disp')])
+    ])
     return wf
+
+
+def merge_and_mean(name='mm'):
+    inputnode = pe.Node(niu.IdentityInterface(fields=['in_files']),
+                        name='inputnode')
+    outputnode = pe.Node(niu.IdentityInterface(fields=['merged', 'mean']),
+                         name='outputnode')
+    merge = pe.MapNode(fsl.Merge(dimension='z'), name='Merge',
+                       iterfield=['in_files'])
+    mean = pe.MapNode(fsl.ImageMaths(op_string='-Tmean'), name='Mean',
+                      iterfield=['in_file'])
+
+    wf = pe.Workflow(name=name)
+    wf.connect([
+        (inputnode, merge,  [(('in_files', transpose), 'in_files')]),
+        (merge, mean,       [('merged_file', 'in_file')]),
+        (merge, outputnode, [('merged_file', 'merged')]),
+        (mean, outputnode,  [('out_file', 'mean')])
+    ])
+    return wf
+
+
+def gen_chunks(dwi, mask, nchunks=6):
+    import nibabel as nb
+    import numpy as np
+
+    mask = nb.load(mask).get_data()
+    mask[mask > 0] = 1
+    mask[mask < 1] = 0
+    nzels = np.nonzero(mask)
+    np.savetxt('nonzeroidx.txt', nzels)
+
+    dshape = mask.shape
+    mask = mask.reshape(-1).astype(np.uint8)
+    els = np.sum(mask)
+    chunkels = round(els / nchunks)
+
+    data = mask.copy()
+
+    out_files = []
+    out_masks = []
+
+    for i in xrange(chunkels):
+        first = i * chunkels
+        last = (i+1) * chunkels + 1
+
+        chunk = data[first:last, ...]
+
+        fname = 'chunk%05d.nii.gz' % i
+        nb.Nifti1Image(chunk, None, None).to_filename(fname)
+        out_files.append(fname)
+
+        mname = 'mask%05d.nii.gz' % i
+        nb.Nifti1Image(mask[first:last],
+                       None, None).to_filename(mname)
+        out_masks.append(mname)
+
+    return out_files, out_masks
+