Updates and a few fixes

Author: erikz

examples/connectivity_tutorial_advanced.py

@@ -53,14 +53,8 @@
 import nibabel as nb
 import os, os.path as op                     # system functions
 import cmp                                   # connectome mapper
-
-"""
-We import several voxel, data, and affine grabbing functions from the Camino workflows, as well as 
-some functions to select the correct files (in our case, aparc+aseg.mgz) from the list outputted by the FreeSurferSource interfaces.
-"""
-
-from nipype.workflows.camino.connectivity_mapping import (get_vox_dims, get_data_dims,
- get_affine, select_aparc, select_aparc_annot, get_first_image)
+from nipype.workflows.camino.connectivity_mapping import select_aparc_annot, get_first_image
+from nipype.workflows.mrtrix.diffusion import get_vox_dims_as_tuple
 
 """
 This needs to point to the freesurfer subjects directory (Recon-all must have been run on subj1 from the FSL course data)
@@ -177,6 +171,14 @@
 
 fsl2mrtrix = pe.Node(interface=mrtrix.FSL2MRTrix(),name='fsl2mrtrix')
 
+"""
+Distortions induced by eddy currents are corrected prior to fitting the tensors.
+The first image is used as a reference for which to warp the others.
+"""
+
+eddycorrect = pe.Node(interface=fsl.EddyCorrect(),name='eddycorrect')
+eddycorrect.inputs.ref_num = 1
+
 """
 Tensors are fitted to each voxel in the diffusion-weighted image and from these three maps are created:
 	* Major eigenvector in each voxel
@@ -252,6 +254,8 @@
 tracks2prob = pe.Node(interface=mrtrix.Tracks2Prob(),name='tracks2prob')
 tracks2prob.inputs.colour = True
 tck2trk = pe.Node(interface=mrtrix.MRTrix2TrackVis(),name='tck2trk')
+tck2trk.inputs.flipy = True
+tck2trk.inputs.flipz = True
 
 """
 Structural segmentation nodes
@@ -283,6 +287,7 @@
 inverse.inputs.interp = ('nearestneighbour')
 inverse.inputs.apply_xfm = True
 inverse_AparcAseg = inverse.clone('inverse_AparcAseg')
+inverseROIsToB0 = inverse.clone('inverseROIsToB0')
 
 """
 Parcellation is performed given the aparc+aseg image from Freesurfer.
@@ -403,7 +408,8 @@
 
 mapping.connect([(inputnode, fsl2mrtrix, [("bvecs", "bvec_file"),
 												("bvals", "bval_file")])])
