WIP: Dmri nipype

examples/dmri_tractography.py

@@ -0,0 +1,75 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 16 20:39:00 2013
+
+@author: bao
+"""
+import os
+import nipype.interfaces.io as nio
+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 EddyCorrection, ResampleVoxelSize, TensorModel, Tracking
+
+path ='/home/bao/tiensy/Nipype_tutorial/data/dmri/temp8/'
+data = path+ 'raw.nii.gz'
+
+###### WORKFLOW DEFINITION #######
+wf=pe.Workflow(name="reconstructing_tractography")
+wf.base_dir= path + 'results'
+wf.config['execution'] = {'remove_unnecessary_outputs': 'True'}
+
+###### NODE DEFINITION #######
+datasink = pe.Node(nio.DataSink(), name='sinker')
+datasink.inputs.base_directory = path + 'out/'
+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.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.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')
+
+wf.connect(brain_extraction_node, 'out_file', datasink, 'raw_bet.@nii.@gz')
+wf.connect(eddy_current_correction_node, 'outputnode.eddy_corrected', datasink, 'raw_bet_ecc.@nii.@gz')
+wf.connect(resample_voxel_size_node, 'resample_file', datasink, 'data_bet_ecc_iso.@nii.@gz')
+wf.connect(tensor_model_node, 'tensor_fa_file', datasink, 'tensor_fa.@nii.@gz')
+wf.connect(tensor_model_node, 'tensor_evecs_file', datasink, 'tensor_evecs.@nii.@gz')
+wf.connect(tracking_node, 'tracks_file', datasink, 'dti_tracks.@dpy')
+
+###### GRAPH and RUN #######
+wf.write_graph()
+wf.run()

nipype/interfaces/dipy/dmri_classInterfaces.py

@@ -0,0 +1,133 @@
+# -*- 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 eddy_correction, resample_voxel_size
+from tracking import tensor_model, tracking
+
+from nipype import logging
+iflogger = logging.getLogger('interface')
+
+from sys import stdout
+import os
+
+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):
+        iflogger.info('Eddy current correction ...')
+        self._out_file =eddy_correction(self.inputs.input_filename,self.inputs.output_filename)
+        iflogger.info('Saving to: {filename}'.format(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):
+        iflogger.info('Resample...')      
+        iflogger.info('Loading data: {filename}'.format(filename=self.inputs.input_filename))
+        self._out_file = resample_voxel_size(self.inputs.input_filename,self.inputs.output_filename)
+        iflogger.info('Saving data after resapling to {filename}'.format(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):
+        iflogger.info('Tensor model ...')    
+        iflogger.info('Loading data ...')
+        (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):
+        iflogger.info('Tracking ...')
+        iflogger.info('Loading data ...')   
+        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)                                
+        iflogger.info('Saving tracks to: {filename}'.format(filename=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/dipy/dmri_preprocessing.py

@@ -0,0 +1,53 @@
+# -*- 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.)
+"""
+
+
+def eddy_correction(input_filename, output_filename=None): 
+
+    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)   
+
+    return filename_save
+
+def resample_voxel_size(input_filename, output_filename=None): 
+
+    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  
+
+
+    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)       
+
+
+    data_img = nib.Nifti1Image(data=data, affine=affine)
+    nib.save(data_img, filename_save)
+
+    return filename_save
+
+
+
+

nipype/interfaces/dipy/dmri_tracking.py

@@ -0,0 +1,117 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 16 14:00:54 2013
+
+@author: bao
+Create and save the tractography from the iso-data (size of voxel must be isotrophy)
+
+"""
+import os
+
+import nibabel as nib
+import numpy as np
+
+from dipy.io.dpy import Dpy
+from dipy.io.gradients import read_bvals_bvecs
+from dipy.data import get_sphere
+from dipy.reconst.dti import quantize_evecs
+from dipy.reconst.dti import TensorModel
+from dipy.reconst.dti import fractional_anisotropy
+from dipy.core.gradients import gradient_table
+from dipy.tracking.eudx import EuDX
+
+
+
+'''
+Create the Tensor Model for data
+input:  data file name, bvecs file name and bvals file name
+output: tensor_fa and tensor_evecs
+'''
+def tensor_model(input_filename_data, input_filename_bvecs, input_filename_bvals, output_filename_fa=None, output_filename_evecs=None): 
+    img = nib.load(input_filename_data)
+    data = img.get_data()
+    affine = img.get_affine()
+
+    bvals, bvecs = read_bvals_bvecs(input_filename_bvals, input_filename_bvecs)
+    gtab = gradient_table(bvals, bvecs)
+
+    mask = data[..., 0] > 50
+    tenmodel = TensorModel(gtab)
+    tenfit = tenmodel.fit(data, mask)
+
+    FA = fractional_anisotropy(tenfit.evals)
+    FA[np.isnan(FA)] = 0
+
+    if output_filename_fa == None:
+        filename_save_fa = input_filename_data.split('.')[0]+'_tensor_fa.nii.gz'
+    else:
+        filename_save_fa = os.path.abspath(output_filename_fa)    
+
+    fa_img = nib.Nifti1Image(FA, img.get_affine())
+    nib.save(fa_img, filename_save_fa)        
+
+    if output_filename_evecs == None:
+        filename_save_evecs = input_filename_data.split('.')[0]+'_tensor_evecs.nii.gz'
+    else:
+        filename_save_evecs = os.path.abspath(output_filename_evecs) 
+
+    evecs_img = nib.Nifti1Image(tenfit.evecs, img.get_affine())
+    nib.save(evecs_img, filename_save_evecs)    
+
+    return filename_save_fa,filename_save_evecs
+
+
+def create_tracks(anisotropy, indices, vertices , seeds, low_thresh):
+
+    eu = EuDX(anisotropy, indices, odf_vertices = vertices, seeds = seeds, a_low=low_thresh)
+    print eu.seed_no    
+
+    tensor_tracks_old = [streamline for streamline in eu]
+    print len(tensor_tracks_old)
+
+    #remove one point tracks
+    tracks = [track for track in tensor_tracks_old if track.shape[0]>1]    
+
+    return tracks
+
+'''
+Tracking to reconstruct the tractography
+input:  tensor_fa and tensor_evecs file name
+        num of seeds        
+output: tractography file name
+'''
+def tracking(input_filename_fa, input_filename_evecs, num_seeds=10000, low_thresh = .2, output_filename=None):    
+
+    fa_img = nib.load(input_filename_fa)
+    FA = fa_img.get_data()
+
+    evecs_img = nib.load(input_filename_evecs)
+    evecs = evecs_img.get_data()
+
+    FA[np.isnan(FA)] = 0    
+
+    sphere = get_sphere('symmetric724')    
+    peak_indices = quantize_evecs(evecs, sphere.vertices)        
+
+    eu = EuDX(FA, peak_indices, seeds = num_seeds, odf_vertices = sphere.vertices,  a_low=low_thresh)
+
+    tensor_tracks_old = [streamline for streamline in eu]   
+
+    #remove one point tracks
+    tracks = [track for track in tensor_tracks_old if track.shape[0]>1]         
+
+    if output_filename == None:
+        filename_save = 'tracks_dti.dpy'
+    else:
+        filename_save = os.path.abspath(output_filename)        
+
+    dpw = Dpy(filename_save, 'w')
+    dpw.write_tracks(tracks)
+    dpw.close()
+
+    return filename_save
+
+
+
+
+