diff --git a/CHANGES b/CHANGES index 60987f1481..937483c71e 100644 --- a/CHANGES +++ b/CHANGES @@ -1,6 +1,8 @@ Next release ============ +* ENH: New interfaces in dipy: RESTORE, EstimateResponseSH, CSD and StreamlineTractography + (https://github.com/nipy/nipype/pull/1090) * ENH: Added interfaces of AFNI (https://github.com/nipy/nipype/pull/1360, https://github.com/nipy/nipype/pull/1361) * ENH: Added support for PETPVC (https://github.com/nipy/nipype/pull/1335) diff --git a/circle.yml b/circle.yml index fb53420362..010f99e3eb 100644 --- a/circle.yml +++ b/circle.yml @@ -21,7 +21,7 @@ dependencies: # Set up python environment - pip install --upgrade pip - pip install -e . - - pip install matplotlib sphinx ipython boto coverage + - pip install matplotlib sphinx ipython boto coverage dipy - gem install fakes3 - if [[ ! -d ~/examples/data ]]; then wget "http://tcpdiag.dl.sourceforge.net/project/nipy/nipype/nipype-0.2/nipype-tutorial.tar.bz2" && tar jxvf nipype-tutorial.tar.bz2 && mv nipype-tutorial/* ~/examples/; fi - if [[ ! -d ~/examples/fsl_course_data ]]; then wget -c "http://fsl.fmrib.ox.ac.uk/fslcourse/fdt1.tar.gz" && wget -c "http://fsl.fmrib.ox.ac.uk/fslcourse/fdt2.tar.gz" && wget -c "http://fsl.fmrib.ox.ac.uk/fslcourse/tbss.tar.gz" && mkdir ~/examples/fsl_course_data && tar zxvf fdt1.tar.gz -C ~/examples/fsl_course_data && tar zxvf fdt2.tar.gz -C ~/examples/fsl_course_data && tar zxvf tbss.tar.gz -C ~/examples/fsl_course_data; fi diff --git a/nipype/interfaces/dipy/__init__.py b/nipype/interfaces/dipy/__init__.py index d45591368c..4047e5a408 100644 --- a/nipype/interfaces/dipy/__init__.py +++ b/nipype/interfaces/dipy/__init__.py @@ -1,4 +1,5 @@ -from .tracks import TrackDensityMap +from .tracks import StreamlineTractography, TrackDensityMap from .tensors import TensorMode, DTI from .preprocess import Resample, Denoise +from .reconstruction import RESTORE, EstimateResponseSH, CSD from .simulate import SimulateMultiTensor diff --git a/nipype/interfaces/dipy/base.py b/nipype/interfaces/dipy/base.py new file mode 100644 index 0000000000..d760172710 --- /dev/null +++ b/nipype/interfaces/dipy/base.py @@ -0,0 +1,87 @@ +# -*- coding: utf-8 -*- +""" Base interfaces for dipy """ +import os.path as op +import numpy as np +from nipype.interfaces.base import (traits, File, isdefined, + BaseInterface, BaseInterfaceInputSpec) +from ... import logging + +IFLOGGER = logging.getLogger('interface') + +HAVE_DIPY = True +try: + import dipy +except ImportError: + HAVE_DIPY = False + + +def no_dipy(): + """ Check if dipy is available """ + global HAVE_DIPY + return not HAVE_DIPY + + +def dipy_version(): + """ Check dipy version """ + if no_dipy(): + return None + + return dipy.__version__ + + +class DipyBaseInterface(BaseInterface): + + """ + A base interface for py:mod:`dipy` computations + """ + def __init__(self, **inputs): + if no_dipy(): + IFLOGGER.error('dipy was not found') + # raise ImportError('dipy was not found') + super(DipyBaseInterface, self).__init__(**inputs) + + +class DipyBaseInterfaceInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc=('input diffusion data')) + in_bval = File(exists=True, mandatory=True, desc=('input b-values table')) + in_bvec = File(exists=True, mandatory=True, desc=('input b-vectors table')) + b0_thres = traits.Int(700, usedefault=True, desc=('b0 threshold')) + out_prefix = traits.Str(desc=('output prefix for file names')) + + +class DipyDiffusionInterface(DipyBaseInterface): + + """ + A base interface for py:mod:`dipy` computations + """ + input_spec = DipyBaseInterfaceInputSpec + + def _get_gradient_table(self): + bval = np.loadtxt(self.inputs.in_bval) + bvec = np.loadtxt(self.inputs.in_bvec).T + try: + from dipy.data import GradientTable + gtab = GradientTable(bvec) + gtab.bvals = bval + except NameError: + from dipy.core.gradients import gradient_table + gtab = gradient_table(bval, bvec) + + gtab.b0_threshold = self.inputs.b0_thres + return gtab + + def _gen_filename(self, name, ext=None): + fname, fext = op.splitext(op.basename(self.inputs.in_file)) + if fext == '.gz': + fname, fext2 = op.splitext(fname) + fext = fext2 + fext + + if not isdefined(self.inputs.out_prefix): + out_prefix = op.abspath(fname) + else: + out_prefix = self.inputs.out_prefix + + if ext is None: + ext = fext + + return out_prefix + '_' + name + ext diff --git a/nipype/interfaces/dipy/preprocess.py b/nipype/interfaces/dipy/preprocess.py index 5ecc49b957..143f239e6c 100644 --- a/nipype/interfaces/dipy/preprocess.py +++ b/nipype/interfaces/dipy/preprocess.py @@ -1,28 +1,20 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -"""Change directory to provide relative paths for doctests +""" +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.path as op -import warnings - import nibabel as nb import numpy as np -from ..base import (traits, TraitedSpec, BaseInterface, File, isdefined) -from ...utils.filemanip import split_filename -from ...utils.misc import package_check +from ..base import (traits, TraitedSpec, File, isdefined) +from .base import DipyBaseInterface from ... import logging -iflogger = logging.getLogger('interface') - -have_dipy = True -try: - package_check('dipy', version='0.6.0') -except Exception as e: - have_dipy = False +IFLOGGER = logging.getLogger('interface') class ResampleInputSpec(TraitedSpec): @@ -33,15 +25,17 @@ class ResampleInputSpec(TraitedSpec): ' is set, then isotropic regridding will ' 'be performed, with spacing equal to the ' 'smallest current zoom.')) - interp = traits.Int(1, mandatory=True, usedefault=True, desc=('order of ' - 'the interpolator (0 = nearest, 1 = linear, etc.')) + interp = traits.Int( + 1, mandatory=True, usedefault=True, + desc=('order of the interpolator (0 = nearest, 1 = linear, etc.')) class ResampleOutputSpec(TraitedSpec): out_file = File(exists=True) -class Resample(BaseInterface): +class Resample(DipyBaseInterface): + """ An interface to reslicing diffusion datasets. See @@ -69,7 +63,7 @@ def _run_interface(self, runtime): resample_proxy(self.inputs.in_file, order=order, new_zooms=vox_size, out_file=out_file) - iflogger.info('Resliced image saved as {i}'.format(i=out_file)) + IFLOGGER.info('Resliced image saved as {i}'.format(i=out_file)) return runtime def _list_outputs(self): @@ -92,17 +86,21 @@ class DenoiseInputSpec(TraitedSpec): noise_model = traits.Enum('rician', 'gaussian', mandatory=True, usedefault=True, desc=('noise distribution model')) + signal_mask = File(desc=('mask in which the mean signal ' + 'will be computed'), exists=True) noise_mask = File(desc=('mask in which the standard deviation of noise ' 'will be computed'), exists=True) patch_radius = traits.Int(1, desc='patch radius') block_radius = traits.Int(5, desc='block_radius') + snr = traits.Float(desc='manually set an SNR') class DenoiseOutputSpec(TraitedSpec): out_file = File(exists=True) -class Denoise(BaseInterface): +class Denoise(DipyBaseInterface): + """ An interface to denoising diffusion datasets [Coupe2008]_. See @@ -140,16 +138,24 @@ def _run_interface(self, runtime): if isdefined(self.inputs.block_radius): settings['block_radius'] = self.inputs.block_radius + snr = None + if isdefined(self.inputs.snr): + snr = self.inputs.snr + + signal_mask = None + if isdefined(self.inputs.signal_mask): + signal_mask = nb.load(self.inputs.signal_mask).get_data() noise_mask = None - if isdefined(self.inputs.in_mask): + if isdefined(self.inputs.noise_mask): noise_mask = nb.load(self.inputs.noise_mask).get_data() - _, s = nlmeans_proxy(self.inputs.in_file, - settings, - noise_mask=noise_mask, + _, s = nlmeans_proxy(self.inputs.in_file, settings, + snr=snr, + smask=signal_mask, + nmask=noise_mask, out_file=out_file) - iflogger.info(('Denoised image saved as {i}, estimated ' - 'sigma={s}').format(i=out_file, s=s)) + IFLOGGER.info(('Denoised image saved as {i}, estimated ' + 'SNR={s}').format(i=out_file, s=str(s))) return runtime def _list_outputs(self): @@ -169,7 +175,7 @@ def resample_proxy(in_file, order=3, new_zooms=None, out_file=None): """ Performs regridding of an image to set isotropic voxel sizes using dipy. """ - from dipy.align.aniso2iso import resample + from dipy.align.reslice import reslice if out_file is None: fname, fext = op.splitext(op.basename(in_file)) @@ -191,7 +197,7 @@ def resample_proxy(in_file, order=3, new_zooms=None, out_file=None): if np.all(im_zooms == new_zooms): return in_file - data2, affine2 = resample(data, affine, im_zooms, new_zooms, order=order) + data2, affine2 = reslice(data, affine, im_zooms, new_zooms, order=order) tmp_zooms = np.array(hdr.get_zooms()) tmp_zooms[:3] = new_zooms[0] hdr.set_zooms(tuple(tmp_zooms)) @@ -203,12 +209,16 @@ def resample_proxy(in_file, order=3, new_zooms=None, out_file=None): def nlmeans_proxy(in_file, settings, - noise_mask=None, out_file=None): + snr=None, + smask=None, + nmask=None, + out_file=None): """ Uses non-local means to denoise 4D datasets """ - package_check('dipy', version='0.8.0.dev') from dipy.denoise.nlmeans import nlmeans + from scipy.ndimage.morphology import binary_erosion + from scipy import ndimage if out_file is None: fname, fext = op.splitext(op.basename(in_file)) @@ -222,13 +232,69 @@ def nlmeans_proxy(in_file, settings, data = img.get_data() aff = img.affine - nmask = data[..., 0] > 80 - if noise_mask is not None: - nmask = noise_mask > 0 - - sigma = np.std(data[nmask == 1]) - den = nlmeans(data, sigma, **settings) + if data.ndim < 4: + data = data[..., np.newaxis] + + data = np.nan_to_num(data) + + if data.max() < 1.0e-4: + raise RuntimeError('There is no signal in the image') + + df = 1.0 + if data.max() < 1000.0: + df = 1000. / data.max() + data *= df + + b0 = data[..., 0] + + if smask is None: + smask = np.zeros_like(b0) + smask[b0 > np.percentile(b0, 85.)] = 1 + + smask = binary_erosion( + smask.astype(np.uint8), iterations=2).astype(np.uint8) + + if nmask is None: + nmask = np.ones_like(b0, dtype=np.uint8) + bmask = settings['mask'] + if bmask is None: + bmask = np.zeros_like(b0) + bmask[b0 > np.percentile(b0[b0 > 0], 10)] = 1 + label_im, nb_labels = ndimage.label(bmask) + sizes = ndimage.sum(bmask, label_im, range(nb_labels + 1)) + maxidx = np.argmax(sizes) + bmask = np.zeros_like(b0, dtype=np.uint8) + bmask[label_im == maxidx] = 1 + nmask[bmask > 0] = 0 + else: + nmask = np.squeeze(nmask) + nmask[nmask > 0.0] = 1 + nmask[nmask < 1] = 0 + nmask = nmask.astype(bool) + + nmask = binary_erosion(nmask, iterations=1).astype(np.uint8) + + den = np.zeros_like(data) + + est_snr = True + if snr is not None: + snr = [snr] * data.shape[-1] + est_snr = False + else: + snr = [] + + for i in range(data.shape[-1]): + d = data[..., i] + if est_snr: + s = np.mean(d[smask > 0]) + n = np.std(d[nmask > 0]) + snr.append(s / n) + + den[..., i] = nlmeans(d, snr[i], **settings) + + den = np.squeeze(den) + den /= df nb.Nifti1Image(den.astype(hdr.get_data_dtype()), aff, hdr).to_filename(out_file) - return out_file, sigma + return out_file, snr diff --git a/nipype/interfaces/dipy/reconstruction.py b/nipype/interfaces/dipy/reconstruction.py new file mode 100644 index 0000000000..6a01eacbf0 --- /dev/null +++ b/nipype/interfaces/dipy/reconstruction.py @@ -0,0 +1,368 @@ +# -*- coding: utf-8 -*- +""" +Interfaces to the reconstruction algorithms in dipy + +""" +import os.path as op + +import numpy as np +import nibabel as nb + +from nipype.interfaces.base import TraitedSpec, File, traits, isdefined +from .base import DipyDiffusionInterface, DipyBaseInterfaceInputSpec + +from nipype import logging +IFLOGGER = logging.getLogger('interface') + + +class RESTOREInputSpec(DipyBaseInterfaceInputSpec): + in_mask = File(exists=True, desc=('input mask in which compute tensors')) + noise_mask = File( + exists=True, desc=('input mask in which compute noise variance')) + + +class RESTOREOutputSpec(TraitedSpec): + fa = File(desc='output fractional anisotropy (FA) map computed from ' + 'the fitted DTI') + md = File(desc='output mean diffusivity (MD) map computed from the ' + 'fitted DTI') + rd = File(desc='output radial diffusivity (RD) map computed from ' + 'the fitted DTI') + mode = File(desc=('output mode (MO) map computed from the fitted DTI')) + trace = File(desc=('output the tensor trace map computed from the ' + 'fitted DTI')) + evals = File(desc=('output the eigenvalues of the fitted DTI')) + evecs = File(desc=('output the eigenvectors of the fitted DTI')) + + +class RESTORE(DipyDiffusionInterface): + + """ + Uses RESTORE [Chang2005]_ to perform DTI fitting with outlier detection. + The interface uses :py:mod:`dipy`, as explained in `dipy's documentation`_. + + .. [Chang2005] Chang, LC, Jones, DK and Pierpaoli, C. RESTORE: robust \ + estimation of tensors by outlier rejection. MRM, 53:1088-95, (2005). + + .. _dipy's documentation: \ + http://nipy.org/dipy/examples_built/restore_dti.html + + + Example + ------- + + >>> from nipype.interfaces import dipy as ndp + >>> dti = ndp.RESTORE() + >>> dti.inputs.in_file = '4d_dwi.nii' + >>> dti.inputs.in_bval = 'bvals' + >>> dti.inputs.in_bvec = 'bvecs' + >>> res = dti.run() # doctest: +SKIP + + + """ + input_spec = RESTOREInputSpec + output_spec = RESTOREOutputSpec + + def _run_interface(self, runtime): + from scipy.special import gamma + from dipy.reconst.dti import TensorModel + import gc + + img = nb.load(self.inputs.in_file) + hdr = img.get_header().copy() + affine = img.get_affine() + data = img.get_data() + gtab = self._get_gradient_table() + + if isdefined(self.inputs.in_mask): + msk = nb.load(self.inputs.in_mask).get_data().astype(np.uint8) + else: + msk = np.ones(data.shape[:3], dtype=np.uint8) + + try_b0 = True + if isdefined(self.inputs.noise_mask): + noise_msk = nb.load(self.inputs.noise_mask).get_data().reshape(-1) + noise_msk[noise_msk > 0.5] = 1 + noise_msk[noise_msk < 1.0] = 0 + noise_msk = noise_msk.astype(np.uint8) + try_b0 = False + elif np.all(data[msk == 0, 0] == 0): + IFLOGGER.info('Input data are masked.') + noise_msk = msk.reshape(-1).astype(np.uint8) + else: + noise_msk = (1 - msk).reshape(-1).astype(np.uint8) + + nb0 = np.sum(gtab.b0s_mask) + dsample = data.reshape(-1, data.shape[-1]) + + if try_b0 and (nb0 > 1): + noise_data = dsample.take(np.where(gtab.b0s_mask), + axis=-1)[noise_msk == 0, ...] + n = nb0 + else: + nodiff = np.where(~gtab.b0s_mask) + nodiffidx = nodiff[0].tolist() + n = 20 if len(nodiffidx) >= 20 else len(nodiffidx) + idxs = np.random.choice(nodiffidx, size=n, replace=False) + noise_data = dsample.take(idxs, axis=-1)[noise_msk == 1, ...] + + # Estimate sigma required by RESTORE + mean_std = np.median(noise_data.std(-1)) + try: + bias = (1. - np.sqrt(2. / (n - 1)) * + (gamma(n / 2.) / gamma((n - 1) / 2.))) + except: + bias = .0 + pass + + sigma = mean_std * (1 + bias) + + if sigma == 0: + IFLOGGER.warn( + ('Noise std is 0.0, looks like data was masked and noise' + ' cannot be estimated correctly. Using default tensor ' + 'model instead of RESTORE.')) + dti = TensorModel(gtab) + else: + IFLOGGER.info(('Performing RESTORE with noise std=%.4f.') % sigma) + dti = TensorModel(gtab, fit_method='RESTORE', sigma=sigma) + + try: + fit_restore = dti.fit(data, msk) + except TypeError: + dti = TensorModel(gtab) + fit_restore = dti.fit(data, msk) + + hdr.set_data_dtype(np.float32) + hdr['data_type'] = 16 + + for k in self._outputs().get(): + scalar = getattr(fit_restore, k) + hdr.set_data_shape(np.shape(scalar)) + nb.Nifti1Image(scalar.astype(np.float32), + affine, hdr).to_filename(self._gen_filename(k)) + + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + for k in outputs.keys(): + outputs[k] = self._gen_filename(k) + return outputs + + +class EstimateResponseSHInputSpec(DipyBaseInterfaceInputSpec): + in_evals = File( + exists=True, mandatory=True, desc=('input eigenvalues file')) + in_mask = File( + exists=True, desc=('input mask in which we find single fibers')) + fa_thresh = traits.Float( + 0.7, usedefault=True, desc=('FA threshold')) + roi_radius = traits.Int( + 10, usedefault=True, desc=('ROI radius to be used in auto_response')) + auto = traits.Bool( + True, usedefault=True, xor=['recursive'], + desc='use the auto_response estimator from dipy') + recursive = traits.Bool( + False, usedefault=True, xor=['auto'], + desc='use the recursive response estimator from dipy') + response = File( + 'response.txt', usedefault=True, desc=('the output response file')) + out_mask = File('wm_mask.nii.gz', usedefault=True, desc='computed wm mask') + + +class EstimateResponseSHOutputSpec(TraitedSpec): + response = File(exists=True, desc=('the response file')) + out_mask = File(exists=True, desc=('output wm mask')) + + +class EstimateResponseSH(DipyDiffusionInterface): + + """ + Uses dipy to compute the single fiber response to be used in spherical + deconvolution methods, in a similar way to MRTrix's command + ``estimate_response``. + + + Example + ------- + + >>> from nipype.interfaces import dipy as ndp + >>> dti = ndp.EstimateResponseSH() + >>> dti.inputs.in_file = '4d_dwi.nii' + >>> dti.inputs.in_bval = 'bvals' + >>> dti.inputs.in_bvec = 'bvecs' + >>> dti.inputs.in_evals = 'dwi_evals.nii' + >>> res = dti.run() # doctest: +SKIP + + + """ + input_spec = EstimateResponseSHInputSpec + output_spec = EstimateResponseSHOutputSpec + + def _run_interface(self, runtime): + from dipy.core.gradients import GradientTable + from dipy.reconst.dti import fractional_anisotropy, mean_diffusivity + from dipy.reconst.csdeconv import recursive_response, auto_response + + img = nb.load(self.inputs.in_file) + affine = img.get_affine() + + if isdefined(self.inputs.in_mask): + msk = nb.load(self.inputs.in_mask).get_data() + msk[msk > 0] = 1 + msk[msk < 0] = 0 + else: + msk = np.ones(imref.get_shape()) + + data = img.get_data().astype(np.float32) + gtab = self._get_gradient_table() + + evals = np.nan_to_num(nb.load(self.inputs.in_evals).get_data()) + FA = np.nan_to_num(fractional_anisotropy(evals)) * msk + indices = np.where(FA > self.inputs.fa_thresh) + S0s = data[indices][:, np.nonzero(gtab.b0s_mask)[0]] + S0 = np.mean(S0s) + + if self.inputs.auto: + response, ratio = auto_response(gtab, data, + roi_radius=self.inputs.roi_radius, + fa_thr=self.inputs.fa_thresh) + response = response[0].tolist() + [S0] + elif self.inputs.recursive: + MD = np.nan_to_num(mean_diffusivity(evals)) * msk + indices = np.logical_or( + FA >= 0.4, (np.logical_and(FA >= 0.15, MD >= 0.0011))) + data = nb.load(self.inputs.in_file).get_data() + response = recursive_response(gtab, data, mask=indices, sh_order=8, + peak_thr=0.01, init_fa=0.08, + init_trace=0.0021, iter=8, + convergence=0.001, + parallel=True) + ratio = abs(response[1] / response[0]) + else: + lambdas = evals[indices] + l01 = np.sort(np.mean(lambdas, axis=0)) + + response = np.array([l01[-1], l01[-2], l01[-2], S0]) + ratio = abs(response[1] / response[0]) + + if ratio > 0.25: + IFLOGGER.warn(('Estimated response is not prolate enough. ' + 'Ratio=%0.3f.') % ratio) + elif ratio < 1.e-5 or np.any(np.isnan(response)): + response = np.array([1.8e-3, 3.6e-4, 3.6e-4, S0]) + IFLOGGER.warn( + ('Estimated response is not valid, using a default one')) + else: + IFLOGGER.info(('Estimated response: %s') % str(response[:3])) + + np.savetxt(op.abspath(self.inputs.response), response) + + wm_mask = np.zeros_like(FA) + wm_mask[indices] = 1 + nb.Nifti1Image( + wm_mask.astype(np.uint8), affine, + None).to_filename(op.abspath(self.inputs.out_mask)) + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + outputs['response'] = op.abspath(self.inputs.response) + outputs['out_mask'] = op.abspath(self.inputs.out_mask) + return outputs + + +class CSDInputSpec(DipyBaseInterfaceInputSpec): + in_mask = File(exists=True, desc=('input mask in which compute tensors')) + response = File(exists=True, desc=('single fiber estimated response')) + sh_order = traits.Int(8, exists=True, usedefault=True, + desc=('maximal shperical harmonics order')) + save_fods = traits.Bool(True, exists=True, usedefault=True, + desc=('save fODFs in file')) + out_fods = File(desc=('fODFs output file name')) + + +class CSDOutputSpec(TraitedSpec): + model = File(desc='Python pickled object of the CSD model fitted.') + out_fods = File(desc=('fODFs output file name')) + + +class CSD(DipyDiffusionInterface): + + """ + Uses CSD [Tournier2007]_ to generate the fODF of DWIs. The interface uses + :py:mod:`dipy`, as explained in `dipy's CSD example + `_. + + .. [Tournier2007] Tournier, J.D., et al. NeuroImage 2007. + Robust determination of the fibre orientation distribution in diffusion + MRI: Non-negativity constrained super-resolved spherical deconvolution + + + Example + ------- + + >>> from nipype.interfaces import dipy as ndp + >>> csd = ndp.CSD() + >>> csd.inputs.in_file = '4d_dwi.nii' + >>> csd.inputs.in_bval = 'bvals' + >>> csd.inputs.in_bvec = 'bvecs' + >>> res = csd.run() # doctest: +SKIP + """ + input_spec = CSDInputSpec + output_spec = CSDOutputSpec + + def _run_interface(self, runtime): + from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel + from dipy.data import get_sphere + # import marshal as pickle + import cPickle as pickle + import gzip + + img = nb.load(self.inputs.in_file) + imref = nb.four_to_three(img)[0] + affine = img.get_affine() + + if isdefined(self.inputs.in_mask): + msk = nb.load(self.inputs.in_mask).get_data() + else: + msk = np.ones(imref.get_shape()) + + data = img.get_data().astype(np.float32) + hdr = imref.get_header().copy() + + gtab = self._get_gradient_table() + resp_file = np.loadtxt(self.inputs.response) + + response = (np.array(resp_file[0:3]), resp_file[-1]) + ratio = response[0][1] / response[0][0] + + if abs(ratio - 0.2) > 0.1: + IFLOGGER.warn(('Estimated response is not prolate enough. ' + 'Ratio=%0.3f.') % ratio) + + csd_model = ConstrainedSphericalDeconvModel( + gtab, response, sh_order=self.inputs.sh_order) + + IFLOGGER.info('Fitting CSD model') + csd_fit = csd_model.fit(data, msk) + + f = gzip.open(self._gen_filename('csdmodel', ext='.pklz'), 'wb') + pickle.dump(csd_model, f, -1) + f.close() + + if self.inputs.save_fods: + sphere = get_sphere('symmetric724') + fods = csd_fit.odf(sphere) + nb.Nifti1Image(fods.astype(np.float32), img.get_affine(), + None).to_filename(self._gen_filename('fods')) + + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + outputs['model'] = self._gen_filename('csdmodel', ext='.pklz') + if self.inputs.save_fods: + outputs['out_fods'] = self._gen_filename('fods') + return outputs diff --git a/nipype/interfaces/dipy/setup.py b/nipype/interfaces/dipy/setup.py new file mode 100644 index 0000000000..6a7b5a30a4 --- /dev/null +++ b/nipype/interfaces/dipy/setup.py @@ -0,0 +1,15 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: + + +def configuration(parent_package='', top_path=None): + from numpy.distutils.misc_util import Configuration + + config = Configuration('dipy', parent_package, top_path) + + # config.add_data_dir('tests') + return config + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(**configuration(top_path='').todict()) diff --git a/nipype/interfaces/dipy/simulate.py b/nipype/interfaces/dipy/simulate.py index 939d7d15a4..f4eb174a22 100644 --- a/nipype/interfaces/dipy/simulate.py +++ b/nipype/interfaces/dipy/simulate.py @@ -7,29 +7,18 @@ """ from __future__ import division +from multiprocessing import (Pool, cpu_count) +import os.path as op from builtins import range -import os.path as op -from multiprocessing import (Pool, cpu_count) import nibabel as nb -from ..base import (traits, TraitedSpec, BaseInterface, BaseInterfaceInputSpec, +from ..base import (traits, TraitedSpec, BaseInterfaceInputSpec, File, InputMultiPath, isdefined) -from ...utils.misc import package_check +from .base import DipyBaseInterface from ... import logging -iflogger = logging.getLogger('interface') - -have_dipy = True -try: - package_check('dipy', version='0.8.0') -except Exception as e: - have_dipy = False -else: - import numpy as np - from dipy.sims.voxel import (multi_tensor, add_noise, - all_tensor_evecs) - from dipy.core.gradients import gradient_table +IFLOGGER = logging.getLogger('interface') class SimulateMultiTensorInputSpec(BaseInterfaceInputSpec): @@ -77,7 +66,7 @@ class SimulateMultiTensorOutputSpec(TraitedSpec): out_bval = File(exists=True, desc='simulated b values') -class SimulateMultiTensor(BaseInterface): +class SimulateMultiTensor(DipyBaseInterface): """ Interface to MultiTensor model simulator in dipy @@ -101,6 +90,8 @@ class SimulateMultiTensor(BaseInterface): output_spec = SimulateMultiTensorOutputSpec def _run_interface(self, runtime): + from dipy.core.gradients import gradient_table + # Gradient table if isdefined(self.inputs.in_bval) and isdefined(self.inputs.in_bvec): # Load the gradient strengths and directions @@ -237,7 +228,7 @@ def _run_interface(self, runtime): pool = Pool(processes=n_proc) # Simulate sticks using dipy - iflogger.info(('Starting simulation of %d voxels, %d diffusion' + IFLOGGER.info(('Starting simulation of %d voxels, %d diffusion' ' directions.') % (len(args), ndirs)) result = np.array(pool.map(_compute_voxel, args)) if np.shape(result)[1] != ndirs: @@ -280,6 +271,7 @@ def _compute_voxel(args): .. [Pierpaoli1996] Pierpaoli et al., Diffusion tensor MR imaging of the human brain, Radiology 201:637-648. 1996. """ + from dipy.sims.voxel import multi_tensor ffs = args['fractions'] gtab = args['gradients'] @@ -297,7 +289,7 @@ def _compute_voxel(args): angles=args['sticks'], fractions=ffs, snr=snr) except Exception as e: pass - # iflogger.warn('Exception simulating dwi signal: %s' % e) + # IFLOGGER.warn('Exception simulating dwi signal: %s' % e) return signal.tolist() diff --git a/nipype/interfaces/dipy/tensors.py b/nipype/interfaces/dipy/tensors.py index edafea1306..c7d02c681c 100644 --- a/nipype/interfaces/dipy/tensors.py +++ b/nipype/interfaces/dipy/tensors.py @@ -5,82 +5,25 @@ >>> datadir = os.path.realpath(os.path.join(filepath, '../../testing/data')) >>> os.chdir(datadir) """ -import os.path as op -import warnings - import nibabel as nb -import numpy as np - -from ..base import (TraitedSpec, BaseInterface, File) -from ...utils.filemanip import split_filename -from ...utils.misc import package_check -from ... import logging -iflogger = logging.getLogger('interface') - -have_dipy = True -try: - package_check('dipy', version='0.6.0') -except Exception as e: - have_dipy = False -else: - import dipy.reconst.dti as dti - from dipy.core.gradients import gradient_table - from dipy.io.utils import nifti1_symmat - - -def tensor_fitting(data, bvals, bvecs, mask_file=None): - """ - Use dipy to fit DTI - - Parameters - ---------- - in_file : str - Full path to a DWI data file. - bvals : str - Full path to a file containing gradient magnitude information (b-values). - bvecs : str - Full path to a file containing gradient direction information (b-vectors). - mask_file : str, optional - Full path to a file containing a binary mask. Defaults to use the entire volume. - - Returns - ------- - TensorFit object, affine - """ - img = nb.load(in_file) - data = img.get_data() - affine = img.affine - if mask_file is not None: - mask = nb.load(self.inputs.mask_file).get_data() - else: - mask = None - # Load information about the gradients: - gtab = grad.gradient_table(self.inputs.bvals, self.inputs.bvecs) +from ..base import TraitedSpec, File, isdefined +from .base import DipyDiffusionInterface, DipyBaseInterfaceInputSpec - # Fit it - tenmodel = dti.TensorModel(gtab) - return tenmodel.fit(data, mask), affine +from ... import logging +IFLOGGER = logging.getLogger('interface') -class DTIInputSpec(TraitedSpec): - in_file = File(exists=True, mandatory=True, - desc='The input 4D diffusion-weighted image file') - bvecs = File(exists=True, mandatory=True, - desc='The input b-vector text file') - bvals = File(exists=True, mandatory=True, - desc='The input b-value text file') +class DTIInputSpec(DipyBaseInterfaceInputSpec): mask_file = File(exists=True, desc='An optional white matter mask') - out_filename = File( - genfile=True, desc='The output filename for the DTI parameters image') class DTIOutputSpec(TraitedSpec): out_file = File(exists=True) -class DTI(BaseInterface): +class DTI(DipyDiffusionInterface): """ Calculates the diffusion tensor model parameters @@ -90,59 +33,52 @@ class DTI(BaseInterface): >>> import nipype.interfaces.dipy as dipy >>> dti = dipy.DTI() >>> dti.inputs.in_file = 'diffusion.nii' - >>> dti.inputs.bvecs = 'bvecs' - >>> dti.inputs.bvals = 'bvals' + >>> dti.inputs.in_bvec = 'bvecs' + >>> dti.inputs.in_bval = 'bvals' >>> dti.run() # doctest: +SKIP """ input_spec = DTIInputSpec output_spec = DTIOutputSpec def _run_interface(self, runtime): - ten_fit, affine = tensor_fitting(self.inputs.in_file, - self.inputs.bvals, - self.inputs.bvecs, - self.inputs.mask_file) + from dipy.reconst import dti + from dipy.io.utils import nifti1_symmat + gtab = self._get_gradient_table() + + img = nb.load(self.inputs.in_file) + data = img.get_data() + affine = img.affine + mask = None + if isdefined(self.inputs.mask_file): + mask = nb.load(self.inputs.mask_file).get_data() + + # Fit it + tenmodel = dti.TensorModel(gtab) + ten_fit = tenmodel.fit(data, mask) lower_triangular = tenfit.lower_triangular() img = nifti1_symmat(lower_triangular, affine) - out_file = op.abspath(self._gen_outfilename()) + out_file = self._gen_filename('dti') nb.save(img, out_file) - iflogger.info('DTI parameters image saved as {i}'.format(i=out_file)) + IFLOGGER.info('DTI parameters image saved as {i}'.format(i=out_file)) return runtime def _list_outputs(self): outputs = self._outputs().get() - outputs['out_file'] = op.abspath(self._gen_outfilename()) + outputs['out_file'] = self._gen_filename('dti') return outputs - def _gen_filename(self, name): - if name is 'out_filename': - return self._gen_outfilename() - else: - return None - - def _gen_outfilename(self): - _, name, _ = split_filename(self.inputs.in_file) - return name + '_dti.nii' - -class TensorModeInputSpec(TraitedSpec): - in_file = File(exists=True, mandatory=True, - desc='The input 4D diffusion-weighted image file') - bvecs = File(exists=True, mandatory=True, - desc='The input b-vector text file') - bvals = File(exists=True, mandatory=True, - desc='The input b-value text file') +class TensorModeInputSpec(DipyBaseInterfaceInputSpec): mask_file = File(exists=True, desc='An optional white matter mask') - out_filename = File( - genfile=True, desc='The output filename for the Tensor mode image') class TensorModeOutputSpec(TraitedSpec): out_file = File(exists=True) -class TensorMode(BaseInterface): +class TensorMode(DipyDiffusionInterface): + """ Creates a map of the mode of the diffusion tensors given a set of diffusion-weighted images, as well as their associated b-values and @@ -160,35 +96,43 @@ class TensorMode(BaseInterface): >>> import nipype.interfaces.dipy as dipy >>> mode = dipy.TensorMode() >>> mode.inputs.in_file = 'diffusion.nii' - >>> mode.inputs.bvecs = 'bvecs' - >>> mode.inputs.bvals = 'bvals' + >>> mode.inputs.in_bvec = 'bvecs' + >>> mode.inputs.in_bval = 'bvals' >>> mode.run() # doctest: +SKIP """ input_spec = TensorModeInputSpec output_spec = TensorModeOutputSpec def _run_interface(self, runtime): - ten_fit = tensor_fitting(self.inputs.in_file, self.inputs.bvals, self.inputs.bvecs, - self.inputs.mask_file) + from dipy.reconst import dti + + # Load the 4D image files + img = nb.load(self.inputs.in_file) + data = img.get_data() + affine = img.get_affine() + + # Load the gradient strengths and directions + gtab = self._get_gradient_table() + + # Mask the data so that tensors are not fit for + # unnecessary voxels + mask = data[..., 0] > 50 + + # Fit the tensors to the data + tenmodel = dti.TensorModel(gtab) + tenfit = tenmodel.fit(data, mask) + + # Calculate the mode of each voxel's tensor + mode_data = tenfit.mode # Write as a 3D Nifti image with the original affine - img = nb.Nifti1Image(tenfit.mode, affine) - out_file = op.abspath(self._gen_outfilename()) + img = nb.Nifti1Image(mode_data, affine) + out_file = self._gen_filename('mode') nb.save(img, out_file) - iflogger.info('Tensor mode image saved as {i}'.format(i=out_file)) + IFLOGGER.info('Tensor mode image saved as {i}'.format(i=out_file)) return runtime def _list_outputs(self): outputs = self._outputs().get() - outputs['out_file'] = op.abspath(self._gen_outfilename()) + outputs['out_file'] = self._gen_filename('mode') return outputs - - def _gen_filename(self, name): - if name is 'out_filename': - return self._gen_outfilename() - else: - return None - - def _gen_outfilename(self): - _, name, _ = split_filename(self.inputs.in_file) - return name + '_mode.nii' diff --git a/nipype/interfaces/dipy/tests/test_auto_CSD.py b/nipype/interfaces/dipy/tests/test_auto_CSD.py new file mode 100644 index 0000000000..d99aac4527 --- /dev/null +++ b/nipype/interfaces/dipy/tests/test_auto_CSD.py @@ -0,0 +1,44 @@ +# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT +from ....testing import assert_equal +from ..reconstruction import CSD + + +def test_CSD_inputs(): + input_map = dict(b0_thres=dict(usedefault=True, + ), + ignore_exception=dict(nohash=True, + usedefault=True, + ), + in_bval=dict(mandatory=True, + ), + in_bvec=dict(mandatory=True, + ), + in_file=dict(mandatory=True, + ), + in_mask=dict(), + out_fods=dict(), + out_prefix=dict(), + response=dict(), + save_fods=dict(exists=True, + usedefault=True, + ), + sh_order=dict(exists=True, + usedefault=True, + ), + ) + inputs = CSD.input_spec() + + for key, metadata in list(input_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(inputs.traits()[key], metakey), value + + +def test_CSD_outputs(): + output_map = dict(model=dict(), + out_fods=dict(), + ) + outputs = CSD.output_spec() + + for key, metadata in list(output_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(outputs.traits()[key], metakey), value diff --git a/nipype/interfaces/dipy/tests/test_auto_DTI.py b/nipype/interfaces/dipy/tests/test_auto_DTI.py index 0abba84540..eb35d6a9e2 100644 --- a/nipype/interfaces/dipy/tests/test_auto_DTI.py +++ b/nipype/interfaces/dipy/tests/test_auto_DTI.py @@ -4,15 +4,19 @@ def test_DTI_inputs(): - input_map = dict(bvals=dict(mandatory=True, + input_map = dict(b0_thres=dict(usedefault=True, ), - bvecs=dict(mandatory=True, + ignore_exception=dict(nohash=True, + usedefault=True, + ), + in_bval=dict(mandatory=True, + ), + in_bvec=dict(mandatory=True, ), in_file=dict(mandatory=True, ), mask_file=dict(), - out_filename=dict(genfile=True, - ), + out_prefix=dict(), ) inputs = DTI.input_spec() diff --git a/nipype/interfaces/dipy/tests/test_auto_Denoise.py b/nipype/interfaces/dipy/tests/test_auto_Denoise.py index dad60ccb5f..6a400231c4 100644 --- a/nipype/interfaces/dipy/tests/test_auto_Denoise.py +++ b/nipype/interfaces/dipy/tests/test_auto_Denoise.py @@ -13,6 +13,8 @@ def test_Denoise_inputs(): usedefault=True, ), patch_radius=dict(), + signal_mask=dict(), + snr=dict(), ) inputs = Denoise.input_spec() diff --git a/nipype/interfaces/dipy/tests/test_auto_DipyBaseInterface.py b/nipype/interfaces/dipy/tests/test_auto_DipyBaseInterface.py new file mode 100644 index 0000000000..ce3bd17584 --- /dev/null +++ b/nipype/interfaces/dipy/tests/test_auto_DipyBaseInterface.py @@ -0,0 +1,16 @@ +# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT +from ....testing import assert_equal +from ..base import DipyBaseInterface + + +def test_DipyBaseInterface_inputs(): + input_map = dict(ignore_exception=dict(nohash=True, + usedefault=True, + ), + ) + inputs = DipyBaseInterface.input_spec() + + for key, metadata in list(input_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(inputs.traits()[key], metakey), value + diff --git a/nipype/interfaces/dipy/tests/test_auto_DipyDiffusionInterface.py b/nipype/interfaces/dipy/tests/test_auto_DipyDiffusionInterface.py new file mode 100644 index 0000000000..e785433355 --- /dev/null +++ b/nipype/interfaces/dipy/tests/test_auto_DipyDiffusionInterface.py @@ -0,0 +1,25 @@ +# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT +from ....testing import assert_equal +from ..base import DipyDiffusionInterface + + +def test_DipyDiffusionInterface_inputs(): + input_map = dict(b0_thres=dict(usedefault=True, + ), + ignore_exception=dict(nohash=True, + usedefault=True, + ), + in_bval=dict(mandatory=True, + ), + in_bvec=dict(mandatory=True, + ), + in_file=dict(mandatory=True, + ), + out_prefix=dict(), + ) + inputs = DipyDiffusionInterface.input_spec() + + for key, metadata in list(input_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(inputs.traits()[key], metakey), value + diff --git a/nipype/interfaces/dipy/tests/test_auto_EstimateResponseSH.py b/nipype/interfaces/dipy/tests/test_auto_EstimateResponseSH.py new file mode 100644 index 0000000000..95e702bd50 --- /dev/null +++ b/nipype/interfaces/dipy/tests/test_auto_EstimateResponseSH.py @@ -0,0 +1,52 @@ +# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT +from ....testing import assert_equal +from ..reconstruction import EstimateResponseSH + + +def test_EstimateResponseSH_inputs(): + input_map = dict(auto=dict(usedefault=True, + xor=['recursive'], + ), + b0_thres=dict(usedefault=True, + ), + fa_thresh=dict(usedefault=True, + ), + ignore_exception=dict(nohash=True, + usedefault=True, + ), + in_bval=dict(mandatory=True, + ), + in_bvec=dict(mandatory=True, + ), + in_evals=dict(mandatory=True, + ), + in_file=dict(mandatory=True, + ), + in_mask=dict(), + out_mask=dict(usedefault=True, + ), + out_prefix=dict(), + recursive=dict(usedefault=True, + xor=['auto'], + ), + response=dict(usedefault=True, + ), + roi_radius=dict(usedefault=True, + ), + ) + inputs = EstimateResponseSH.input_spec() + + for key, metadata in list(input_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(inputs.traits()[key], metakey), value + + +def test_EstimateResponseSH_outputs(): + output_map = dict(out_mask=dict(), + response=dict(), + ) + outputs = EstimateResponseSH.output_spec() + + for key, metadata in list(output_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(outputs.traits()[key], metakey), value diff --git a/nipype/interfaces/dipy/tests/test_auto_RESTORE.py b/nipype/interfaces/dipy/tests/test_auto_RESTORE.py new file mode 100644 index 0000000000..c06bc74573 --- /dev/null +++ b/nipype/interfaces/dipy/tests/test_auto_RESTORE.py @@ -0,0 +1,42 @@ +# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT +from ....testing import assert_equal +from ..reconstruction import RESTORE + + +def test_RESTORE_inputs(): + input_map = dict(b0_thres=dict(usedefault=True, + ), + ignore_exception=dict(nohash=True, + usedefault=True, + ), + in_bval=dict(mandatory=True, + ), + in_bvec=dict(mandatory=True, + ), + in_file=dict(mandatory=True, + ), + in_mask=dict(), + noise_mask=dict(), + out_prefix=dict(), + ) + inputs = RESTORE.input_spec() + + for key, metadata in list(input_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(inputs.traits()[key], metakey), value + + +def test_RESTORE_outputs(): + output_map = dict(evals=dict(), + evecs=dict(), + fa=dict(), + md=dict(), + mode=dict(), + rd=dict(), + trace=dict(), + ) + outputs = RESTORE.output_spec() + + for key, metadata in list(output_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(outputs.traits()[key], metakey), value diff --git a/nipype/interfaces/dipy/tests/test_auto_StreamlineTractography.py b/nipype/interfaces/dipy/tests/test_auto_StreamlineTractography.py new file mode 100644 index 0000000000..b4c4dae679 --- /dev/null +++ b/nipype/interfaces/dipy/tests/test_auto_StreamlineTractography.py @@ -0,0 +1,54 @@ +# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT +from ....testing import assert_equal +from ..tracks import StreamlineTractography + + +def test_StreamlineTractography_inputs(): + input_map = dict(gfa_thresh=dict(mandatory=True, + usedefault=True, + ), + ignore_exception=dict(nohash=True, + usedefault=True, + ), + in_file=dict(mandatory=True, + ), + in_model=dict(), + in_peaks=dict(), + min_angle=dict(mandatory=True, + usedefault=True, + ), + multiprocess=dict(mandatory=True, + usedefault=True, + ), + num_seeds=dict(mandatory=True, + usedefault=True, + ), + out_prefix=dict(), + peak_threshold=dict(mandatory=True, + usedefault=True, + ), + save_seeds=dict(mandatory=True, + usedefault=True, + ), + seed_coord=dict(), + seed_mask=dict(), + tracking_mask=dict(), + ) + inputs = StreamlineTractography.input_spec() + + for key, metadata in list(input_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(inputs.traits()[key], metakey), value + + +def test_StreamlineTractography_outputs(): + output_map = dict(gfa=dict(), + odf_peaks=dict(), + out_seeds=dict(), + tracks=dict(), + ) + outputs = StreamlineTractography.output_spec() + + for key, metadata in list(output_map.items()): + for metakey, value in list(metadata.items()): + yield assert_equal, getattr(outputs.traits()[key], metakey), value diff --git a/nipype/interfaces/dipy/tests/test_auto_TensorMode.py b/nipype/interfaces/dipy/tests/test_auto_TensorMode.py index dc2f0a06ff..53f77d5d33 100644 --- a/nipype/interfaces/dipy/tests/test_auto_TensorMode.py +++ b/nipype/interfaces/dipy/tests/test_auto_TensorMode.py @@ -4,15 +4,19 @@ def test_TensorMode_inputs(): - input_map = dict(bvals=dict(mandatory=True, + input_map = dict(b0_thres=dict(usedefault=True, ), - bvecs=dict(mandatory=True, + ignore_exception=dict(nohash=True, + usedefault=True, + ), + in_bval=dict(mandatory=True, + ), + in_bvec=dict(mandatory=True, ), in_file=dict(mandatory=True, ), mask_file=dict(), - out_filename=dict(genfile=True, - ), + out_prefix=dict(), ) inputs = TensorMode.input_spec() diff --git a/nipype/interfaces/dipy/tests/test_auto_TrackDensityMap.py b/nipype/interfaces/dipy/tests/test_auto_TrackDensityMap.py index 5c2e2d0b8d..187c4c0d49 100644 --- a/nipype/interfaces/dipy/tests/test_auto_TrackDensityMap.py +++ b/nipype/interfaces/dipy/tests/test_auto_TrackDensityMap.py @@ -5,6 +5,9 @@ def test_TrackDensityMap_inputs(): input_map = dict(data_dims=dict(), + ignore_exception=dict(nohash=True, + usedefault=True, + ), in_file=dict(mandatory=True, ), out_filename=dict(usedefault=True, diff --git a/nipype/interfaces/dipy/tracks.py b/nipype/interfaces/dipy/tracks.py index 7daff72589..1b0ad6f381 100644 --- a/nipype/interfaces/dipy/tracks.py +++ b/nipype/interfaces/dipy/tracks.py @@ -7,55 +7,43 @@ >>> os.chdir(datadir) """ import os.path as op -import warnings - +import numpy as np import nibabel as nb import nibabel.trackvis as nbt -from ..base import (TraitedSpec, BaseInterface, BaseInterfaceInputSpec, +from ..base import (TraitedSpec, BaseInterfaceInputSpec, File, isdefined, traits) -from ...utils.filemanip import split_filename -from ...utils.misc import package_check +from .base import DipyBaseInterface from ... import logging -iflogger = logging.getLogger('interface') - -have_dipy = True -try: - package_check('dipy', version='0.6.0') -except Exception as e: - have_dipy = False -else: - from dipy.tracking.utils import density_map +IFLOGGER = logging.getLogger('interface') -class TrackDensityMapInputSpec(TraitedSpec): +class TrackDensityMapInputSpec(BaseInterfaceInputSpec): in_file = File(exists=True, mandatory=True, desc='The input TrackVis track file') reference = File(exists=True, desc='A reference file to define RAS coordinates space') points_space = traits.Enum('rasmm', 'voxel', None, usedefault=True, desc='coordinates of trk file') - voxel_dims = traits.List(traits.Float, minlen=3, maxlen=3, desc='The size of each voxel in mm.') data_dims = traits.List(traits.Int, minlen=3, maxlen=3, desc='The size of the image in voxels.') out_filename = File('tdi.nii', usedefault=True, - desc=('The output filename for the tracks in TrackVis ' - '(.trk) format')) + desc='The output filename for the tracks in TrackVis ' + '(.trk) format') class TrackDensityMapOutputSpec(TraitedSpec): out_file = File(exists=True) -class TrackDensityMap(BaseInterface): +class TrackDensityMap(DipyBaseInterface): """ Creates a tract density image from a TrackVis track file using functions from dipy - Example ------- @@ -70,6 +58,8 @@ class TrackDensityMap(BaseInterface): def _run_interface(self, runtime): from numpy import min_scalar_type + from dipy.tracking.utils import density_map + tracks, header = nbt.read(self.inputs.in_file) streams = ((ii[0]) for ii in tracks) @@ -79,9 +69,9 @@ def _run_interface(self, runtime): data_dims = refnii.shape[:3] kwargs = dict(affine=affine) else: - iflogger.warn(('voxel_dims and data_dims are deprecated' - 'as of dipy 0.7.1. Please use reference ' - 'input instead')) + IFLOGGER.warn( + 'voxel_dims and data_dims are deprecated as of dipy 0.7.1. Please use reference ' + 'input instead') if not isdefined(self.inputs.data_dims): data_dims = header['dim'] @@ -101,13 +91,214 @@ def _run_interface(self, runtime): out_file = op.abspath(self.inputs.out_filename) nb.save(img, out_file) - iflogger.info( - ('Track density map saved as {i}, size={d}, ' - 'dimensions={v}').format(i=out_file, d=img.shape, - v=img.header.get_zooms())) + IFLOGGER.info( + 'Track density map saved as %s, size=%s, dimensions=%s', + out_file, img.shape, img.header.get_zooms()) + return runtime def _list_outputs(self): outputs = self._outputs().get() outputs['out_file'] = op.abspath(self.inputs.out_filename) return outputs + + +class StreamlineTractographyInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc=('input diffusion data')) + in_model = File(exists=True, desc=('input f/d-ODF model extracted from.')) + tracking_mask = File(exists=True, + desc=('input mask within which perform tracking')) + seed_mask = File(exists=True, + desc=('input mask within which perform seeding')) + in_peaks = File(exists=True, desc=('peaks computed from the odf')) + seed_coord = File(exists=True, + desc=('file containing the list of seed voxel ' + 'coordinates (N,3)')) + gfa_thresh = traits.Float(0.2, mandatory=True, usedefault=True, + desc=('GFA threshold to compute tracking mask')) + peak_threshold = traits.Float( + 0.5, mandatory=True, usedefault=True, + desc=('threshold to consider peaks from model')) + min_angle = traits.Float(25.0, mandatory=True, usedefault=True, + desc=('minimum separation angle')) + multiprocess = traits.Bool(True, mandatory=True, usedefault=True, + desc=('use multiprocessing')) + save_seeds = traits.Bool(False, mandatory=True, usedefault=True, + desc=('save seeding voxels coordinates')) + num_seeds = traits.Int(10000, mandatory=True, usedefault=True, + desc=('desired number of tracks in tractography')) + out_prefix = traits.Str(desc=('output prefix for file names')) + + +class StreamlineTractographyOutputSpec(TraitedSpec): + tracks = File(desc='TrackVis file containing extracted streamlines') + gfa = File(desc=('The resulting GFA (generalized FA) computed using the ' + 'peaks of the ODF')) + odf_peaks = File(desc=('peaks computed from the odf')) + out_seeds = File(desc=('file containing the (N,3) *voxel* coordinates used' + ' in seeding.')) + + +class StreamlineTractography(DipyBaseInterface): + + """ + Streamline tractography using EuDX [Garyfallidis12]_. + + .. [Garyfallidis12] Garyfallidis E., “Towards an accurate brain + tractography”, PhD thesis, University of Cambridge, 2012 + + Example + ------- + + >>> from nipype.interfaces import dipy as ndp + >>> track = ndp.StreamlineTractography() + >>> track.inputs.in_file = '4d_dwi.nii' + >>> track.inputs.in_model = 'model.pklz' + >>> track.inputs.tracking_mask = 'dilated_wm_mask.nii' + >>> res = track.run() # doctest: +SKIP + """ + input_spec = StreamlineTractographyInputSpec + output_spec = StreamlineTractographyOutputSpec + + def _run_interface(self, runtime): + from dipy.reconst.peaks import peaks_from_model + from dipy.tracking.eudx import EuDX + from dipy.data import get_sphere + # import marshal as pickle + import cPickle as pickle + import gzip + + if (not (isdefined(self.inputs.in_model) or + isdefined(self.inputs.in_peaks))): + raise RuntimeError(('At least one of in_model or in_peaks should ' + 'be supplied')) + + img = nb.load(self.inputs.in_file) + imref = nb.four_to_three(img)[0] + affine = img.get_affine() + + data = img.get_data().astype(np.float32) + hdr = imref.get_header().copy() + hdr.set_data_dtype(np.float32) + hdr['data_type'] = 16 + + sphere = get_sphere('symmetric724') + + self._save_peaks = False + if isdefined(self.inputs.in_peaks): + IFLOGGER.info('Peaks file found, skipping ODF peaks search...') + f = gzip.open(self.inputs.in_peaks, 'rb') + peaks = pickle.load(f) + f.close() + else: + self._save_peaks = True + IFLOGGER.info('Loading model and computing ODF peaks') + f = gzip.open(self.inputs.in_model, 'rb') + odf_model = pickle.load(f) + f.close() + + peaks = peaks_from_model( + model=odf_model, + data=data, + sphere=sphere, + relative_peak_threshold=self.inputs.peak_threshold, + min_separation_angle=self.inputs.min_angle, + parallel=self.inputs.multiprocess) + + f = gzip.open(self._gen_filename('peaks', ext='.pklz'), 'wb') + pickle.dump(peaks, f, -1) + f.close() + + hdr.set_data_shape(peaks.gfa.shape) + nb.Nifti1Image(peaks.gfa.astype(np.float32), affine, + hdr).to_filename(self._gen_filename('gfa')) + + IFLOGGER.info('Performing tractography') + + if isdefined(self.inputs.tracking_mask): + msk = nb.load(self.inputs.tracking_mask).get_data() + msk[msk > 0] = 1 + msk[msk < 0] = 0 + else: + msk = np.ones(imref.get_shape()) + + gfa = peaks.gfa * msk + seeds = self.inputs.num_seeds + + if isdefined(self.inputs.seed_coord): + seeds = np.loadtxt(self.inputs.seed_coord) + + elif isdefined(self.inputs.seed_mask): + seedmsk = nb.load(self.inputs.seed_mask).get_data() + assert(seedmsk.shape == data.shape[:3]) + seedmsk[seedmsk > 0] = 1 + seedmsk[seedmsk < 1] = 0 + seedps = np.array(np.where(seedmsk == 1), dtype=np.float32).T + vseeds = seedps.shape[0] + nsperv = (seeds // vseeds) + 1 + IFLOGGER.info(('Seed mask is provided (%d voxels inside ' + 'mask), computing seeds (%d seeds/voxel).') % + (vseeds, nsperv)) + if nsperv > 1: + IFLOGGER.info(('Needed %d seeds per selected voxel ' + '(total %d).') % (nsperv, vseeds)) + seedps = np.vstack(np.array([seedps] * nsperv)) + voxcoord = seedps + np.random.uniform(-1, 1, size=seedps.shape) + nseeds = voxcoord.shape[0] + seeds = affine.dot(np.vstack((voxcoord.T, + np.ones((1, nseeds)))))[:3, :].T + + if self.inputs.save_seeds: + np.savetxt(self._gen_filename('seeds', ext='.txt'), seeds) + + if isdefined(self.inputs.tracking_mask): + tmask = msk + a_low = 0.1 + else: + tmask = gfa + a_low = self.inputs.gfa_thresh + + eu = EuDX(tmask, + peaks.peak_indices[..., 0], + seeds=seeds, + affine=affine, + odf_vertices=sphere.vertices, + a_low=a_low) + + ss_mm = [np.array(s) for s in eu] + + trkfilev = nb.trackvis.TrackvisFile( + [(s, None, None) for s in ss_mm], points_space='rasmm', affine=np.eye(4)) + trkfilev.to_file(self._gen_filename('tracked', ext='.trk')) + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + outputs['tracks'] = self._gen_filename('tracked', ext='.trk') + outputs['gfa'] = self._gen_filename('gfa') + if self._save_peaks: + outputs['odf_peaks'] = self._gen_filename('peaks', ext='.pklz') + if self.inputs.save_seeds: + if isdefined(self.inputs.seed_coord): + outputs['out_seeds'] = self.inputs.seed_coord + else: + outputs['out_seeds'] = self._gen_filename('seeds', + ext='.txt') + + return outputs + + def _gen_filename(self, name, ext=None): + fname, fext = op.splitext(op.basename(self.inputs.in_file)) + if fext == '.gz': + fname, fext2 = op.splitext(fname) + fext = fext2 + fext + + if not isdefined(self.inputs.out_prefix): + out_prefix = op.abspath(fname) + else: + out_prefix = self.inputs.out_prefix + + if ext is None: + ext = fext + + return out_prefix + '_' + name + ext diff --git a/nipype/testing/data/dilated_wm_mask.nii b/nipype/testing/data/dilated_wm_mask.nii new file mode 100644 index 0000000000..e69de29bb2 diff --git a/nipype/testing/data/dwi_evals.nii b/nipype/testing/data/dwi_evals.nii new file mode 100644 index 0000000000..e69de29bb2 diff --git a/nipype/testing/data/model.pklz b/nipype/testing/data/model.pklz new file mode 100644 index 0000000000..e69de29bb2