diff --git a/examples/diffusion_toolkit_tutorial.py b/examples/dtk_dti_tutorial.py similarity index 62% rename from examples/diffusion_toolkit_tutorial.py rename to examples/dtk_dti_tutorial.py index fa3adfa935..b553c53157 100644 --- a/examples/diffusion_toolkit_tutorial.py +++ b/examples/dtk_dti_tutorial.py @@ -2,7 +2,7 @@ """ A pipeline example that uses several interfaces to perform analysis on diffusion weighted images using -FSL FDT tools. +Diffusion Toolkit and FSL. This tutorial is based on the 2010 FSL course and uses data freely available at the FSL website at: @@ -69,15 +69,7 @@ info = dict(dwi=[['subject_id', 'data']], bvecs=[['subject_id','bvecs']], - bvals=[['subject_id','bvals']], - seed_file = [['subject_id','MASK_average_thal_right']], - target_masks = [['subject_id',['MASK_average_M1_right', - 'MASK_average_S1_right', - 'MASK_average_occipital_right', - 'MASK_average_pfc_right', - 'MASK_average_pmc_right', - 'MASK_average_ppc_right', - 'MASK_average_temporal_right']]]) + bvals=[['subject_id','bvals']]) infosource = pe.Node(interface=util.IdentityInterface(fields=['subject_id']), name="infosource") @@ -111,9 +103,7 @@ # http://www.fmrib.ox.ac.uk/fslcourse/fsl_course_data2.tar.gz datasource.inputs.base_directory = os.path.abspath('fsl_course_data/fdt/') -datasource.inputs.field_template = dict(dwi='%s/%s.nii.gz', - seed_file="%s.bedpostX/%s.nii.gz", - target_masks="%s.bedpostX/%s.nii.gz") +datasource.inputs.field_template = dict(dwi='%s/%s.nii.gz') datasource.inputs.template_args = info @@ -160,7 +150,7 @@ computeTensor.connect([ (fslroi,bet,[('roi_file','in_file')]), - (eddycorrect,dtifit,[('eddy_corrected','dwi')]) + (eddycorrect,dtifit,[('eddy_corrected','DWI')]) ]) @@ -168,61 +158,23 @@ """ Setup for Tracktography ----------------------- -Here we will create a workflow to enable probabilistic tracktography -and hard segmentation of the seed region +Here we will create a workflow to enable deterministic tracktography """ tractography = pe.Workflow(name='tractography') -tractography.base_dir = os.path.abspath('diffusion_toolkit_tutorial') - -""" -estimate the diffusion parameters: phi, theta, and so on -""" - -bedpostx = pe.Node(interface=fsl.BEDPOSTX(),name='bedpostx') -bedpostx.inputs.fibres = 1 - -bedpostx_2f = pe.Node(interface=fsl.BEDPOSTX(),name='bedpostx_2f') -bedpostx_2f.inputs.fibres = 2 - - -flirt = pe.Node(interface=fsl.FLIRT(), name='flirt') -flirt.inputs.reference = fsl.Info.standard_image('MNI152_T1_2mm_brain.nii.gz') -flirt.inputs.dof = 12 - -""" -perform probabilistic tracktography -""" - -probtrackx = pe.Node(interface=fsl.ProbTrackX(),name='probtrackx') -probtrackx.inputs.mode='seedmask' -probtrackx.inputs.loop_check=True -probtrackx.inputs.c_thresh = 0.2 -probtrackx.inputs.n_steps=2000 -probtrackx.inputs.step_length=0.5 -probtrackx.inputs.n_samples=5000 -probtrackx.inputs.force_dir=True -probtrackx.inputs.opd=True -probtrackx.inputs.os2t=True - - -""" -perform hard segmentation on the output of probtrackx -""" - -findthebiggest = pe.Node(interface=fsl.FindTheBiggest(),name='findthebiggest') +dtk_tracker = pe.Node(interface=dtk.DTITracker(), name="dtk_tracker") +dtk_tracker.inputs.invert_x = True +smooth_trk = pe.Node(interface=dtk.SplineFilter(), name="smooth_trk") +smooth_trk.inputs.step_length = 0.5 """ connect all the nodes for this workflow """ tractography.connect([ - (bedpostx,probtrackx,[('bpx_out_directory','bpx_directory')]), - (bedpostx,probtrackx,[('bpx_out_directory','out_dir')]), - (probtrackx,findthebiggest,[('targets','in_files')]), - (flirt, probtrackx, [('out_matrix_file','xfm')]) - ]) + (dtk_tracker, smooth_trk, [('track_file', 'track_file')]) + ]) """ @@ -242,24 +194,16 @@ def getstripdir(subject_id): """ dwiproc = pe.Workflow(name="dwiproc") -dwiproc.base_dir = os.path.abspath('dti_tutorial') +dwiproc.base_dir = os.path.abspath('dtk_dti_tutorial') dwiproc.connect([ (infosource,datasource,[('subject_id', 'subject_id')]), (datasource,computeTensor,[('dwi','fslroi.in_file'), ('bvals','dtifit.bvals'), ('bvecs','dtifit.bvecs'), ('dwi','eddycorrect.in_file')]), - (datasource,tractography,[('bvals','bedpostx.bvals'), - ('bvecs','bedpostx.bvecs'), - ('seed_file','probtrackx.seed_file'), - ('target_masks','probtrackx.target_masks')]), - (computeTensor,tractography,[('eddycorrect.eddy_corrected','bedpostx.dwi'), - ('bet.mask_file','bedpostx.mask'), - ('bet.mask_file','probtrackx.mask'), - ('fslroi.roi_file','flirt.in_file')]), - (infosource, datasink,[('subject_id','container'), - (('subject_id', getstripdir),'strip_dir')]), - (tractography,datasink,[('findthebiggest.out_file','fbiggest.@biggestsegmentation')]) + (computeTensor,tractography,[('bet.mask_file','dtk_tracker.mask1_file'), + ('dtifit.tensor','dtk_tracker.tensor_file') + ]) ]) dwiproc.run() diff --git a/examples/dtk_odf_tutorial.py b/examples/dtk_odf_tutorial.py new file mode 100644 index 0000000000..47d47d2860 --- /dev/null +++ b/examples/dtk_odf_tutorial.py @@ -0,0 +1,203 @@ + +""" +A pipeline example that uses several interfaces to +perform analysis on diffusion weighted images using +Diffusion Toolkit and FSL. + +This tutorial uses data from out nipype-tutorial package. +""" + + +""" +Tell python where to find the appropriate functions. +""" + +import nipype.interfaces.io as nio # Data i/o +import nipype.interfaces.fsl as fsl # fsl +import nipype.interfaces.diffusion_toolkit as dtk +import nipype.interfaces.utility as util # utility +import nipype.pipeline.engine as pe # pypeline engine +import os # system functions + +""" +Confirm package dependencies are installed. (This is only for the +tutorial, rarely would you put this in your own code.) +""" + +from nipype.utils.misc import package_check + +package_check('numpy', '1.3', 'tutorial1') +package_check('scipy', '0.7', 'tutorial1') +package_check('networkx', '1.0', 'tutorial1') +package_check('IPython', '0.10', 'tutorial1') + + +""" +Setting up workflows +-------------------- +This is a generic workflow for DTI data analysis using the FSL +""" + +""" +Data specific components +------------------------ + +The nipype tutorial contains data for two subjects. Subject data +is in two subdirectories, ``dwis1`` and ``dwis2``. Each subject directory +contains each of the following files: bvec, bval, diffusion weighted data, a set of target masks, +a seed file, and a transformation matrix. + +Below we set some variables to inform the ``datasource`` about the +layout of our data. We specify the location of the data, the subject +sub-directories and a dictionary that maps each run to a mnemonic (or +field) for the run type (``dwi`` or ``bvals``). These fields become +the output fields of the ``datasource`` node in the pipeline. + +Specify the subject directories +""" + +subject_list = ['siemens_hardi_test'] + + +""" +Map field names to individual subject runs +""" + +info = dict(dwi=[['subject_id', 'siemens_hardi_test_data']], + bvecs=[['subject_id','siemens_hardi_test_data.bvec']], + bvals=[['subject_id','siemens_hardi_test_data.bval']]) + +infosource = pe.Node(interface=util.IdentityInterface(fields=['subject_id']), + name="infosource") + +"""Here we set up iteration over all the subjects. The following line +is a particular example of the flexibility of the system. The +``datasource`` attribute ``iterables`` tells the pipeline engine that +it should repeat the analysis on each of the items in the +``subject_list``. In the current example, the entire first level +preprocessing and estimation will be repeated for each subject +contained in subject_list. +""" + +infosource.iterables = ('subject_id', subject_list) + +""" +Now we create a :class:`nipype.interfaces.io.DataGrabber` object and +fill in the information from above about the layout of our data. The +:class:`nipype.pipeline.engine.Node` module wraps the interface object +and provides additional housekeeping and pipeline specific +functionality. +""" + +datasource = pe.Node(interface=nio.DataGrabber(infields=['subject_id'], + outfields=info.keys()), + name = 'datasource') + +datasource.inputs.template = "%s/%s" + +# This needs to point to the fdt folder you can find after extracting +# http://www.fmrib.ox.ac.uk/fslcourse/fsl_course_data2.tar.gz +datasource.inputs.base_directory = os.path.abspath('data') + +datasource.inputs.field_template = dict(dwi='%s/%s.nii') +datasource.inputs.template_args = info + + +""" +Setup for ODF Computation +-------------------------------------- +Here we will create a generic workflow for ODF computation +""" + +compute_ODF = pe.Workflow(name='compute_ODF') + +""" +extract the volume with b=0 (nodif_brain) +""" + +fslroi = pe.Node(interface=fsl.ExtractROI(),name='fslroi') +fslroi.inputs.t_min=0 +fslroi.inputs.t_size=1 + +""" +create a brain mask from the nodif_brain +""" + +bet = pe.Node(interface=fsl.BET(),name='bet') +bet.inputs.mask=True +bet.inputs.frac=0.34 + +""" +correct the diffusion weighted images for eddy_currents +""" + +eddycorrect = pe.Node(interface=fsl.EddyCorrect(),name='eddycorrect') +eddycorrect.inputs.ref_num=0 + + +hardi_mat = pe.Node(interface=dtk.HARDIMat(),name='hardi_mat') + +odf_recon = pe.Node(interface=dtk.ODFRecon(),name='odf_recon') + +""" +connect all the nodes for this workflow +""" + +compute_ODF.connect([ + (fslroi,bet,[('roi_file','in_file')]), + (eddycorrect, odf_recon,[('eddy_corrected','DWI')]), + (eddycorrect, hardi_mat,[('eddy_corrected','reference_file')]), + (hardi_mat, odf_recon, [('out_file', 'matrix')]) + ]) + + + +""" +Setup for Tracktography +----------------------- +Here we will create a workflow to enable deterministic tracktography +""" + +tractography = pe.Workflow(name='tractography') + +odf_tracker = pe.Node(interface=dtk.ODFTracker(), name="odf_tracker") + +smooth_trk = pe.Node(interface=dtk.SplineFilter(), name="smooth_trk") +smooth_trk.inputs.step_length = 1 +""" +connect all the nodes for this workflow +""" + +tractography.connect([ + (odf_tracker, smooth_trk, [('track_file', 'track_file')]) + ]) + + +""" +Setup the pipeline that combines the two workflows: tractography and compute_ODF +---------------------------------------------------------------------------------- +""" + +dwiproc = pe.Workflow(name="dwiproc") +dwiproc.base_dir = os.path.abspath('dtk_odf_tutorial') +dwiproc.connect([ + (infosource,datasource,[('subject_id', 'subject_id')]), + (datasource,compute_ODF,[('dwi','fslroi.in_file'), + ('bvals','hardi_mat.bvals'), + ('bvecs','hardi_mat.bvecs'), + ('dwi','eddycorrect.in_file')]), + (compute_ODF,tractography,[('bet.mask_file','odf_tracker.mask1_file'), + ('odf_recon.ODF','odf_tracker.ODF'), + ('odf_recon.max','odf_tracker.max') + ]) + ]) + +dwiproc.inputs.compute_ODF.hardi_mat.oblique_correction = True +dwiproc.inputs.compute_ODF.odf_recon.n_directions = 31 +dwiproc.inputs.compute_ODF.odf_recon.n_b0 = 5 +dwiproc.inputs.compute_ODF.odf_recon.n_output_directions = 181 + +dwiproc.run() +dwiproc.write_graph() + + diff --git a/examples/dti_tutorial.py b/examples/fsl_dti_tutorial.py similarity index 98% rename from examples/dti_tutorial.py rename to examples/fsl_dti_tutorial.py index 79d59de471..bd8d797b13 100644 --- a/examples/dti_tutorial.py +++ b/examples/fsl_dti_tutorial.py @@ -174,7 +174,7 @@ """ tractography = pe.Workflow(name='tractography') -tractography.base_dir = os.path.abspath('dti_tutorial') +tractography.base_dir = os.path.abspath('fsl_dti_tutorial') """ estimate the diffusion parameters: phi, theta, and so on @@ -243,7 +243,7 @@ def getstripdir(subject_id): """ dwiproc = pe.Workflow(name="dwiproc") -dwiproc.base_dir = os.path.abspath('dti_tutorial') +dwiproc.base_dir = os.path.abspath('fsl_dti_tutorial') dwiproc.connect([ (infosource,datasource,[('subject_id', 'subject_id')]), (datasource,computeTensor,[('dwi','fslroi.in_file'), diff --git a/nipype/interfaces/diffusion_toolkit/__init__.py b/nipype/interfaces/diffusion_toolkit/__init__.py index b7af039552..2a086ecb29 100644 --- a/nipype/interfaces/diffusion_toolkit/__init__.py +++ b/nipype/interfaces/diffusion_toolkit/__init__.py @@ -1 +1,3 @@ -from nipype.interfaces.diffusion_toolkit.preproc import DTIRecon \ No newline at end of file +from nipype.interfaces.diffusion_toolkit.postproc import SplineFilter +from nipype.interfaces.diffusion_toolkit.dti import (DTIRecon, DTITracker) +from nipype.interfaces.diffusion_toolkit.odf import (HARDIMat, ODFRecon, ODFTracker) diff --git a/nipype/interfaces/diffusion_toolkit/dti.py b/nipype/interfaces/diffusion_toolkit/dti.py new file mode 100644 index 0000000000..b508828f1e --- /dev/null +++ b/nipype/interfaces/diffusion_toolkit/dti.py @@ -0,0 +1,166 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Provides interfaces to various commands provided by diffusion toolkit + + Change directory to provide relative paths for doctests + >>> import os + >>> filepath = os.path.dirname( os.path.realpath( __file__ ) ) + >>> datadir = os.path.realpath(os.path.join(filepath, '../../testing/data')) + >>> os.chdir(datadir) + +""" +import re +from nipype.utils.filemanip import fname_presuffix, split_filename, copyfile +import os +from nipype.utils.misc import isdefined +__docformat__ = 'restructuredtext' + +from nipype.interfaces.base import (TraitedSpec, File, traits, CommandLine, + CommandLineInputSpec) + +class DTIReconInputSpec(CommandLineInputSpec): + DWI = File(desc='Input diffusion volume', argstr='%s',exists=True, mandatory=True,position=1) + out_prefix = traits.Str("dti", desc='Output file prefix', argstr='%s', usedefault=True,position=2) + output_type = traits.Enum('nii', 'analyze', 'ni1', 'nii.gz', argstr='-ot %s', desc='output file type', usedefault=True) + bvecs = File(exists=True, desc = 'b vectors file', + argstr='-gm %s', mandatory=True) + bvals = File(exists=True,desc = 'b values file', mandatory=True) + n_averages = traits.Int(desc='Number of averages', argstr='-nex %s') + image_orientation_vectors = traits.List(traits.Float(), minlen=6, maxlen=6, desc="""specify image orientation vectors. if just one argument given, +will treat it as filename and read the orientation vectors from +the file. if 6 arguments are given, will treat them as 6 float +numbers and construct the 1st and 2nd vector and calculate the 3rd +one automatically. +this information will be used to determine image orientation, +as well as to adjust gradient vectors with oblique angle when""", argstr="-iop %f") + oblique_correction = traits.Bool(desc="""when oblique angle(s) applied, some SIEMENS dti protocols do not +adjust gradient accordingly, thus it requires adjustment for correct +diffusion tensor calculation""", argstr="-oc") + b0_threshold = traits.Float(desc="""program will use b0 image with the given threshold to mask out high +background of fa/adc maps. by default it will calculate threshold +automatically. but if it failed, you need to set it manually.""", argstr="-b0_th") + + +class DTIReconOutputSpec(TraitedSpec): + ADC = File(exists=True) + B0 = File(exists=True) + L1 = File(exists=True) + L2 = File(exists=True) + L3 = File(exists=True) + exp = File(exists=True) + FA = File(exists=True) + FA_color = File(exists=True) + tensor = File(exists=True) + V1 = File(exists=True) + V2 = File(exists=True) + V3 = File(exists=True) + +class DTIRecon(CommandLine): + """Use dti_recon to generate tensors and other maps + """ + + input_spec=DTIReconInputSpec + output_spec=DTIReconOutputSpec + + _cmd = 'dti_recon' + + def _create_gradient_matrix(self, bvecs_file, bvals_file): + _gradient_matrix_file = 'gradient_matrix.txt' + bvals = [val for val in re.split('\s+', open(bvals_file).readline().strip())] + bvecs_f = open(bvecs_file) + bvecs_x = [val for val in re.split('\s+', bvecs_f.readline().strip())] + bvecs_y = [val for val in re.split('\s+', bvecs_f.readline().strip())] + bvecs_z = [val for val in re.split('\s+', bvecs_f.readline().strip())] + bvecs_f.close() + gradient_matrix_f = open(_gradient_matrix_file, 'w') + for i in range(len(bvals)): + gradient_matrix_f.write("%s, %s, %s, %s\n"%(bvecs_x[i], bvecs_y[i], bvecs_z[i], bvals[i])) + gradient_matrix_f.close() + return _gradient_matrix_file + + def _format_arg(self, name, spec, value): + if name == "bvecs": + new_val = self._create_gradient_matrix(self.inputs.bvecs, self.inputs.bvals) + return super(DTIRecon, self)._format_arg("bvecs", spec, new_val) + return super(DTIRecon, self)._format_arg(name, spec, value) + + def _list_outputs(self): + out_prefix = self.inputs.out_prefix + output_type = self.inputs.output_type + + outputs = self.output_spec().get() + outputs['ADC'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_adc.'+ output_type)) + outputs['B0'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_b0.'+ output_type)) + outputs['L1'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_e1.'+ output_type)) + outputs['L2'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_e2.'+ output_type)) + outputs['L3'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_e3.'+ output_type)) + outputs['exp'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_exp.'+ output_type)) + outputs['FA'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_fa.'+ output_type)) + outputs['FA_color'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_fa_color.'+ output_type)) + outputs['tensor'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_tensor.'+ output_type)) + outputs['V1'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_v1.'+ output_type)) + outputs['V2'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_v2.'+ output_type)) + outputs['V3'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_v3.'+ output_type)) + + return outputs + +class DTITrackerInputSpec(CommandLineInputSpec): + tensor_file = File(exists=True, desc="reconstructed tensor file") + input_type = traits.Enum('nii', 'analyze', 'ni1', 'nii.gz', desc="""input and output file type. accepted values are: +analyze -> analyze format 7.5 +ni1 -> nifti format saved in seperate .hdr and .img file +nii -> nifti format with one .nii file +nii.gz -> nifti format with compression +default type is 'nii'""", argstr = "-it %s") + tracking_method = traits.Enum('fact', 'rk2', 'tl', 'sl', desc="""fact -> use FACT method for tracking. this is the default method. +rk2 -> use 2nd order runge-kutta method for tracking. +tl -> use tensorline method for tracking. +sl -> use interpolated streamline method with fixed step-length""", argstr="-%s") + step_length = traits.Float(desc="""set step length, in the unit of minimum voxel size. +default value is 0.5 for interpolated streamline method +and 0.1 for other methods""", argstr="-l %f") + angle_threshold = traits.Float(desc="set angle threshold. default value is 35 degree", argstr="-at %f") + angle_threshold_weight = traits.Float(desc="set angle threshold weighting factor. weighting will be be applied \ +on top of the angle_threshold", argstr = "-atw %f") + random_seed = traits.Int(desc = "use random location in a voxel instead of the center of the voxel \ + to seed. can also define number of seed per voxel. default is 1", argstr="-rseed") + invert_x = traits.Bool(desc="invert x component of the vector", argstr = "-ix") + invert_y = traits.Bool(desc="invert y component of the vector", argstr = "-iy") + invert_z = traits.Bool(desc="invert z component of the vector", argstr = "-iz") + swap_xy = traits.Bool(desc="swap x & y vectors while tracking", argstr = "-sxy") + swap_yz = traits.Bool(desc="swap y & z vectors while tracking", argstr = "-syz") + swap_zx = traits.Bool(desc="swap x & z vectors while tracking", argstr = "-szx") + mask1_file = File(desc="first mask image", mandatory=True, argstr="-m %s", position=2) + mask1_threshold = traits.Float(desc="threshold value for the first mask image, if not given, the program will \ +try automatically find the threshold", position=3) + mask2_file = File(desc="second mask image", argstr="-m2 %s", position=4) + mask2_threshold = traits.Float(desc="threshold value for the second mask image, if not given, the program will \ +try automatically find the threshold", position=5) + input_data_prefix = traits.Str("dti", desc="for internal naming use only", position=0, argstr="%s", usedefault=True) + output_file = File("tracks.trk", "file containing tracks", argstr="%s", position=1, usedefault=True) + output_mask = File(desc="output a binary mask file in analyze format", argstr="-om %s") + primary_vector = traits.Enum('v2', 'v3', desc = "which vector to use for fibre tracking: v2 or v3. If not set use v1", argstr="-%s") + +class DTITrackerOutputSpec(TraitedSpec): + track_file = File(exists=True) + mask_file = File(exists=True) + +class DTITracker(CommandLine): + input_spec=DTITrackerInputSpec + output_spec=DTITrackerOutputSpec + + _cmd = 'dti_tracker' + + def _run_interface(self, runtime): + _, _, ext = split_filename(self.inputs.tensor_file) + copyfile(self.inputs.tensor_file, os.path.abspath(self.inputs.input_data_prefix + "_tensor" + ext), copy=False) + + return super(DTITracker, self)._run_interface(runtime) + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['track_file'] = os.path.abspath(self.inputs.output_file) + if isdefined(self.inputs.output_mask) and self.inputs.output_mask: + outputs['mask_file'] = os.path.abspath(self.inputs.output_mask) + + return outputs \ No newline at end of file diff --git a/nipype/interfaces/diffusion_toolkit/odf.py b/nipype/interfaces/diffusion_toolkit/odf.py new file mode 100644 index 0000000000..ce1be469ba --- /dev/null +++ b/nipype/interfaces/diffusion_toolkit/odf.py @@ -0,0 +1,228 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Provides interfaces to various commands provided by diffusion toolkit + + Change directory to provide relative paths for doctests + >>> import os + >>> filepath = os.path.dirname( os.path.realpath( __file__ ) ) + >>> datadir = os.path.realpath(os.path.join(filepath, '../../testing/data')) + >>> os.chdir(datadir) + +""" +import re +from nipype.utils.filemanip import fname_presuffix, split_filename, copyfile +import os +from nipype.utils.misc import isdefined +__docformat__ = 'restructuredtext' + +from nipype.interfaces.base import (TraitedSpec, File, traits, CommandLine, + CommandLineInputSpec) + +class HARDIMatInputSpec(CommandLineInputSpec): + bvecs = File(exists=True, desc = 'b vectors file', + argstr='%s', position=1, mandatory=True) + bvals = File(exists=True,desc = 'b values file', mandatory=True) + out_file = File("recon_mat.dat", desc = 'output matrix file', argstr='%s', usedefault=True, position=2) + order = traits.Int(argsstr='-order %s', desc="""maximum order of spherical harmonics. must be even number. default +is 4""") + odf_file = File(exists=True, argstr='-odf %s', desc="""filename that contains the reconstruction points on a HEMI-sphere. +use the pre-set 181 points by default""") + reference_file = File(exists=True, argstr='-ref %s', desc="""provide a dicom or nifti image as the reference for the program to +figure out the image orientation information. if no such info was +found in the given image header, the next 5 options -info, etc., +will be used if provided. if image orientation info can be found +in the given reference, all other 5 image orientation options will +be IGNORED""") + image_info = File(exists=True, argstr='-info %s', desc="""specify image information file. the image info file is generated +from original dicom image by diff_unpack program and contains image +orientation and other information needed for reconstruction and +tracking. by default will look into the image folder for .info file""") + image_orientation_vectors = traits.List(traits.Float(), minlen=6, maxlen=6, desc="""specify image orientation vectors. if just one argument given, +will treat it as filename and read the orientation vectors from +the file. if 6 arguments are given, will treat them as 6 float +numbers and construct the 1st and 2nd vector and calculate the 3rd +one automatically. +this information will be used to determine image orientation, +as well as to adjust gradient vectors with oblique angle when""", argstr="-iop %f") + oblique_correction = traits.Bool(desc="""when oblique angle(s) applied, some SIEMENS dti protocols do not +adjust gradient accordingly, thus it requires adjustment for correct +diffusion tensor calculation""", argstr="-oc") + +class HARDIMatOutputSpec(TraitedSpec): + out_file = File(exists=True, desc='output matrix file') + + +class HARDIMat(CommandLine): + """Use hardi_mat to calculate a reconstruction matrix from a gradient table + """ + input_spec=HARDIMatInputSpec + output_spec=HARDIMatOutputSpec + + _cmd = 'hardi_mat' + + def _create_gradient_matrix(self, bvecs_file, bvals_file): + _gradient_matrix_file = 'gradient_matrix.txt' + bvals = [val for val in re.split('\s+', open(bvals_file).readline().strip())] + bvecs_f = open(bvecs_file) + bvecs_x = [val for val in re.split('\s+', bvecs_f.readline().strip())] + bvecs_y = [val for val in re.split('\s+', bvecs_f.readline().strip())] + bvecs_z = [val for val in re.split('\s+', bvecs_f.readline().strip())] + bvecs_f.close() + gradient_matrix_f = open(_gradient_matrix_file, 'w') + for i in range(len(bvals)): + if int(bvals[i]) == 0: + continue + gradient_matrix_f.write("%s %s %s\n"%(bvecs_x[i], bvecs_y[i], bvecs_z[i])) + gradient_matrix_f.close() + return _gradient_matrix_file + + def _format_arg(self, name, spec, value): + if name == "bvecs": + new_val = self._create_gradient_matrix(self.inputs.bvecs, self.inputs.bvals) + return super(HARDIMat, self)._format_arg("bvecs", spec, new_val) + return super(HARDIMat, self)._format_arg(name, spec, value) + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['out_file'] = os.path.abspath(self.inputs.out_file) + return outputs + +class ODFReconInputSpec(CommandLineInputSpec): + DWI = File(desc='Input raw data', argstr='%s',exists=True, mandatory=True,position=1) + n_directions = traits.Int(desc='Number of directions', argstr='%s', mandatory=True, position=2) + n_output_directions = traits.Int(desc='Number of output directions', argstr='%s', mandatory=True, position=3) + out_prefix = traits.Str("odf", desc='Output file prefix', argstr='%s', usedefault=True, position=4) + matrix = File(argstr='-mat %s', exists=True, desc="""use given file as reconstruction matrix.""", mandatory=True) + n_b0 = traits.Int(argstr='-b0 %s', desc="""number of b0 scans. by default the program gets this information +from the number of directions and number of volumes in +the raw data. useful when dealing with incomplete raw +data set or only using part of raw data set to reconstruct""", mandatory=True) + output_type = traits.Enum('nii', 'analyze', 'ni1', 'nii.gz', argstr='-ot %s', desc='output file type', usedefault=True) + sharpness = traits.Float(desc="""smooth or sharpen the raw data. factor > 0 is smoothing. +factor < 0 is sharpening. default value is 0 +NOTE: this option applies to DSI study only""", argstr='-s %f') + filter = traits.Bool(desc="""apply a filter (e.g. high pass) to the raw image""", argstr='-f') + subtract_background = traits.Bool(desc="""subtract the background value before reconstruction""", argstr='-bg') + dsi = traits.Bool(desc="""indicates that the data is dsi""", argstr='-dsi') + output_entropy = traits.Bool(desc="""output entropy map""", argstr='-oe') + image_orientation_vectors = traits.List(traits.Float(), minlen=6, maxlen=6, desc="""specify image orientation vectors. if just one argument given, +will treat it as filename and read the orientation vectors from +the file. if 6 arguments are given, will treat them as 6 float +numbers and construct the 1st and 2nd vector and calculate the 3rd +one automatically. +this information will be used to determine image orientation, +as well as to adjust gradient vectors with oblique angle when""", argstr="-iop %f") + oblique_correction = traits.Bool(desc="""when oblique angle(s) applied, some SIEMENS dti protocols do not +adjust gradient accordingly, thus it requires adjustment for correct +diffusion tensor calculation""", argstr="-oc") + + +class ODFReconOutputSpec(TraitedSpec): + B0 = File(exists=True) + DWI = File(exists=True) + max = File(exists=True) + ODF = File(exists=True) + entropy = File() + +class ODFRecon(CommandLine): + """Use odf_recon to generate tensors and other maps + """ + + input_spec=ODFReconInputSpec + output_spec=ODFReconOutputSpec + + _cmd = 'odf_recon' + + def _list_outputs(self): + out_prefix = self.inputs.out_prefix + output_type = self.inputs.output_type + + outputs = self.output_spec().get() + outputs['B0'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_b0.'+ output_type)) + outputs['DWI'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_dwi.'+ output_type)) + outputs['max'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_max.'+ output_type)) + outputs['ODF'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_odf.'+ output_type)) + if isdefined(self.inputs.output_entropy): + outputs['entropy'] = os.path.abspath(fname_presuffix("", prefix=out_prefix, suffix='_entropy.'+ output_type)) + + return outputs + +class ODFTrackerInputSpec(CommandLineInputSpec): + max = File(exists=True, mandatory=True) + ODF = File(exists=True, mandatory=True) + input_data_prefix = traits.Str("odf", desc='recon data prefix', argstr='%s', usedefault=True, position=0) + out_file = File("tracks.trk", desc = 'output track file', argstr='%s', usedefault=True, position=1) + input_output_type = traits.Enum('nii', 'analyze', 'ni1', 'nii.gz', argstr='-it %s', desc='input and output file type', usedefault=True) + runge_kutta2 = traits.Bool(argstr='-rk2', desc="""use 2nd order runge-kutta method for tracking. +default tracking method is non-interpolate streamline""") + step_length = traits.Float(argstr='-l %f', desc="""set step length, in the unit of minimum voxel size. +default value is 0.1.""") + angle_threshold = traits.Float(argstr='-at %f',desc="""set angle threshold. default value is 35 degree for +default tracking method and 25 for rk2""") + random_seed = traits.Int(argstr='-rseed %s', desc="""use random location in a voxel instead of the center of the voxel +to seed. can also define number of seed per voxel. default is 1""") + invert_x = traits.Bool(argstr='-ix', desc='invert x component of the vector') + invert_y = traits.Bool(argstr='-iy', desc='invert y component of the vector') + invert_z = traits.Bool(argstr='-iz', desc='invert z component of the vector') + swap_xy = traits.Bool(argstr='-sxy', desc='swap x and y vectors while tracking') + swap_yz = traits.Bool(argstr='-syz', desc='swap y and z vectors while tracking') + swap_zx = traits.Bool(argstr='-szx', desc='swap x and z vectors while tracking') + disc = traits.Bool(argstr='-disc', desc='use disc tracking') + mask1_file = File(desc="first mask image", mandatory=True, argstr="-m %s", position=2) + mask1_threshold = traits.Float(desc="threshold value for the first mask image, if not given, the program will \ +try automatically find the threshold", position=3) + mask2_file = File(desc="second mask image", argstr="-m2 %s", position=4) + mask2_threshold = traits.Float(desc="threshold value for the second mask image, if not given, the program will \ +try automatically find the threshold", position=5) + limit = traits.Int(argstr='-limit %d', desc="""in some special case, such as heart data, some track may go into +infinite circle and take long time to stop. this option allows +setting a limit for the longest tracking steps (voxels)""") + dsi = traits.Bool(argstr='-dsi', desc=""" specify the input odf data is dsi. because dsi recon uses fixed +pre-calculated matrix, some special orientation patch needs to +be applied to keep dti/dsi/q-ball consistent.""") + image_orientation_vectors = traits.List(traits.Float(), minlen=6, maxlen=6, desc="""specify image orientation vectors. if just one argument given, +will treat it as filename and read the orientation vectors from +the file. if 6 arguments are given, will treat them as 6 float +numbers and construct the 1st and 2nd vector and calculate the 3rd +one automatically. +this information will be used to determine image orientation, +as well as to adjust gradient vectors with oblique angle when""", argstr="-iop %f") + slice_order = traits.Int(argstr='-sorder %d', desc='set the slice order. 1 means normal, -1 means reversed. default value is 1') + voxel_order = traits.Enum('RAS', 'RPS', 'RAI', 'RPI', 'LAI', 'LAS', 'LPS', 'LPI', argstr='-vorder %s', desc="""specify the voxel order in RL/AP/IS (human brain) reference. must be +3 letters with no space in between. +for example, RAS means the voxel row is from L->R, the column +is from P->A and the slice order is from I->S. +by default voxel order is determined by the image orientation +(but NOT guaranteed to be correct because of various standards). +for example, siemens axial image is LPS, coronal image is LIP and +sagittal image is PIL. +this information also is NOT needed for tracking but will be saved +in the track file and is essential for track display to map onto +the right coordinates""") + +class ODFTrackerOutputSpec(TraitedSpec): + track_file = File(exists=True, desc='output track file') + +class ODFTracker(CommandLine): + """Use odf_tracker to generate track file + """ + + input_spec=ODFTrackerInputSpec + output_spec=ODFTrackerOutputSpec + + _cmd = 'odf_tracker' + + def _run_interface(self, runtime): + _, _, ext = split_filename(self.inputs.max) + copyfile(self.inputs.max, os.path.abspath(self.inputs.input_data_prefix + "_max" + ext), copy=False) + + _, _, ext = split_filename(self.inputs.ODF) + copyfile(self.inputs.ODF, os.path.abspath(self.inputs.input_data_prefix + "_odf" + ext), copy=False) + + return super(ODFTracker, self)._run_interface(runtime) + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['track_file'] = os.path.abspath(self.inputs.out_file) + return outputs + diff --git a/nipype/interfaces/diffusion_toolkit/postproc.py b/nipype/interfaces/diffusion_toolkit/postproc.py new file mode 100644 index 0000000000..20633ccf94 --- /dev/null +++ b/nipype/interfaces/diffusion_toolkit/postproc.py @@ -0,0 +1,35 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Provides interfaces to various commands provided by diffusion toolkit + + Change directory to provide relative paths for doctests + >>> import os + >>> filepath = os.path.dirname( os.path.realpath( __file__ ) ) + >>> datadir = os.path.realpath(os.path.join(filepath, '../../testing/data')) + >>> os.chdir(datadir) + +""" +import os +__docformat__ = 'restructuredtext' + +from nipype.interfaces.base import (TraitedSpec, File, traits, CommandLine, + CommandLineInputSpec) + +class SplineFilterInputSpec(CommandLineInputSpec): + track_file = File(exists=True, desc="file containing tracks to be filtered", position=0, argstr="%s", mandatory=True) + step_length = traits.Float(desc="in the unit of minimum voxel size", position=1, argstr="%f", mandatory=True) + output_file = File("spline_tracks.trk", desc="target file for smoothed tracks", position=2, argstr="%s", usedefault=True) + +class SplineFilterOutputSpec(TraitedSpec): + smoothed_track_file = File(exists=True) + +class SplineFilter(CommandLine): + input_spec=SplineFilterInputSpec + output_spec=SplineFilterOutputSpec + + _cmd = "spline_filter" + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['smoothed_track_file'] = os.path.abspath(self.inputs.output_file) + return outputs \ No newline at end of file diff --git a/nipype/interfaces/diffusion_toolkit/preproc.py b/nipype/interfaces/diffusion_toolkit/preproc.py deleted file mode 100644 index 4d9eb080f5..0000000000 --- a/nipype/interfaces/diffusion_toolkit/preproc.py +++ /dev/null @@ -1,105 +0,0 @@ -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -"""Provides interfaces to various commands provided by diffusion toolkit - - Change directory to provide relative paths for doctests - >>> import os - >>> filepath = os.path.dirname( os.path.realpath( __file__ ) ) - >>> datadir = os.path.realpath(os.path.join(filepath, '../../testing/data')) - >>> os.chdir(datadir) - -""" -import re -from nipype.utils.filemanip import fname_presuffix -__docformat__ = 'restructuredtext' - -from nipype.interfaces.base import (TraitedSpec, File, traits, CommandLine, - CommandLineInputSpec) - - -class DTIReconInputSpec(CommandLineInputSpec): - dwi = File(desc='Input diffusion volume', argstr='%s',exists=True, mandatory=True,position=1) - out_prefix = traits.Str("dti", desc='Output file prefix', argstr='%s', usedefault=True,position=2) - output_type = traits.Enum('nii', 'analyze', 'ni1', 'nii.gz', argstr='-ot %s', desc='output file type', usedefault=True) - bvecs = File(exists=True, desc = 'b vectors file', - argstr='-gm %s', mandatory=True) - bvals = File(exists=True,desc = 'b values file', mandatory=True) - n_averages = traits.Int(desc='Number of averages', argstr='-nex %s') - image_orientation_vectors = traits.List(traits.Float(), minlen=6, maxlen=6, desc="""specify image orientation vectors. if just one argument given, -will treat it as filename and read the orientation vectors from -the file. if 6 arguments are given, will treat them as 6 float -numbers and construct the 1st and 2nd vector and calculate the 3rd -one automatically. -this information will be used to determine image orientation, -as well as to adjust gradient vectors with oblique angle when""", argstr="-iop %f") - oblique_correction = traits.Bool(desc="""when oblique angle(s) applied, some SIEMENS dti protocols do not -adjust gradient accordingly, thus it requires adjustment for correct -diffusion tensor calculation""", argstr="-oc") - b0_threshold = traits.Float(desc="""program will use b0 image with the given threshold to mask out high -background of fa/adc maps. by default it will calculate threshold -automatically. but if it failed, you need to set it manually.""", argstr="-b0_th") - - -class DTIReconOutputSpec(TraitedSpec): - ADC = File(exists=True) - B0 = File(exists=True) - L1 = File(exists=True) - L2 = File(exists=True) - L3 = File(exists=True) - exp = File(exists=True) - FA = File(exists=True) - FA_color = File(exists=True) - tensor = File(exists=True) - V1 = File(exists=True) - V2 = File(exists=True) - V3 = File(exists=True) - -class DTIRecon(CommandLine): - """Use dti_recon to generate tensors and other maps - """ - - input_spec=DTIReconInputSpec - output_spec=DTIReconOutputSpec - - _gradient_matrix_file = 'gradient_matrix.txt' - _cmd = 'dti_recon' - - def _format_arg(self, name, spec, value): - if name == "bvecs": - new_val = self._create_gradient_matrix(self.inputs.bvecs, self.inputs.bvals) - return super(DTIRecon, self)._format_arg("bvecs", spec, new_val) - return super(DTIRecon, self)._format_arg(name, spec, value) - - def _create_gradient_matrix(self, bvecs_file, bvals_file): - bvals = [val for val in re.split('\s+', open(bvals_file).readline().strip())] - bvecs_f = open(bvecs_file) - bvecs_x = [val for val in re.split('\s+', bvecs_f.readline().strip())] - bvecs_y = [val for val in re.split('\s+', bvecs_f.readline().strip())] - bvecs_z = [val for val in re.split('\s+', bvecs_f.readline().strip())] - bvecs_f.close() - gradient_matrix_f = open(self._gradient_matrix_file, 'w') - print len(bvals), len(bvecs_x), len(bvecs_y), len(bvecs_z) - for i in range(len(bvals)): - gradient_matrix_f.write("%s, %s, %s, %s\n"%(bvecs_x[i], bvecs_y[i], bvecs_z[i], bvals[i])) - gradient_matrix_f.close() - return self._gradient_matrix_file - - def _list_outputs(self): - out_prefix = self.inputs.out_prefix - output_type = self.inputs.output_type - - outputs = self.output_spec().get() - outputs['ADC'] = fname_presuffix("", prefix=out_prefix, suffix='_adc.'+ output_type) - outputs['B0'] = fname_presuffix("", prefix=out_prefix, suffix='_b0.'+ output_type) - outputs['L1'] = fname_presuffix("", prefix=out_prefix, suffix='_e1.'+ output_type) - outputs['L2'] = fname_presuffix("", prefix=out_prefix, suffix='_e2.'+ output_type) - outputs['L3'] = fname_presuffix("", prefix=out_prefix, suffix='_e3.'+ output_type) - outputs['exp'] = fname_presuffix("", prefix=out_prefix, suffix='_exp.'+ output_type) - outputs['FA'] = fname_presuffix("", prefix=out_prefix, suffix='_fa.'+ output_type) - outputs['FA_color'] = fname_presuffix("", prefix=out_prefix, suffix='_fa_color.'+ output_type) - outputs['tensor'] = fname_presuffix("", prefix=out_prefix, suffix='_tensor.'+ output_type) - outputs['V1'] = fname_presuffix("", prefix=out_prefix, suffix='_v1.'+ output_type) - outputs['V2'] = fname_presuffix("", prefix=out_prefix, suffix='_v2.'+ output_type) - outputs['V3'] = fname_presuffix("", prefix=out_prefix, suffix='_v3.'+ output_type) - - return outputs \ No newline at end of file