diff --git a/examples/dmri_tractography.py b/examples/dmri_tractography.py new file mode 100644 index 0000000000..7f79ba1977 --- /dev/null +++ b/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() diff --git a/nipype/interfaces/dipy/dmri_classInterfaces.py b/nipype/interfaces/dipy/dmri_classInterfaces.py new file mode 100644 index 0000000000..7a8c150e0f --- /dev/null +++ b/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 + +### + diff --git a/nipype/interfaces/dipy/dmri_preprocessing.py b/nipype/interfaces/dipy/dmri_preprocessing.py new file mode 100644 index 0000000000..18fa4e1dc5 --- /dev/null +++ b/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 + + + + diff --git a/nipype/interfaces/dipy/dmri_tracking.py b/nipype/interfaces/dipy/dmri_tracking.py new file mode 100644 index 0000000000..7a27208030 --- /dev/null +++ b/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 + + + + + diff --git a/nipype/interfaces/dmri_classInterfaces.py b/nipype/interfaces/dmri_classInterfaces.py new file mode 100644 index 0000000000..048d76cc03 --- /dev/null +++ b/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 + +### \ No newline at end of file diff --git a/nipype/interfaces/dmri_preprocessing.py b/nipype/interfaces/dmri_preprocessing.py new file mode 100644 index 0000000000..75e7f4ade7 --- /dev/null +++ b/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 + + + + \ No newline at end of file diff --git a/nipype/interfaces/dmri_tracking.py b/nipype/interfaces/dmri_tracking.py new file mode 100644 index 0000000000..a1f2649f98 --- /dev/null +++ b/nipype/interfaces/dmri_tracking.py @@ -0,0 +1,129 @@ +# -*- 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): + + print 'Tensor model ...' + + print 'Loading data ...' + 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) + print "Saving fa to:", 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) + print "Saving evecs to:", 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 + #stop + 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): + + print 'Tracking ...' + + print 'Loading data ...' + + 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] + #tracks = create_tracks(FA, peak_indices,sphere.vertices, num_seeds, low_thresh) + + if output_filename == None: + filename_save = 'tracks_dti.dpy' + else: + filename_save = os.path.abspath(output_filename) + + print 'Saving tracks to:', filename_save + dpw = Dpy(filename_save, 'w') + dpw.write_tracks(tracks) + dpw.close() + + return filename_save + + + + +