Camino group tutorial runs properly

erikz · Erik Ziegler · commit 7b80b0054da9 · 2012-03-06T00:43:07.000+01:00
diff --git a/examples/dmri_group_connectivity_camino.py b/examples/dmri_group_connectivity_camino.py
@@ -14,7 +14,7 @@
 
 We perform this analysis using one healthy subject and two subjects who suffer from Parkinson's disease.
 
-The whole package (754 mb as .tar.gz / 1 gb uncompressed) including the Freesurfer directories for these subjects, can be acquired from here:
+The whole package (960 mb as .tar.gz / 1.3 gb uncompressed) including the Freesurfer directories for these subjects, can be acquired from here:
 
     * http://db.tt/b6F1t0QV
 
@@ -79,7 +79,7 @@
 
 group_list = {}
 group_list['controls'] = ['cont17']
-group_list['parkinsons'] = ['pat07', 'pat20']
+group_list['parkinsons'] = ['pat10', 'pat20']
 
 """
 The output directory must be named as well.
@@ -131,10 +131,10 @@
     Define the parcellation scheme to use.
     """
 
-    parcellation_name = 'scale500'
+    parcellation_scheme = 'NativeFreesurfer'
     cmp_config = cmp.configuration.PipelineConfiguration()
-    cmp_config.parcellation_scheme = "Lausanne2008"
-    l1pipeline.inputs.connectivity.inputnode.resolution_network_file = cmp_config._get_lausanne_parcellation('Lausanne2008')[parcellation_name]['node_information_graphml']
+    cmp_config.parcellation_scheme = parcellation_scheme
+    l1pipeline.inputs.connectivity.inputnode.resolution_network_file = cmp_config._get_lausanne_parcellation(parcellation_scheme)['freesurferaparc']['node_information_graphml']
 
     """
     The first level pipeline we have tweaked here is run within the for loop.
diff --git a/examples/dmri_group_connectivity_mrtrix.py b/examples/dmri_group_connectivity_mrtrix.py
@@ -14,7 +14,7 @@
 
 We perform this analysis using one healthy subject and two subjects who suffer from Parkinson's disease.
 
-The whole package (754 mb as .tar.gz / 1 gb uncompressed) including the Freesurfer directories for these subjects, can be acquired from here:
+The whole package (960 mb as .tar.gz / 1.3 gb uncompressed) including the Freesurfer directories for these subjects, can be acquired from here:
 
     * http://db.tt/b6F1t0QV
 
@@ -77,7 +77,7 @@
 
 group_list = {}
 group_list['controls'] = ['cont17']
-group_list['parkinsons'] = ['pat07', 'pat20']
+group_list['parkinsons'] = ['pat10', 'pat20']
 
 """
 The output directory must be named as well.
