WIP: Dmri nipype

examples/dmri_tractography.py

@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 16 20:39:00 2013
+
+@author: bao
+"""
+import os
+import nipype.interfaces.utility as util
+import nipype.pipeline.engine as pe
+
+from nipype.interfaces import fsl
+from nipype.workflows.dmri.fsl.dti import create_eddy_correct_pipeline
+
+from dmri_classInterfaces import BrainExtraction, EddyCorrection, ResampleVoxelSize, TensorModel, Tracking
+
+path ='/home/bao/tiensy/Nipype_tutorial/data/dmri/temp6/'
+data = path+ 'raw.nii.gz'
+
+###### WORKFLOW DEFINITION #######
+wf=pe.Workflow(name="reconstructing_tractography")
+wf.base_dir= path + 'results'
+wf.config['execution'] = {'remove_unnecessary_outputs': 'False',
+                            }
+
+
+###### NODE DEFINITION #######
+brain_extraction_node = pe.Node(fsl.BET(), name="brain_extraction_node")
+eddy_current_correction_node = create_eddy_correct_pipeline("nipype_eddycorrect_wkf")
+resample_voxel_size_node = pe.Node(ResampleVoxelSize(), name='resample_voxel_size_node')
+tensor_model_node = pe.Node(TensorModel(), name='tensor_model_node')
+tracking_node = pe.Node(Tracking(), name='tracking_node')
+
+
+###### INPUT NODE DEFINITION #######
+#inputs: brain_extraction_node
+brain_extraction_node.inputs.in_file=data
+brain_extraction_node.inputs.frac = 0.2
+brain_extraction_node.inputs.functional = True
+#brain_extraction_node.inputs.robust = True
+brain_extraction_node.inputs.vertical_gradient = 0
+brain_extraction_node.inputs.out_file = path + 'raw_bet.nii.gz'
+
+
+#inputs: eddy_current_correction_node
+#eddy_current_correction_node.inputs.inputnode.in_file = path + 'raw_bet.nii.gz'
+eddy_current_correction_node.inputs.inputnode.ref_num = 0
+
+#inputs: resample_voxel_size_node
+resample_voxel_size_node.inputs.output_filename = path + 'data_bet_ecc_iso.nii.gz'
+
+
+#inputs: tensor_model_node
+tensor_model_node.inputs.input_filename_bvecs = path + 'raw.bvec'
+tensor_model_node.inputs.input_filename_bvals = path + 'raw.bval'
+
+tensor_model_node.inputs.output_filename_fa = path + 'tensor_fa.nii.gz'
+tensor_model_node.inputs.output_filename_evecs = path + 'tensor_evecs.nii.gz'
+
+#inputs: tracking_node
+tracking_node.inputs.num_seeds = 1000
+tracking_node.inputs.low_thresh = 0.2
+tracking_node.inputs.output_filename = path + 'dti_tracks.dpy'
+
+###### NODE CONNECTIONS #######
+wf.connect(brain_extraction_node,'out_file', eddy_current_correction_node, 'inputnode.in_file')
+wf.connect(eddy_current_correction_node,'outputnode.eddy_corrected', resample_voxel_size_node ,'input_filename')
+wf.connect(resample_voxel_size_node,'resample_file',tensor_model_node ,'input_filename_data')
+wf.connect(tensor_model_node, 'tensor_fa_file', tracking_node,'input_filename_fa')
+wf.connect(tensor_model_node, 'tensor_evecs_file', tracking_node , 'input_filename_evecs')
+
+
+###### GRAPH and RUN #######
+wf.write_graph()
+wf.run()

nipype/interfaces/dmri_classInterfaces.py

@@ -0,0 +1,147 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 16 18:43:39 2013
+
+@author: bao
+"""
+
+from nipype.interfaces.base import BaseInterface, BaseInterfaceInputSpec, traits, File, TraitedSpec
+from nipype.utils.filemanip import split_filename
+
+from preprocessing import brain_extraction, eddy_correction, resample_voxel_size
+from tracking import tensor_model, tracking
+
+from nipype import logging
+iflogger = logging.getLogger('interface')
+
+import nibabel as nb
+import numpy as np
+from sys import stdout
+import os
+
+class BrainExtractionInputSpec(BaseInterfaceInputSpec):                                                     
+    input_filename = File(exists=True,desc="Nifti file to be processed",mandatory=True)
+    output_filename = File(exists=False,desc="Output file name",mandatory=False)
+
+class BrainExtractionOutputSpec(TraitedSpec):
+    bet_file = File(exists=True ,desc="Output file name")
+
+
+class BrainExtraction(BaseInterface):
+    input_spec = BrainExtractionInputSpec
+    output_spec = BrainExtractionOutputSpec
+
+    def _run_interface(self, runtime):
+        self._out_file = brain_extraction(self.inputs.input_filename,self.inputs.output_filename)
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs["bet_file"] = os.path.abspath(self._out_file)
+        return outputs
+
+###
+
+class EddyCorrectionInputSpec(BaseInterfaceInputSpec):                                                     
+    input_filename = File(exists=True,desc="Nifti file to be processed",mandatory=True)
+    output_filename = File(exists=False,desc="Output file name",mandatory=False)
+
+class EddyCorrectionOutputSpec(TraitedSpec):
+    eddy_current_correction_file = File(exists=True ,desc="Output file name")
+
+
+class EddyCorrection(BaseInterface):
+    input_spec = EddyCorrectionInputSpec
+    output_spec = EddyCorrectionOutputSpec
+
+    def _run_interface(self, runtime):
+        self._out_file =eddy_correction(self.inputs.input_filename,self.inputs.output_filename)
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs["eddy_current_correction_file"] = os.path.abspath(self._out_file)
+        return outputs
+
+###
+
+class ResampleVoxelSizeInputSpec(BaseInterfaceInputSpec):                                                     
+    input_filename = File(exists=True,desc="Nifti file to be processed",mandatory=True)
+    output_filename = File(exists=False,desc="Output file name",mandatory=False)
+
+class ResampleVoxelSizeOutputSpec(TraitedSpec):
+    resample_file = File(exists=True ,desc="Output file name")
+
+
+class ResampleVoxelSize(BaseInterface):
+    input_spec = ResampleVoxelSizeInputSpec
+    output_spec = ResampleVoxelSizeOutputSpec
+
+    def _run_interface(self, runtime):
+        self._out_file = resample_voxel_size(self.inputs.input_filename,self.inputs.output_filename)
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs["resample_file"] = os.path.abspath(self._out_file)
+        return outputs
+
+###
+
+class TensorModelInputSpec(BaseInterfaceInputSpec):                                                     
+    input_filename_data = File(exists=True,desc="Nifti file to be processed",mandatory=True)
+    input_filename_bvecs = File(exists=True,desc="bvec file",mandatory=True)
+    input_filename_bvals = File(exists=True,desc="bval file",mandatory=True) 
+    output_filename_fa = File(exists=False,desc="Output fa file name",mandatory=False)
+    output_filename_evecs = File(exists=False,desc="Output evecs file name",mandatory=False)
+
+class TensorModelOutputSpec(TraitedSpec):
+    tensor_fa_file = File(exists=True ,desc="Output fa file name")
+    tensor_evecs_file = File(exists=True ,desc="Output evecs file name")
+
+
+class TensorModel(BaseInterface):
+    input_spec = TensorModelInputSpec
+    output_spec = TensorModelOutputSpec
+
+    def _run_interface(self, runtime):
+        (self.fa_file,self.evecs_file)  = tensor_model(self.inputs.input_filename_data, self.inputs.input_filename_bvecs,
+                                      self.inputs.input_filename_bvals, self.inputs.output_filename_fa,
+                                      self.inputs.output_filename_evecs)
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs['tensor_fa_file'] = os.path.abspath(self.fa_file)
+        outputs['tensor_evecs_file'] = os.path.abspath(self.evecs_file)
+        return outputs
+
+###
+
+class TrackingInputSpec(BaseInterfaceInputSpec):                                                     
+    input_filename_fa = File(exists=True,desc="FA file to be processed",mandatory=True)
+    input_filename_evecs = File(exists=True,desc="Evecs file to be processed",mandatory=True)
+    num_seeds = traits.Long(desc="Num of seeds for initializing the position of tracks",mandatory=False)
+    low_thresh = traits.Float(desc="Lower threshold for  FA, typical 0.2 ",mandatory=False)
+    output_filename = File(exists=False,desc="Output file name",mandatory=False)
+
+class TrackingOutputSpec(TraitedSpec):
+    tracks_file = File(exists=True ,desc="Output file name")
+
+
+class Tracking(BaseInterface):
+    input_spec = TrackingInputSpec
+    output_spec = TrackingOutputSpec
+
+    def _run_interface(self, runtime):
+        self._out_file = tracking(self.inputs.input_filename_fa, self.inputs.input_filename_evecs,
+                                  self.inputs.num_seeds, self.inputs.low_thresh, 
+                                  self.inputs.output_filename)
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs["tracks_file"] = os.path.abspath(self._out_file)
+        return outputs
+
+###