-mapping.connect([(inputnode, dwi2tensor,[("dwi","in_file")])])
+mapping.connect([(inputnode, eddycorrect,[("dwi","in_file")])])
+mapping.connect([(eddycorrect, dwi2tensor,[("eddy_corrected","in_file")])])
 mapping.connect([(fsl2mrtrix, dwi2tensor,[("encoding_file","encoding_file")])])
 
 mapping.connect([(dwi2tensor, tensor2vector,[['tensor','in_file']]),
@@ -417,7 +423,7 @@
 fractional anisotropy image, and thresholds it to get the single-fiber voxels.
 """
 
-mapping.connect([(inputnode, MRconvert,[("dwi","in_file")])])
+mapping.connect([(eddycorrect, MRconvert,[("eddy_corrected","in_file")])])
 mapping.connect([(MRconvert, threshold_b0,[("converted","in_file")])])
 mapping.connect([(threshold_b0, median3d,[("out_file","in_file")])])
 mapping.connect([(median3d, erode_mask_firstpass,[("out_file","in_file")])])
@@ -430,8 +436,8 @@
 Here the thresholded white matter mask is created for seeding the tractography.
 """
 
-mapping.connect([(inputnode, bet,[("dwi","in_file")])])
-mapping.connect([(inputnode, gen_WM_mask,[("dwi","in_file")])])
+mapping.connect([(eddycorrect, bet,[("eddy_corrected","in_file")])])
+mapping.connect([(eddycorrect, gen_WM_mask,[("eddy_corrected","in_file")])])
 mapping.connect([(bet, gen_WM_mask,[("mask_file","binary_mask")])])
 mapping.connect([(fsl2mrtrix, gen_WM_mask,[("encoding_file","encoding_file")])])
 mapping.connect([(gen_WM_mask, threshold_wmmask,[("WMprobabilitymap","in_file")])])
@@ -440,15 +446,15 @@
 Next we estimate the fiber response distribution.
 """
 
-mapping.connect([(inputnode, estimateresponse,[("dwi","in_file")])])
+mapping.connect([(eddycorrect, estimateresponse,[("eddy_corrected","in_file")])])
 mapping.connect([(fsl2mrtrix, estimateresponse,[("encoding_file","encoding_file")])])
 mapping.connect([(threshold_FA, estimateresponse,[("out_file","mask_image")])])
 
 """
 Run constrained spherical deconvolution.
 """
 
-mapping.connect([(inputnode, csdeconv,[("dwi","in_file")])])
+mapping.connect([(eddycorrect, csdeconv,[("eddy_corrected","in_file")])])
 mapping.connect([(gen_WM_mask, csdeconv,[("WMprobabilitymap","mask_image")])])
 mapping.connect([(estimateresponse, csdeconv,[("response","response_file")])])
 mapping.connect([(fsl2mrtrix, csdeconv,[("encoding_file","encoding_file")])])
@@ -460,42 +466,45 @@
 mapping.connect([(threshold_wmmask, probCSDstreamtrack,[("out_file","seed_file")])])
 mapping.connect([(csdeconv, probCSDstreamtrack,[("spherical_harmonics_image","in_file")])])
 mapping.connect([(probCSDstreamtrack, tracks2prob,[("tracked","in_file")])])
-mapping.connect([(inputnode, tracks2prob,[("dwi","template_file")])])
+mapping.connect([(eddycorrect, tracks2prob,[("eddy_corrected","template_file")])])
 
 """
 Structural Processing
 ---------------------
 First, we coregister the structural image to the diffusion image and then obtain the inverse of transformation.
 """
 
-mapping.connect([(inputnode, coregister,[('dwi','in_file')])])
+mapping.connect([(eddycorrect, coregister,[("eddy_corrected","in_file")])])
 mapping.connect([(mri_convert_Brain, coregister,[('out_file','reference')])])
 mapping.connect([(coregister, convertxfm,[('out_matrix_file','in_file')])])
-mapping.connect([(inputnode, inverse,[('dwi','reference')])])
+mapping.connect([(eddycorrect, inverse,[("eddy_corrected","reference")])])
 mapping.connect([(convertxfm, inverse,[('out_file','in_matrix_file')])])
 mapping.connect([(mri_convert_Brain, inverse,[('out_file','in_file')])])
 
 """
 The b0 image is upsampled to the same dimensions as the parcellated structural image to improve their coregistration.
 """
 
-mapping.connect([(inputnode, resampleb0,[(('dwi', get_first_image), 'in_file')])])
+mapping.connect([(eddycorrect, resampleb0,[(('eddy_corrected', get_first_image), 'in_file')])])
 mapping.connect([(resampleb0, inverse_AparcAseg,[('out_file','reference')])])
 mapping.connect([(convertxfm, inverse_AparcAseg,[('out_file','in_matrix_file')])])
+mapping.connect([(eddycorrect, inverseROIsToB0,[(('eddy_corrected', get_first_image), 'reference')])])
+mapping.connect([(convertxfm, inverseROIsToB0,[('out_file','in_matrix_file')])])
 
 """
 The parcellation is connected for transformation into diffusion space.
 """
 
 mapping.connect([(inputnode, parcellate,[("subject_id","subject_id")])])
 mapping.connect([(parcellate, inverse_AparcAseg,[('roi_file','in_file')])])
+mapping.connect([(parcellate, inverseROIsToB0,[('roi_file','in_file')])])
 
 """
 The MRtrix-tracked fibers are converted to TrackVis format (with voxel and data dimensions grabbed from the DWI).
 The connectivity matrix is created with the .trk fibers and the coregistered parcellation file.
 """
 
-mapping.connect([(inputnode, tck2trk,[("dwi","image_file")])])
+mapping.connect([(eddycorrect, tck2trk,[("eddy_corrected","image_file")])])
 mapping.connect([(probCSDstreamtrack, tck2trk,[("tracked","in_file")])])
 mapping.connect([(tck2trk, creatematrix,[("out_file","tract_file")])])
 mapping.connect([(inputnode, creatematrix,[("subject_id","out_matrix_file")])])
@@ -519,7 +528,8 @@
 mapping.connect([(mris_convertLHlabels, giftiLabels,[("converted","in1")])])
 mapping.connect([(mris_convertRHlabels, giftiLabels,[("converted","in2")])])
 
-mapping.connect([(inputnode, niftiVolumes,[("dwi","in2")])])
+mapping.connect([(eddycorrect, niftiVolumes,[("eddy_corrected","in2")])])
+
 mapping.connect([(mri_convert_Brain, niftiVolumes,[("out_file","in3")])])
 mapping.connect([(inverse_AparcAseg, niftiVolumes,[("out_file","in1")])])
 

nipype/interfaces/cmtk/cmtk.py

@@ -175,7 +175,7 @@ def cmat(track_file, roi_file, resolution_network_file, matrix_name, matrix_mat_
 		G.add_node(int(u), d)
 		# compute a position for the node based on the mean position of the
 		# ROI in voxel coordinates (segmentation volume )
-		G.node[int(u)]['dn_position'] = str(tuple(np.mean( np.where(roiData== int(d["dn_correspondence_id"]) ) , axis = 1)))
+		G.node[int(u)]['dn_position'] = tuple(np.mean( np.where(roiData== int(d["dn_correspondence_id"]) ) , axis = 1))
 
     dis = 0
     for i in range(endpoints.shape[0]):
@@ -275,6 +275,10 @@ def cmat(track_file, roi_file, resolution_network_file, matrix_name, matrix_mat_
     fibdev_dict['cmatrix'] = fibdev_mlab
 
     path, name, ext = split_filename(matrix_mat_name)
+    if not ext == '.mat':
+        ext = '.mat'
+        matrix_mat_name = matrix_mat_name + ext
+
     print 'Writing matlab matrix as {mat}'.format(mat=matrix_mat_name)
     sio.savemat(matrix_mat_name, numfib_dict)
 
@@ -306,16 +310,12 @@ def save_fibers(oldhdr, oldfib, fname, indices):
     """ Stores a new trackvis file fname using only given indices """
     import nibabel.trackvis as tv
     import cPickle
-    
     hdrnew = oldhdr.copy()
-
     outstreams = []
     for i in indices:
         outstreams.append( oldfib[i] )
-
     n_fib_out = len(outstreams)
     hdrnew['n_count'] = n_fib_out
-
     print "Writing final no orphan fibers as %s" % fname 
     tv.write(fname,outstreams,hdrnew)
 
@@ -325,10 +325,10 @@ class CreateMatrixInputSpec(TraitedSpec):
     tract_file = File(exists=True, mandatory=True, desc='Trackvis tract file')
     resolution_network_file = File(exists=True, mandatory=True, desc='Parcellation files from Connectome Mapping Toolkit')
     out_matrix_file = File(genfile = True, desc='NetworkX graph describing the connectivity')
-    out_matrix_mat_file = File(genfile = True, desc='Matlab matrix describing the connectivity')
-    out_mean_fiber_length_matrix_mat_file = File(genfile = True, desc='Matlab matrix describing the mean fiber lengths between each node.')
-    out_fiber_length_std_matrix_mat_file = File(genfile = True, desc='Matlab matrix describing the deviation in fiber lengths connecting each node.')
-    out_endpoint_array_name = File(genfile = True, desc='Name for the generated endpoint arrays')
+    out_matrix_mat_file = File('cmatrix.mat', usedefault=True, desc='Matlab matrix describing the connectivity')
+    out_mean_fiber_length_matrix_mat_file = File(genfile=True, desc='Matlab matrix describing the mean fiber lengths between each node.')
+    out_fiber_length_std_matrix_mat_file = File(genfile=True, desc='Matlab matrix describing the deviation in fiber lengths connecting each node.')
+    out_endpoint_array_name = File(genfile=True, desc='Name for the generated endpoint arrays')
 
 class CreateMatrixOutputSpec(TraitedSpec):
     matrix_file = File(desc='NetworkX graph describing the connectivity')
@@ -349,8 +349,8 @@ class CreateMatrix(BaseInterface):
     Example
     -------
 
-    >>> import nipype.interfaces.cmtk.cmtk as ck
-    >>> conmap = ck.CreateMatrix()
+    >>> import nipype.interfaces.cmtk as cmtk
+    >>> conmap = cmtk.CreateMatrix()
     >>> conmap.roi_file = 'fsLUT_aparc+aseg.nii'
     >>> conmap.dict_file = 'fsLUT_aparc+aseg.pck'
     >>> conmap.tract_file = 'fibers.trk'
@@ -365,17 +365,29 @@ def _run_interface(self, runtime):
             path, name, _ = split_filename(self.inputs.out_matrix_file)
             matrix_file = op.abspath(name + '.pck')
         else:
-            matrix_file = self._gen_outfilename('pck')
-        if isdefined(self.inputs.out_matrix_mat_file):
-            path, name, _ = split_filename(self.inputs.out_matrix_mat_file)
-            matrix_mat_file = op.abspath(name + '.mat')
+            matrix_file = self._gen_outfilename('.pck')
+
+        matrix_mat_file = op.abspath(self.inputs.out_matrix_mat_file)
+        path, name, ext = split_filename(matrix_mat_file)
+        if not ext == '.mat':
+            ext = '.mat'
+            matrix_mat_file = matrix_mat_file + ext
+            
+        if isdefined(self.inputs.out_mean_fiber_length_matrix_mat_file):
+            mean_fiber_length_matrix_mat_file = op.abspath(self.inputs.out_mean_fiber_length_matrix_mat_file)
         else:
-            matrix_mat_file = self._gen_outfilename('mat')
-
+            mean_fiber_length_matrix_name = op.abspath(self._gen_outfilename('_mean_fiber_length'))
+        
+        if isdefined(self.inputs.out_fiber_length_std_matrix_mat_file):
+            fiber_length_std_matrix_mat_file = op.abspath(self.inputs.out_fiber_length_std_matrix_mat_file)
+        else:
+            fiber_length_std_matrix_name = op.abspath(self._gen_outfilename('_fiber_length_std.mat'))
+        
         if not isdefined(self.inputs.out_endpoint_array_name):
             _, endpoint_name , _ = split_filename(self.inputs.tract_file)
+            endpoint_name = op.abspath(endpoint_name)
         else:
-            endpoint_name = self.inputs.out_endpoint_array_name
+            endpoint_name = op.abspath(self.inputs.out_endpoint_array_name)
 
         if isdefined(self.inputs.dict_file):
             cmat(self.inputs.tract_file, self.inputs.roi_file, self.inputs.resolution_network_file,
@@ -390,26 +402,27 @@ def _list_outputs(self):
         outputs = self.output_spec().get()
         if isdefined(self.inputs.out_matrix_file):
             path, name, _ = split_filename(self.inputs.out_matrix_file)
-            outputs['matrix_file']= op.abspath(name + '.pck')
+            outputs['matrix_file'] = op.abspath(name + '.pck')
         else:
-            outputs['matrix_file']=op.abspath(self._gen_outfilename('pck'))
-        if isdefined(self.inputs.out_matrix_mat_file):
-            path, name, _ = split_filename(self.inputs.out_matrix_mat_file)
-            outputs['matrix_mat_file']= op.abspath(name + '.mat')
-            outputs['mean_fiber_length_matrix_mat_file']= op.abspath(name + '_mean_fiber_length.mat')
-            outputs['fiber_length_std_matrix_mat_file']= op.abspath(name + '_fiber_length_std.mat')
+            outputs['matrix_file'] = op.abspath(self._gen_outfilename('.pck'))
+
+        matrix_mat_file = op.abspath(self.inputs.out_matrix_mat_file)
+        path, name, ext = split_filename(matrix_mat_file)
+        if not ext == '.mat':
+            ext = '.mat'
+            matrix_mat_file = matrix_mat_file + ext
+
+        outputs['matrix_mat_file'] = matrix_mat_file        
+        if isdefined(self.inputs.out_mean_fiber_length_matrix_mat_file):
+            outputs['mean_fiber_length_matrix_mat_file']= op.abspath(self.inputs.out_mean_fiber_length_matrix_mat_file)
         else:
-            name = op.abspath(self._gen_outfilename('mat'))
-            outputs['matrix_mat_file'] = name
-            outputs['mean_fiber_length_matrix_mat_file'] = op.abspath(name + '_mean_fiber_length.mat')
-            outputs['fiber_length_std_matrix_mat_file'] = op.abspath(name + '_fiber_length_std.mat')
-
-        if isdefined(self.inputs.out_matrix_mat_file):
-            path, name, _ = split_filename(self.inputs.out_matrix_mat_file)
-            outputs['matrix_mat_file']= op.abspath(name + '.mat')
+            outputs['mean_fiber_length_matrix_mat_file']=op.abspath(self._gen_outfilename('_mean_fiber_length.mat'))
+        
+        if isdefined(self.inputs.out_fiber_length_std_matrix_mat_file):
+            outputs['fiber_length_std_matrix_mat_file']= op.abspath(self.inputs.out_fiber_length_std_matrix_mat_file)
         else:
-            outputs['matrix_mat_file']=op.abspath(self._gen_outfilename('mat'))
-
+            outputs['fiber_length_std_matrix_mat_file']= op.abspath(self._gen_outfilename('_fiber_length_std.mat'))
+        
         if isdefined(self.inputs.out_endpoint_array_name):
             outputs['endpoint_file'] = op.abspath(self.inputs.out_endpoint_array_name + '_endpoints.npy')
             outputs['endpoint_file_mm'] = op.abspath(self.inputs.out_endpoint_array_name + '_endpointsmm.npy')
@@ -429,8 +442,11 @@ def _list_outputs(self):
         return outputs
 
     def _gen_outfilename(self, ext):
-        _, name , _ = split_filename(self.inputs.tract_file)
-        return name + '.' + ext
+		if isdefined(self.inputs.out_matrix_file):
+			_, name , _ = split_filename(self.inputs.out_matrix_file)
+		else:
+			_, name , _ = split_filename(self.inputs.tract_file)
+		return name + ext
 
 class ROIGenInputSpec(BaseInterfaceInputSpec):
     aparc_aseg_file = File(exists=True, mandatory=True, desc='Freesurfer aparc+aseg file')
@@ -451,8 +467,8 @@ class ROIGen(BaseInterface):
     Example
     -------
 
-    >>> import nipype.interfaces.cmtk.cmtk as ck
-    >>> rg = ck.ROIGen()
+    >>> import nipype.interfaces.cmtk as cmtk
+    >>> rg = cmtk.ROIGen()
     >>> rg.inputs.aparc_aseg_file = 'aparc+aseg.nii'
     >>> rg.inputs.use_freesurfer_LUT = True
     >>> rg.inputs.freesurfer_dir = '/usr/local/freesurfer'