diff --git a/nipype/interfaces/cmtk/cmtk.py b/nipype/interfaces/cmtk/cmtk.py
@@ -212,7 +212,7 @@ def cmat(track_file, roi_file, resolution_network_file, matrix_name, matrix_mat_
             # compute a position for the node based on the mean position of the
             # ROI in voxel coordinates (segmentation volume )
             xyz = tuple(np.mean(np.where(np.flipud(roiData) == int(d["dn_correspondence_id"])) , axis=1))
-            G.node[int(u)]['dn_position'] = tuple([xyz[0], xyz[1], xyz[2]])
+            G.node[int(u)]['dn_position'] = tuple([xyz[0], xyz[2], -xyz[1]])
 
     if intersections:
         iflogger.info("Filtering tractography from intersections")
diff --git a/nipype/interfaces/cmtk/parcellation.py b/nipype/interfaces/cmtk/parcellation.py
@@ -18,6 +18,10 @@
 import shutil
 from nipype.utils.misc import package_check
 import warnings
+import logging
+
+logging.basicConfig()
+iflogger = logging.getLogger('interface')
 
 have_cmp = True
 try:
@@ -30,8 +34,8 @@
     from cmp.util import runCmd
 
 def create_annot_label(subject_id, subjects_dir, fs_dir, parcellation_name):
-    print "Create the cortical labels necessary for our ROIs"
-    print "================================================="
+    iflogger.info("Create the cortical labels necessary for our ROIs")
+    iflogger.info("=================================================")
     fs_label_dir = op.join(op.join(subjects_dir,subject_id), 'label')
     output_dir = op.abspath(op.curdir)
     paths = []
@@ -68,15 +72,15 @@ def create_annot_label(subject_id, subjects_dir, fs_dir, parcellation_name):
         mris_cmd = 'mris_ca_label %s %s "%s/surf/%s.sphere.reg" "%s" "%s" ' % (subject_id, out[0],
 		    op.join(subjects_dir,subject_id), out[0], cmp_config.get_lausanne_atlas(out[1]), op.join(fs_label_dir, out[2]))
         runCmd( mris_cmd, log )
-        print '-----------'
+        iflogger.info('-----------')
 
         annot = '--annotation "%s"' % out[4]
 
         mri_an_cmd = 'mri_annotation2label --subject %s --hemi %s --outdir "%s" %s' % (subject_id, out[0], op.join(output_dir, out[3]), annot)
-        print mri_an_cmd
+        iflogger.info(mri_an_cmd)
         runCmd( mri_an_cmd, log )
-        print '-----------'
-        print os.environ['SUBJECTS_DIR']
+        iflogger.info('-----------')
+        iflogger.info(os.environ['SUBJECTS_DIR'])
     # extract cc and unknown to add to tractography mask, we do not want this as a region of interest
     # in FS 5.0, unknown and corpuscallosum are not available for the 35 scale (why?),
     # but for the other scales only, take the ones from _60
@@ -97,12 +101,12 @@ def create_annot_label(subject_id, subjects_dir, fs_dir, parcellation_name):
     mri_cmd = 'mri_convert -i "%s/mri/aseg.mgz" -o "%s/mri/aseg.nii.gz"' % (op.join(subjects_dir,subject_id), op.join(subjects_dir,subject_id))
     runCmd( mri_cmd, log )
 
-    print "[ DONE ]"
+    iflogger.info("[ DONE ]")
 
 def create_roi(subject_id, subjects_dir, fs_dir, parcellation_name):
     """ Creates the ROI_%s.nii.gz files using the given parcellation information
     from networks. Iteratively create volume. """
-    print "Create the ROIs:"
+    iflogger.info("Create the ROIs:")
     output_dir = op.abspath(op.curdir)
     fs_dir = op.join(subjects_dir, subject_id)
     cmp_config = cmp.configuration.PipelineConfiguration()
@@ -113,9 +117,9 @@ def create_roi(subject_id, subjects_dir, fs_dir, parcellation_name):
     aseg = nb.load(op.join(fs_dir, 'mri', 'aseg.nii.gz'))
     asegd = aseg.get_data()
 
-    print "Working on parcellation: "
-    print cmp_config._get_lausanne_parcellation('Lausanne2008')[parcellation_name]
-    print "========================"
+    iflogger.info("Working on parcellation: ")
+    iflogger.info(cmp_config._get_lausanne_parcellation('Lausanne2008')[parcellation_name])
+    iflogger.info("========================")
     pg = nx.read_graphml(pgpath)
     # each node represents a brain region
     # create a big 256^3 volume for storage of all ROIs
@@ -124,30 +128,30 @@ def create_roi(subject_id, subjects_dir, fs_dir, parcellation_name):
     count = 0
     for brk, brv in pg.nodes_iter(data=True):
         count = count + 1
-        print brv
-        print brk
+        iflogger.info(brv)
+        iflogger.info(brk)
         if brv['dn_hemisphere'] == 'left':
             hemi = 'lh'
         elif brv['dn_hemisphere'] == 'right':
             hemi = 'rh'
         if brv['dn_region'] == 'subcortical':
-            print brv
-            print "---------------------"
-            print "Work on brain region: %s" % (brv['dn_region'])
-            print "Freesurfer Name: %s" %  brv['dn_fsname']
-            print "Region %s of %s " % (count, pg.number_of_nodes())
-            print "---------------------"
+            iflogger.info(brv)
+            iflogger.info("---------------------")
+            iflogger.info("Work on brain region: %s" % (brv['dn_region']))
+            iflogger.info("Freesurfer Name: %s" %  brv['dn_fsname'])
+            iflogger.info("Region %s of %s " % (count, pg.number_of_nodes()))
+            iflogger.info("---------------------")
             # if it is subcortical, retrieve roi from aseg
             idx = np.where(asegd == int(brv['dn_fs_aseg_val']))
             rois[idx] = int(brv['dn_correspondence_id'])
 
         elif brv['dn_region'] == 'cortical':
-            print brv
-            print "---------------------"
-            print "Work on brain region: %s" % (brv['dn_region'])
-            print "Freesurfer Name: %s" %  brv['dn_fsname']
-            print "Region %s of %s " % (count, pg.number_of_nodes())
-            print "---------------------"
+            iflogger.info(brv)
+            iflogger.info("---------------------")
+            iflogger.info("Work on brain region: %s" % (brv['dn_region']))
+            iflogger.info("Freesurfer Name: %s" %  brv['dn_fsname'])
+            iflogger.info("Region %s of %s " % (count, pg.number_of_nodes()))
+            iflogger.info("---------------------")
 
             labelpath = op.join(output_dir, parval['fs_label_subdir_name'] % hemi)
             # construct .label file name
@@ -180,11 +184,11 @@ def create_roi(subject_id, subjects_dir, fs_dir, parcellation_name):
         img = nb.Nifti1Image(rois, aseg.get_affine(), hdr2)
         nb.save(img, out_roi)
 
-    print "[ DONE ]"
+    iflogger.info("[ DONE ]")
 
 
 def create_wm_mask(subject_id, subjects_dir, fs_dir, parcellation_name):
-    print "Create white matter mask"
+    iflogger.info("Create white matter mask")
     fs_dir = op.join(subjects_dir, subject_id)
     cmp_config = cmp.configuration.PipelineConfiguration()
     cmp_config.parcellation_scheme = "Lausanne2008"
@@ -291,32 +295,32 @@ def create_wm_mask(subject_id, subjects_dir, fs_dir, parcellation_name):
     # now remove all the structures from the white matter
     idx = np.where( (csfA != 0) | (csfB != 0) | (gr_ncl != 0) | (remaining != 0) )
     wmmask[idx] = 0
-    print "Removing lateral ventricles and eroded grey nuclei and brainstem from white matter mask"
+    iflogger.info("Removing lateral ventricles and eroded grey nuclei and brainstem from white matter mask")
 
     # ADD voxels from 'cc_unknown.nii.gz' dataset
     ccun = nb.load(op.join(fs_dir, 'label', 'cc_unknown.nii.gz'))
     ccund = ccun.get_data()
     idx = np.where(ccund != 0)
-    print "Add corpus callosum and unknown to wm mask"
+    iflogger.info("Add corpus callosum and unknown to wm mask")
     wmmask[idx] = 1
 
     # check if we should subtract the cortical rois from this parcellation
     parval = cmp_config._get_lausanne_parcellation('Lausanne2008')[parcellation_name]
-    print "Loading %s to subtract cortical ROIs from white matter mask" % ('ROI_%s.nii.gz' % parcellation_name)
+    iflogger.info("Loading %s to subtract cortical ROIs from white matter mask" % ('ROI_%s.nii.gz' % parcellation_name))
     roi = nb.load(op.join(op.curdir, 'ROI_%s.nii.gz' % parcellation_name))
     roid = roi.get_data()
     assert roid.shape[0] == wmmask.shape[0]
     pg = nx.read_graphml(pgpath)
     for brk, brv in pg.nodes_iter(data=True):
         if brv['dn_region'] == 'cortical':
-            print "Subtracting region %s with intensity value %s" % (brv['dn_region'], brv['dn_correspondence_id'])
+            iflogger.info("Subtracting region %s with intensity value %s" % (brv['dn_region'], brv['dn_correspondence_id']))
             idx = np.where(roid == int(brv['dn_correspondence_id']))
             wmmask[idx] = 0
 
     # output white matter mask. crop and move it afterwards
     wm_out = op.join(fs_dir, 'mri', 'fsmask_1mm.nii.gz')
     img = nb.Nifti1Image(wmmask, fsmask.get_affine(), fsmask.get_header() )
-    print "Save white matter mask: %s" % wm_out
+    iflogger.info("Save white matter mask: %s" % wm_out)
     nb.save(img, wm_out)
 
 def crop_and_move_datasets(subject_id, subjects_dir, fs_dir, parcellation_name, out_roi_file):
@@ -328,7 +332,7 @@ def crop_and_move_datasets(subject_id, subjects_dir, fs_dir, parcellation_name,
     reg_path = out_roi_file
     output_dir = op.abspath(op.curdir)
 
-    print "Cropping and moving datasets to %s" % output_dir
+    iflogger.info("Cropping and moving datasets to %s" % output_dir)
     ds = [
           (op.join(fs_dir, 'mri', 'aseg.nii.gz'), op.join(output_dir, 'aseg.nii.gz') ),
           (op.join(fs_dir, 'mri', 'ribbon.nii.gz'), op.join(output_dir, 'ribbon.nii.gz') ),
@@ -339,7 +343,7 @@ def crop_and_move_datasets(subject_id, subjects_dir, fs_dir, parcellation_name,
     ds.append( (op.join(op.curdir, 'ROI_%s.nii.gz' % parcellation_name), op.join(op.curdir, 'ROI_HR_th.nii.gz')) )
     orig = op.join(fs_dir, 'mri', 'orig', '001.mgz')
     for d in ds:
-        print "Processing %s:" % d[0]
+        iflogger.info("Processing %s:" % d[0])
         if not op.exists(d[0]):
             raise Exception('File %s does not exist.' % d[0])
         # reslice to original volume because the roi creation with freesurfer
@@ -382,8 +386,8 @@ class Parcellate(BaseInterface):
     def _run_interface(self, runtime):
         if self.inputs.subjects_dir:
            os.environ.update({'SUBJECTS_DIR': self.inputs.subjects_dir})
-        print "ROI_HR_th.nii.gz / fsmask_1mm.nii.gz CREATION"
-        print "============================================="
+        iflogger.info("ROI_HR_th.nii.gz / fsmask_1mm.nii.gz CREATION")
+        iflogger.info("=============================================")
         create_annot_label(self.inputs.subject_id, self.inputs.subjects_dir, self.inputs.freesurfer_dir, self.inputs.parcellation_name)
         create_roi(self.inputs.subject_id, self.inputs.subjects_dir, self.inputs.freesurfer_dir, self.inputs.parcellation_name)
         create_wm_mask(self.inputs.subject_id, self.inputs.subjects_dir, self.inputs.freesurfer_dir, self.inputs.parcellation_name)
@@ -396,6 +400,11 @@ def _list_outputs(self):
             outputs['roi_file'] = op.abspath(self.inputs.out_roi_file)
         else:
             outputs['roi_file'] = op.abspath(self._gen_outfilename('nii.gz', 'ROI'))
+        outputs['white_matter_mask_file'] = op.abspath('fsmask_1mm.nii.gz')
+        outputs['cc_unknown_file'] = op.abspath('cc_unknown.nii.gz')
+        outputs['ribbon_file'] = op.abspath('ribbon.nii.gz')
+        outputs['aseg_file'] = op.abspath('aseg.nii.gz')
+        outputs['ROI_HR_th_file'] = op.abspath('ROI_HR_th.nii.gz')
         return outputs
 
     def _gen_outfilename(self, ext, prefix='ROI'):
diff --git a/nipype/workflows/dmri/camino/connectivity_mapping.py b/nipype/workflows/dmri/camino/connectivity_mapping.py
@@ -221,6 +221,7 @@ def create_connectivity_pipeline(name="connectivity"):
     """
 
     roigen = pe.Node(interface=cmtk.ROIGen(), name="ROIGen")
+    roigen_structspace = roigen.clone("ROIGen_structspace")
 
     """
     The CreateMatrix interface takes in the remapped aparc+aseg image as well as the label dictionary and fiber tracts
@@ -231,6 +232,7 @@ def create_connectivity_pipeline(name="connectivity"):
     specific tracts that connect between user-selected regions.
     """
 
+    createnodes = pe.Node(interface=cmtk.CreateNodes(), name="CreateNodes")
     creatematrix = pe.Node(interface=cmtk.CreateMatrix(), name="CreateMatrix")
     creatematrix.inputs.count_region_intersections = True
 
@@ -382,12 +384,16 @@ def create_connectivity_pipeline(name="connectivity"):
     """
 
     mapping.connect(inputnode_within, 'resolution_network_file',
+                    createnodes, 'resolution_network_file')
+    mapping.connect(createnodes, 'node_network',
                     creatematrix, 'resolution_network_file')
     mapping.connect([(FreeSurferSource, mri_convert_AparcAseg, [(('aparc_aseg', select_aparc), 'in_file')])])
 
     mapping.connect([(b0Strip, inverse_AparcAseg,[('out_file','reference')])])
     mapping.connect([(convertxfm, inverse_AparcAseg,[('out_file','in_matrix_file')])])
     mapping.connect([(mri_convert_AparcAseg, inverse_AparcAseg,[('out_file','in_file')])])
+    mapping.connect([(mri_convert_AparcAseg, roigen_structspace,[('out_file','aparc_aseg_file')])])
+    mapping.connect([(roigen_structspace, createnodes,[("roi_file","roi_file")])])
 
     mapping.connect([(inverse_AparcAseg, roigen,[("out_file","aparc_aseg_file")])])
     mapping.connect([(roigen, creatematrix,[("roi_file","roi_file")])])
@@ -432,13 +438,13 @@ def create_connectivity_pipeline(name="connectivity"):
     CFFConverter.inputs.script_files = op.abspath(inspect.getfile(inspect.currentframe()))
     mapping.connect([(giftiSurfaces, CFFConverter,[("out","gifti_surfaces")])])
     mapping.connect([(giftiLabels, CFFConverter,[("out","gifti_labels")])])
-    mapping.connect([(creatematrix, CFFConverter,[("matrix_file","gpickled_networks")])])
+    mapping.connect([(creatematrix, CFFConverter,[("matrix_files","gpickled_networks")])])
 
     mapping.connect([(niftiVolumes, CFFConverter,[("out","nifti_volumes")])])
     mapping.connect([(fiberDataArrays, CFFConverter,[("out","data_files")])])
     mapping.connect([(camino2trackvis, CFFConverter,[("trackvis","tract_files")])])
     mapping.connect([(inputnode_within, CFFConverter,[("subject_id","title")])])
-
+	
     """
     Finally, we create another higher-level workflow to connect our mapping workflow with the info and datagrabbing nodes
     declared at the beginning. Our tutorial can is now extensible to any arbitrary number of subjects by simply adding
@@ -453,7 +459,7 @@ def create_connectivity_pipeline(name="connectivity"):
                                                                 "tracts",
                                                                 "connectome",
                                                                 "cmatrix",
-                                                                "gpickled_network",
+                                                                "networks",
                                                                 "rois",
                                                                 "mean_fiber_length",
                                                                 "fiber_length_std",
@@ -477,7 +483,7 @@ def create_connectivity_pipeline(name="connectivity"):
         ("CreateMatrix.mean_fiber_length_matrix_mat_file", "mean_fiber_length"),
         ("CreateMatrix.fiber_length_std_matrix_mat_file", "fiber_length_std"),
         ("fa2nii.nifti_file", "fa"),
-        ("CreateMatrix.matrix_file", "gpickled_network"),
+        ("CreateMatrix.matrix_files", "networks"),
         ("ROIGen.roi_file", "rois"),
         ("mri_convert_Brain.out_file", "struct"),
         ("trace2nii.nifti_file", "trace"),
diff --git a/nipype/workflows/dmri/camino/group_connectivity.py b/nipype/workflows/dmri/camino/group_connectivity.py
@@ -90,7 +90,7 @@ def create_group_connectivity_pipeline(group_list, group_id, data_dir, subjects_
                                               ("outputnode.cmatrix", "@l1output.cmatrix"),
                                               ("outputnode.rois", "@l1output.rois"),
                                               ("outputnode.struct", "@l1output.struct"),
-                                              ("outputnode.gpickled_network", "@l1output.gpickled_network"),
+                                              ("outputnode.networks", "@l1output.networks"),
                                               ("outputnode.mean_fiber_length", "@l1output.mean_fiber_length"),
                                               ("outputnode.fiber_length_std", "@l1output.fiber_length_std"),
                                               ])])
diff --git a/nipype/workflows/dmri/connectivity/nx.py b/nipype/workflows/dmri/connectivity/nx.py