nipype/interfaces/dmri_preprocessing.py

@@ -0,0 +1,78 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 16 13:23:24 2013
+
+@author: bao
+Pre-processing dMRI data, including 03 steps
+1. Brain extraction
+2. Eddy current correction
+3. Resample data to isotrophy, usually to voxel size (2.,2.,2.)
+"""
+import os
+
+import nibabel as nib
+from dipy.external.fsl import bet, eddy_correct
+from dipy.align.aniso2iso import resample
+
+
+
+def brain_extraction(input_filename, output_filename=None):
+
+    print 'Brain extraction ...'
+
+    if output_filename == None:
+        filename_save = input_filename.split('.')[0]+'_bet.nii.gz'
+    else:
+        filename_save = os.path.abspath(output_filename)
+
+    bet(input_filename, filename_save,options=' -R -F -f .2 -g 0')
+
+    print "Saving to:", filename_save    
+
+    return filename_save
+
+def eddy_correction(input_filename, output_filename=None): 
+
+    print 'Eddy current correction ...'
+
+    if output_filename == None:
+        filename_save = input_filename.split('.')[0]+'_ecc.nii.gz'
+    else:
+        filename_save = os.path.abspath(output_filename)
+
+    eddy_correct(input_filename,filename_save)
+
+    print "Saving to:", filename_save    
+
+    return filename_save
+
+def resample_voxel_size(input_filename, output_filename=None): 
+
+    print("Loading data: %s" % input_filename)
+    img = nib.load(input_filename)
+
+    old_data = img.get_data()
+    old_affine = img.get_affine()
+
+    zooms=img.get_header().get_zooms()[:3]
+    print 'Old zooms:', zooms
+    new_zooms=(2.,2.,2.)
+    print 'New zoom', new_zooms  
+
+    print 'Resample data and affine ...'      
+    data,affine=resample(old_data,old_affine,zooms,new_zooms)
+
+    if output_filename == None:
+        filename_save = input_filename.split('.')[0]+'_iso.nii.gz'
+    else:
+        filename_save = os.path.abspath(output_filename)       
+
+    print "Saving data after resapling to ", filename_save
+    data_img = nib.Nifti1Image(data=data, affine=affine)
+    nib.save(data_img, filename_save)
+
+    return filename_save
+
+
+
+