From 674c96ea4b8e472362e268dfe1e106c1484f91da Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 28 Feb 2013 10:15:32 +0100 Subject: [PATCH 01/11] First implementation for MBIS --- nipype/interfaces/mbis/__init__.py | 0 nipype/interfaces/mbis/postproc.py | 236 +++++++++++++++++++ nipype/interfaces/mbis/segmentation.py | 303 +++++++++++++++++++++++++ 3 files changed, 539 insertions(+) create mode 100644 nipype/interfaces/mbis/__init__.py create mode 100644 nipype/interfaces/mbis/postproc.py create mode 100644 nipype/interfaces/mbis/segmentation.py diff --git a/nipype/interfaces/mbis/__init__.py b/nipype/interfaces/mbis/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/nipype/interfaces/mbis/postproc.py b/nipype/interfaces/mbis/postproc.py new file mode 100644 index 0000000000..ef47c30098 --- /dev/null +++ b/nipype/interfaces/mbis/postproc.py @@ -0,0 +1,236 @@ + +def loadParameters( fname ): + import csv + import numpy as np + + data = [] + with open( fname ) as csvfile: + for row in csvfile: + row = row.replace("[","").replace("]","").replace( " ", "") + data.append( [ float(val) for val in row.split(',') ] ) + + nparam = len( data[0] ) + nclasses = int( 0.5 * (nparam*4+1)**0.5 - 0.5 ) + # Load parameters + mu = [] + sigma = [] + + for row in data: + mu.append( np.matrix( [ float(val) for val in row[0:nclasses] ] ) ) + sigma.append( np.matrix( np.reshape( [ float(val) for val in row[nclasses:] ], (nclasses,nclasses) ) ) ) + + # Return arrays + return ( mu, sigma ) + +def distancePV ( sample, mask, params_tissue1, params_tissue2, distance ): + from scipy.spatial.distance import mahalanobis,euclidean + import numpy as np + + # Direction vector between pure tissues + d_vect = np.ravel(params_tissue2[0] - params_tissue1[0]).T + mu1 = np.ravel(params_tissue1[0]) + mu2 = np.ravel(params_tissue2[0]) + SI1 = params_tissue1[1].getI() + SI2 = params_tissue2[1].getI() + + if distance=='mahalanobis': + norm = np.array( [ 1/(1+ mahalanobis(pix,mu2,SI2)/ mahalanobis(pix,mu1,SI1)) for pix in sample[mask==1] ] ) + elif distance=='dummy': + norm = mask*0.5 + else: + norm = np.array( [ 1/(1+ euclidean(pix,mu2)/ euclidean(pix,mu1)) for pix in sample[mask==1] ] ) + result = np.zeros( np.shape( mask ) ) + result[mask==1] = norm + return result + +def computeMAP( tpms, normalize=True ): + import numpy as np + normalizer = np.sum( tpms, axis=0 ) + if normalize: + for tpm in tpms: + tpm[normalizer>0]/= normalizer[normalizer>0] + + out_seg = np.zeros( shape=tpms[0].shape, dtype=np.uint8 ) + out_seg[normalizer>0] = np.argmax( tpms, axis=0 )[normalizer>0] + out_seg[normalizer>0] += 1 + + return out_seg, tpms + +def fusePV2( in_files, in_maps, parameters, pt_list=( (0,1,2), (2,3), (4,) ), distance='mahalanobis', reorder=True, prefix='./' ): + import nibabel as nib + import numpy as np + import os + import sys + import collections + from scipy.spatial.distance import mahalanobis,euclidean + + nmaps = len( in_maps ) + nclasses = np.shape( parameters )[1] + assert nmaps == nclasses + + idxs = np.array(range( 0, nclasses )) + corder = range( 0, nclasses ) + npts = len( pt_list ) + + pv_task = np.array([ len(t)>1 for t in pt_list ]) # Check what tissues have pv actions + + if not pv_task.any(): + assert( nmaps == npts ) + return in_maps + + # Load probability maps + initmaps = [ nib.load(f).get_data() for f in in_maps ] + + # If reorder is True, find unordered tissue signatures + if reorder: + means = parameters[0] + firstmeans = np.array( [ np.ravel(val)[0] for val in means ] ) + idxs = np.argsort(firstmeans) + corder = np.take(corder, idxs) + + # Load images + channels = [ nib.load(c) for c in in_files ] + + # Prepare sample + data_shape = np.shape( channels[0] ) + data_sample = [ channel.get_data().reshape(-1) for channel in channels ] + data_sample = np.swapaxes( data_sample, 0, 1) + + t_flatten = y=collections.Counter( [element for tupl in pt_list for element in tupl] ) + pt_mapping = [ [ i, [] ] for i in range(0,npts) ] + #pv_mapping = [ [ np.where(idxs==i)[0] , [], i ] for i in t_flatten if t_flatten[i]>1 ] + pv_mapping = [ [ idxs[i] , [] ] for i in t_flatten if t_flatten[i]>1 ] + pt_maps = [ np.zeros( shape=data_shape, dtype=float) for i in range(0,npts)] + for pt_tuple,out_i in zip( pt_list, range(0,npts) ): + for t in pt_tuple: + in_i = idxs[t] + if t_flatten[t]==1: + pt_mapping[out_i][1].append( in_i ) + pt_maps[out_i] += initmaps[ in_i ] + else: + pvm_id = pv_mapping[:][0].index(in_i) + pv_mapping[pvm_id][1].append( out_i ) + + for m in pv_mapping: + in_id = m[0] + neighs = [ pt_mapping[val][1] for val in m[1] ] + dist = [] + mask = np.zeros(data_shape) + mask[initmaps[in_id]>0.001] = 1 + mu = [] + SI = [] + out_ids = [] + for n,l in zip(neighs, m[1][:]): + out_ids.append(l) + mu.append( parameters[0][n[-1]] ) + SI.append( parameters[1][n[-1]].getI() ) + + diff = np.zeros( shape=np.shape(mask.reshape(-1)), dtype=float) + norm = np.array([1/(1+ mahalanobis(pix,mu[1],SI[1])/ mahalanobis(pix,mu[0],SI[0])) for pix in data_sample[mask.reshape(-1)==1]]) + diff[mask.reshape(-1)==1] = norm + pt_maps[out_ids[0]]+= diff.reshape(data_shape) * initmaps[in_id] + pt_maps[out_ids[1]]+= (1-diff.reshape(data_shape)) * initmaps[in_id] + + # Normalize tpms and compute MAP + normalizer = np.sum( pt_maps, axis=0 ) + for pt_map in pt_maps: + pt_map[normalizer>0]/= normalizer[normalizer>0] + + # Generate output names + ppmnames = [ '%s_pvseg%02d.nii.gz' % ( prefix, i ) for i in range(0,len(pt_maps)) ] + + # Save + ref = nib.load( in_maps[0] ) + ref = channels[0] + for i in range(0,len(pt_maps)): + nii = nib.Nifti1Image( np.reshape( pt_maps[i], data_shape), ref.get_affine(), ref.get_header() ) + nib.save( nii, ppmnames[i] ) + + return ppmnames + + +def fusePV( in_files, in_maps, parameters, pt_list=[ 0, 2, 4 ], distance='mahalanobis', reorder=True, prefix='./' ): + import nibabel as nib + import numpy as np + import os + + nmaps = len( in_maps ) + nclasses = np.shape( parameters )[1] + assert nmaps == nclasses + + corder = range( 0, nclasses ) + npts = len( pt_list ) + + if npts == nmaps: + return in_maps + + pt_list = np.sort( pt_list ) + + # Load probability maps + initmaps = [ nib.load(f) for f in in_maps ] + + # If reorder is True, find unordered tissue signatures + if reorder: + means = parameters[0] + firstmeans = np.array( [ np.ravel(val)[0] for val in means ] ) + m = np.argsort(firstmeans) + corder = np.take(corder, m) + new_idx = np.take( corder, pt_list ) + + # Load images + channels = [ nib.load(c) for c in in_files ] + + # Prepare sample + data_shape = np.shape( channels[0] ) + data_sample = [ channel.get_data().reshape(-1) for channel in channels ] + data_sample = np.swapaxes( data_sample, 0, 1) + + + # Split between pv (partial volume) and pt (pure tissue) maps + pt_niis = [] + pv_niis = [] + pt_param = [] + + for tmap,i in zip(initmaps,range(0,nclasses)): + idx = np.where( new_idx==i )[0] + if len(idx) == 1: + pt_niis.append( tmap ) + pt_param.append( [ parameters[0][i], parameters[1][i] ] ) + else: + pv_niis.append( tmap ) + + # Compute the steps required + steps = [ val-pt_list[i-1]-1 for val,i in zip( pt_list[1:], range(1,len(pt_list[1:])+1 ) ) ] + + # Extract data and initialize normalizer + pt_maps = [ m.get_data().reshape(-1) for m in pt_niis ] + pv_maps = [ m.get_data().reshape(-1) for m in pv_niis ] + + # Process + for pt_map,i in zip( pt_maps[:-1],range(0,len(pt_maps[:-1])) ): + curr_steps = steps[i] + for step in range(1,curr_steps+1): + pv_idx = (step-1)+i + pv_map = pv_maps[pv_idx] + mask = np.zeros( np.shape( pt_map ) ) + mask[pv_map>0.001] = 1 + + if not step == curr_steps: # Direct addition of the pv map to the last pure tissue map + pt_map+= pv_map + else: # Split pv fraction proportionally to the distance to a contiguous pure tissue + dist = distancePV( data_sample, mask, pt_param[i], pt_param[i+1], distance ) + pt_map+= dist * pv_map + pt_maps[i+1]+= (1-dist) * pv_map + + # Normalize tpms and compute MAP + normalizer = np.sum( pt_maps, axis=0 ) + + # Generate output names + ppmnames = [ '%s_pvseg%02d.nii.gz' % ( prefix[:-6], i ) for i in range(0,len(pt_niis)) ] + + # Save + for i in range(0,len(pt_niis)): + nii = nib.Nifti1Image( np.reshape( pt_maps[i], data_shape) , pt_niis[0].get_affine(), pt_niis[0].get_header() ) + nib.save( nii, ppmnames[i] ) + + return ppmnames diff --git a/nipype/interfaces/mbis/segmentation.py b/nipype/interfaces/mbis/segmentation.py new file mode 100644 index 0000000000..d60d2d9c2a --- /dev/null +++ b/nipype/interfaces/mbis/segmentation.py @@ -0,0 +1,303 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" MBIS nipype interface definition + Based on FAST interface definition + + 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, os.path as op +import warnings +import numpy as np +import nibabel as nib +from nipype.interfaces.base import (TraitedSpec, File, InputMultiPath, + OutputMultiPath, Undefined, traits, + isdefined, OutputMultiPath, + CommandLineInputSpec, CommandLine, + BaseInterface, BaseInterfaceInputSpec, + traits ) +from nipype.utils.filemanip import split_filename,fname_presuffix +import postproc as pp +#import csv + +warn = warnings.warn +warnings.filterwarnings('always', category=UserWarning) + + + +class MBISInputSpec( CommandLineInputSpec ): + in_files = InputMultiPath( File(exists=True), copyfile=False, + desc='image, or multi-channel set of images, ' \ + 'to be segmented', + argstr='-C %s', position=-1, mandatory=True ) + + mask = File(exists=True, desc='binary mask file', argstr='-x %s' ) + + mask_auto = traits.Bool( desc='channels are implicitly masked', argstr='-M' ) + + out_prefix = File('outpath', desc='base name of output files', + argstr='-o %s', genfile=True) # uses in_file name as basename if none given + + initialization = traits.Enum( 'none','em','km','km+em', argstr='-I %s',desc='Initialization mode' ) + + number_classes = traits.Range(low=2, high=10, argstr='-n %d', + desc='number of tissue-type classes', value=3) + + output_steps = traits.Bool(desc='output intermediate steps', + argstr='--output-steps') + + + output_biasfield = traits.Bool(desc='output estimated bias field', + argstr='--bias-output') + + output_biascorrected = traits.Bool(desc='output restored image ' \ + '(bias-corrected image)', + argstr='--bias-corrected-output') + + output_stats = File( 'outcsvfile', desc='output file containing mixture parameters', argstr='--output-stats %s' ) + + probability_maps = traits.Bool(desc='outputs a separate binary image for each ' \ + 'tissue type', + argstr='-g' ) + + priors = InputMultiPath(File(exist=True), desc='initialize with prior images', + argstr='-P %s', minlen=3, maxlen=10) + + no_bias = traits.Bool(desc='do not remove bias field', + argstr='--bias-skip' ) + + no_normalize = traits.Bool( desc='do not normalize input channels', argstr='-N', + value=False ) + + em_iters = traits.Range(low=1, high=50, value=3, + desc='number of EM iterations', + argstr='--em-iterations %d') + mrf_iters = traits.Range(low=1, high=10, + desc='number of MRF iterations', + argstr='--mrf-iterations %d') + + mrf_lambda = traits.Range(low=0.01, high=3.0, + desc='MRF lambda parameter (segmentation spatial smoothness)', + argstr='-l %.3f') + + manual_init = File(exists=True, desc='Filename containing intensities', + argstr='-f %s') + manual_init_means = File(exists=True, desc='Filename containing intensities', + argstr='--init-means %s') + +# no_pve = traits.Bool(desc='turn off PVE (partial volume estimation)', +# argstr='--nopve') + +# use_priors = traits.Bool(desc='use priors throughout', +# argstr='-P') # must also set -a!, +# # mutually inclusive?? +# # No, conditional +# # mandatory... need to +# # figure out how to +# # handle with traits. +# verbose = traits.Bool(desc='switch on diagnostic messages', +# argstr='-v') + + + +class MBISOutputSpec( TraitedSpec ): + out_classified = File(desc='path/name of binary segmented volume file' \ + ' one val for each class _mrf') + out_parameters = traits.File(desc='csv file with tissue purameters') + bias_field = OutputMultiPath(File(desc='Estimated bias field _bias')) + probability_maps = OutputMultiPath(File(desc='filenames, one for each class, for each ' \ + 'input, mrf_x')) + restored_image = OutputMultiPath(File(desc='restored images (one for each input image) ' \ + 'named according to the input images _corrected_chXX')) + + normalized_inputs = OutputMultiPath( File( desc='normalized input channels' ) ) + + +class MBIS(CommandLine): + """ Use MBIS for segmentation and bias correction. + + Examples + -------- + >>> from nipype.interfaces import mbis + >>> from nipype.testing import example_data + + Assign options through the ``inputs`` attribute: + + >>> mbisr = mbis.MBIS() + >>> mbisr.inputs.in_files = example_data('structural.nii') + >>> out = mbisr.run() #doctest: +SKIP + + """ + _cmd = 'brain_seg' + input_spec = MBISInputSpec + output_spec = MBISOutputSpec + + def _list_outputs(self): + outputs = self.output_spec().get() + if not isdefined(self.inputs.number_classes): + nclasses = 3 + else: + nclasses = self.inputs.number_classes + + outprefix = self._getdefaultprefix() + + outputs['out_classified'] = '%s_mrf.nii.gz' % outprefix + + if self.inputs.probability_maps: + outputs['probability_maps'] = [] + for i in range(1,nclasses+1): + outputs['probability_maps'].append('%s_mrf_seg%02d.nii.gz' % (outprefix,i) ) + + if (not (isdefined(self.inputs.no_bias))) and self.inputs.output_biascorrected: + outputs['restored_image'] = [] + for val, f in enumerate(self.inputs.in_files): + # image numbering is 1-based + outputs['restored_image'].append('%s_corrected_ch%02d.nii.gz' % (outprefix,val+1) ) + + if (not (isdefined(self.inputs.no_bias))) and self.inputs.output_biasfield: + outputs['bias_field'] = [] + for val, f in enumerate(self.inputs.in_files): + # image numbering is 1-based + outputs['bias_field'].append('%s_bias_ch%02d.nii.gz' % (outprefix,val) ) + + if ((not self.inputs.no_normalize) and self.inputs.output_steps ): + outputs['normalized_inputs' ] = [] + for val, f in enumerate(self.inputs.in_files): + outputs['normalized_inputs'].append('%s_normin_%0d.nii.gz' % (outprefix,val) ) + + if isdefined(self.inputs.output_stats): + fname= outprefix + '_stats_final' + self.inputs.output_stats + outputs['out_parameters'] = fname +# outputs['out_parameters'] = np.loadtxt( fname, delimiter='[],' ) +# with open( self.inputs.output_stats ) as csvfile: +# dataReader = csv.reader( csvfile ) +# outputs['out_parameters'] = np.array( [ [ row ] for row in dataReader ] ) +# csvfile.close() + + return outputs + + def _getdefaultprefix( self, name='mbis' ): + if not isdefined(self.inputs.out_prefix): + return os.path.abspath( os.path.join( os.getcwd(), name ) ) + else: + return self.inputs.out_prefix + + def _gen_filename(self,name): + if name == 'out_prefix': + return self._getdefaultprefix() + return None + + + def _gen_fname(self, prefix, suffix=None, ext='.nii.gz', cwd=None): + """Generate a filename based on the given parameters. + + The filename will take the form: preffix. + + Parameters + ---------- + prefix : str + Filename to base the new filename on. + + suffix : str + Suffix to add to the `basename`. (defaults is '' ) + + ext : str + Desired extension (default is nii.gz) + + Returns + ------- + fname : str + New filename based on given parameters. + + """ + if ((prefix == '') or (prefix is None) ): + prefix = './' + if suffix is None: + suffix = '' + if cwd is None: + cwd = os.getcwd() + + suffix = ''.join((suffix,ext)) + fname = fname_presuffix(prefix, suffix=suffix, use_ext=False, newpath=cwd ) + return fname + + +class PVMergeInputSpec(BaseInterfaceInputSpec): + in_files = InputMultiPath( File(exists=True), copyfile=False, desc='image, or multi-channel set of images,' \ + 'used for segmentation', mandatory=True ) + in_maps = InputMultiPath( File(exists=True), copyfile=False, desc='Resulting tissue probability maps', mandatory=True ) + in_labelling = File( exists=True, desc='Hard labelling resulting from segmentation', mandatory=True) + parameters = File( exists=True, desc='CSV file containing parameters for all the classes', mandatory=True) + pure_tissues = traits.Tuple( value=( (0,), (1,), (2,) ), minlen=1, desc='Identifiers of the pure tissue classes' ) + dist = traits.String( value="euclidean", desc="Distance definition to be used" ) + reorder = traits.Bool( value=True, desc='Reorder classes if the classification is not ordered by first contrast means' ) + +class PVMergeOutputSpec(TraitedSpec): + out_files = OutputMultiPath( File( desc='filenames, one for each pure tissue according to prefix_msegXX.nii.gz' ) ) + +class PVMerge(BaseInterface): + input_spec = PVMergeInputSpec + output_spec = PVMergeOutputSpec + _ppmnames = [] + + def _run_interface(self, runtime): + params = pp.loadParameters(self.inputs.parameters) + path, base, _ = split_filename( self.inputs.in_maps[0] ) + + self._ppmnames = pp.fusePV2( self.inputs.in_files, + self.inputs.in_maps, + params, + self.inputs.pure_tissues, + self.inputs.dist, + self.inputs.reorder, + os.path.abspath( os.path.join( os.getcwd(), base ) ) ) + + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + outputs["out_files"] = self._ppmnames + return outputs + +class MAPInputSpec(BaseInterfaceInputSpec): + in_files = InputMultiPath( File(exists=True), copyfile=False, desc='Input tissue probability maps', mandatory=True ) + normalize = traits.Bool( desc='apply normalization to input tpms', value=True ) + prefix = traits.String( desc='prefix path and basename for output', value='' ) + +class MAPOutputSpec(TraitedSpec): + out_tpms = OutputMultiPath( File( desc='Normalized tpms' ) ) + out_file = File( desc='output labelling after MAP' ) + +class MAP(BaseInterface): + input_spec = MAPInputSpec + output_spec = MAPOutputSpec + + def _run_interface(self, runtime): + import nibabel as nib + import numpy as np + path, base, _ = split_filename( self.inputs.in_files[0] ) + prefix=self.inputs.prefix + if not isdefined(self.inputs.prefix): + prefix = os.getcwd() + "/"+ base + self._out_fname = '%s_map.nii.gz' % prefix + assert( len(self.inputs.in_files) > 1 ) + tpm_data = np.array( [ nib.load(tpm).get_data() for tpm in self.inputs.in_files ] ) + (seg_file,tpms_norm) = pp.computeMAP( tpm_data, self.inputs.normalize ) + ref = nib.load(self.inputs.in_files[0] ) + nib.save( nib.Nifti1Image( seg_file, ref.get_affine(), ref.get_header() ), self._out_fname ) + if self.inputs.normalize: + self._tpmnames = [ '%s_tpm%02d.nii.gz' % (prefix,i) for i in range(1,len(self.inputs.in_files)+1) ] + for tpm,name in zip(tpms_norm,self._tpmnames): + nib.save( nib.Nifti1Image( tpm, ref.get_affine(), ref.get_header() ), name ) + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + if self.inputs.normalize: + outputs["out_tpms"] = self._tpmnames + outputs["out_file"] = self._out_fname + return outputs From ad8f4f6754c9373c3601aa7b5a4029206834e488 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 28 Feb 2013 10:21:46 +0100 Subject: [PATCH 02/11] Standardized import to add all utilities at once --- nipype/interfaces/mbis/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/nipype/interfaces/mbis/__init__.py b/nipype/interfaces/mbis/__init__.py index e69de29bb2..931c115a0b 100644 --- a/nipype/interfaces/mbis/__init__.py +++ b/nipype/interfaces/mbis/__init__.py @@ -0,0 +1,2 @@ + +from .segmentation import (MBIS,PVMerge,MAP) From 449d6dc7e5abb79ababcf96e1d3d89a74aef1f2b Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 28 Feb 2013 10:32:27 +0100 Subject: [PATCH 03/11] Added MBIS to sub-package configuration list --- nipype/interfaces/setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nipype/interfaces/setup.py b/nipype/interfaces/setup.py index b4fcb6dc36..430947e3de 100644 --- a/nipype/interfaces/setup.py +++ b/nipype/interfaces/setup.py @@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None): config.add_subpackage('dipy') config.add_subpackage('freesurfer') config.add_subpackage('fsl') + config.add_subpackage('mbis') config.add_subpackage('mne') config.add_subpackage('mrtrix') config.add_subpackage('nipy') From f5c1be61c5af10730db9a4ffc061933e2ada835c Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 28 Feb 2013 11:52:40 +0100 Subject: [PATCH 04/11] Minor fixes * Added pipeline to __init__.py so its loadable the same way the others are * Polished code in pipeline creation --- nipype/workflows/dmri/fsl/__init__.py | 4 ++-- nipype/workflows/dmri/fsl/dti.py | 20 +------------------- 2 files changed, 3 insertions(+), 21 deletions(-) diff --git a/nipype/workflows/dmri/fsl/__init__.py b/nipype/workflows/dmri/fsl/__init__.py index 80c1ecc4f2..e25c80bdc0 100644 --- a/nipype/workflows/dmri/fsl/__init__.py +++ b/nipype/workflows/dmri/fsl/__init__.py @@ -1,4 +1,4 @@ -from dti import create_bedpostx_pipeline, create_eddy_correct_pipeline +from dti import create_bedpostx_pipeline, create_eddy_correct_pipeline, create_epi_correct_pipeline from tbss import (create_tbss_1_preproc, create_tbss_2_reg, create_tbss_3_postreg, create_tbss_4_prestats, - create_tbss_all, create_tbss_non_FA) \ No newline at end of file + create_tbss_all, create_tbss_non_FA) diff --git a/nipype/workflows/dmri/fsl/dti.py b/nipype/workflows/dmri/fsl/dti.py index 1cce5a5a52..3945eed6ea 100644 --- a/nipype/workflows/dmri/fsl/dti.py +++ b/nipype/workflows/dmri/fsl/dti.py @@ -319,7 +319,7 @@ def create_epi_correct_pipeline(name="epi_correct", register_to_ref=False ): ,(mask_mag_dil, prelude, [('out_file', 'mask_file' ) ]) ,(prelude, fill_phase, [('unwrapped_phase_file','in_file')] ) ,(inputnode, vsm, [('fieldmap_mag', 'in_file')]) - ,(fill_phase, vsm, [('out_file','phasemap_file')]) #, (('out_file',_genvsmpath),'shift_out_file') ] ) + ,(fill_phase, vsm, [('out_file','phasemap_file')]) ,(inputnode, vsm, [(('te_diff',_ms2sec),'asym_se_time'),('vsm_sigma','smooth2d')]) ,(dwell_time, vsm, [(('dwell_time',_ms2sec),'dwell_time') ]) ,(mask_mag_dil, vsm, [('out_file','mask_file')] ) @@ -444,27 +444,9 @@ def _vsm_remove_mean( in_file, mask_file, in_unwarped ): nib.save( img, out_file ) return out_file -#def _generate_fwdmap( in_file, mask_file, shift_in_file, in_unwarped): -# import os -# out_file = os.path.abspath( './fwdmap.nii.gz' ) -# cmd = 'fugue -i %s -w %s --loadshift=%s --mask=%s' % ( in_file, out_file, shift_in_file, mask_file ) -# print " Running:",cmd -# os.system(cmd) -# return out_file - -def _genvsmpath( in_file ): - import os - name,fext = os.path.splitext( os.path.basename(in_file) ) - if fext == '.gz': name,_ = os.path.splitext(name) - out_file = os.path.abspath( os.path.join( os.path.dirname(in_file), '%s_vsm.nii.gz' % name )) - return out_file - def _ms2sec( val ): return val*1e-3; -def _yield( in_file, in_block ): - return in_file - def _split_dwi( in_file ): import nibabel as nib import os From e2718e9b4629833e942f01815fef3663e962a38d Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Mon, 4 Mar 2013 12:47:26 +0100 Subject: [PATCH 05/11] Polished code Removed some useless functions and improved susceptibility pipeline --- nipype/workflows/dmri/fsl/dti.py | 1 - 1 file changed, 1 deletion(-) diff --git a/nipype/workflows/dmri/fsl/dti.py b/nipype/workflows/dmri/fsl/dti.py index 3945eed6ea..3d321cfbca 100644 --- a/nipype/workflows/dmri/fsl/dti.py +++ b/nipype/workflows/dmri/fsl/dti.py @@ -339,7 +339,6 @@ def create_epi_correct_pipeline(name="epi_correct", register_to_ref=False ): # Select reference volume from EPI (B0 in dMRI and a middle frame in fMRI) select_epi = pe.Node( interface=fsl.utils.ExtractROI(t_size=1), name="select_epi" ) - # fugue -i %s -w %s --loadshift=%s --mask=%s % ( mag_name, magfw_name, vsmmag_name, mask_name ), log ) # Forward Map vsm_fwd = pe.Node( interface=fsl.FUGUE(save_warped=True), name="vsm_fwd") vsm_reg = pe.Node( interface=fsl.FLIRT( bins=256, cost='corratio', dof=6, interp='trilinear', searchr_x=[-10,10], searchr_y=[-10,10], searchr_z=[-10,10] ), name="vsm_registration" ) From 660852479fb99fc5ada3be097c874c20a4573862 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Mon, 4 Mar 2013 16:33:57 +0100 Subject: [PATCH 06/11] Modified FSL FLIRT interface to provide log support In order that supporting logging for fdt_rotate_bvecs script, the verbose output of FLIRT is saved to file if the save_log input is provided. --- nipype/interfaces/fsl/preprocess.py | 35 +++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/nipype/interfaces/fsl/preprocess.py b/nipype/interfaces/fsl/preprocess.py index c16649ac61..362818010e 100644 --- a/nipype/interfaces/fsl/preprocess.py +++ b/nipype/interfaces/fsl/preprocess.py @@ -391,6 +391,9 @@ class FLIRTInputSpec(FSLCommandInputSpec): out_matrix_file = File(argstr='-omat %s', desc='output affine matrix in 4x4 asciii format', genfile=True, position=3, hash_files=False) + + out_log = File( desc='output log' ) + in_matrix_file = File(argstr='-init %s', desc='input 4x4 affine matrix') apply_xfm = traits.Bool(argstr='-applyxfm', requires=['in_matrix_file'], desc='apply transformation supplied by in_matrix_file') @@ -461,6 +464,8 @@ class FLIRTInputSpec(FSLCommandInputSpec): desc='do not use blurring on downsampling') rigid2D = traits.Bool(argstr='-2D', desc='use 2D rigid body mode - ignores dof') + + save_log = traits.Bool(desc='save to log file') verbose = traits.Int(argstr='-verbose %d', desc='verbose mode, 0 is least') @@ -471,6 +476,7 @@ class FLIRTOutputSpec(TraitedSpec): out_matrix_file = File(exists=True, desc='path/name of calculated affine transform ' \ '(if generated)') + out_log = File(desc='path/name of output log (if generated)' ) class FLIRT(FSLCommand): @@ -508,14 +514,43 @@ def _list_outputs(self): outputs['out_file'] = os.path.abspath(outputs['out_file']) outputs['out_matrix_file'] = self.inputs.out_matrix_file + # Generate an out_matrix file if one is not provided if not isdefined(outputs['out_matrix_file']): outputs['out_matrix_file'] = self._gen_fname(self.inputs.in_file, suffix='_flirt.mat', change_ext=False) outputs['out_matrix_file'] = os.path.abspath(outputs['out_matrix_file']) + + if isdefined( self.inputs.save_log ) and self.inputs.save_log: + outputs['out_log'] = self.inputs.out_log + # Generate an out_log file if one is not provided + if not isdefined( outputs['out_log']): + outputs['out_log'] = self._gen_fname( self.inputs.in_file, + suffix='_flirt.log', + change_ext=False) + outputs['out_log'] = os.path.abspath( outputs['out_log'] ) + return outputs + + def aggregate_outputs(self, runtime=None, needed_outputs=None): + outputs = super(FLIRT,self).aggregate_outputs(runtime=runtime, needed_outputs=needed_outputs) + if isdefined( self.inputs.save_log ) and self.inputs.save_log: + with open(outputs.out_log, "a") as text_file: + text_file.write(runtime.stdout + '\n' ) return outputs + def _parse_inputs(self, skip=None): + skip = [] + + if isdefined( self.inputs.save_log ) and self.inputs.save_log: + if not isdefined( self.inputs.verbose ) or self.inputs.verbose==0: + self.inputs.verbose=1 + + skip.append( 'save_log' ) + + return super(FLIRT,self)._parse_inputs(skip=skip) + + def _gen_filename(self, name): if name in ('out_file', 'out_matrix_file'): return self._list_outputs()[name] From 04c568ced928de0ee4855b7ba6e7502917a21559 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 5 Mar 2013 19:14:26 +0100 Subject: [PATCH 07/11] Added motion correction pipeline * Added motion correction pipeline, rotating the b-matrix * Some minor fixes on the rest of pipelines --- nipype/workflows/dmri/fsl/__init__.py | 2 +- nipype/workflows/dmri/fsl/dti.py | 154 ++++++++++++++++++++++---- 2 files changed, 136 insertions(+), 20 deletions(-) diff --git a/nipype/workflows/dmri/fsl/__init__.py b/nipype/workflows/dmri/fsl/__init__.py index e25c80bdc0..db5679e7da 100644 --- a/nipype/workflows/dmri/fsl/__init__.py +++ b/nipype/workflows/dmri/fsl/__init__.py @@ -1,4 +1,4 @@ -from dti import create_bedpostx_pipeline, create_eddy_correct_pipeline, create_epi_correct_pipeline +from dti import create_bedpostx_pipeline, create_eddy_correct_pipeline, create_susceptibility_correct_pipeline from tbss import (create_tbss_1_preproc, create_tbss_2_reg, create_tbss_3_postreg, create_tbss_4_prestats, create_tbss_all, create_tbss_non_FA) diff --git a/nipype/workflows/dmri/fsl/dti.py b/nipype/workflows/dmri/fsl/dti.py index 3d321cfbca..2ac6a7403f 100644 --- a/nipype/workflows/dmri/fsl/dti.py +++ b/nipype/workflows/dmri/fsl/dti.py @@ -12,6 +12,41 @@ def transpose(samples_over_fibres): a = a.reshape(-1,1) return a.T.tolist() +def create_dmri_preprocessing(name="dMRI_preprocessing"): + """Creates a workflow that chains the necessary pipelines to + correct for motion, eddy currents, and susceptibility + artifacts in EPI dMRI sequences. + + Example + ------- + + >>> nipype_dmri_preprocess = create_dmri_preprocessing("nipype_dmri_prep") + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.run() # doctest: +SKIP + + Inputs:: + + inputnode. + + Outputs:: + + outputnode.dmri_corrected + + """ + + pipeline = pe.Workflow(name=name) + + return pipeline + + + def create_bedpostx_pipeline(name="bedpostx"): """Creates a pipeline that does the same as bedpostx script from FSL - calculates diffusion model parameters (distributions not MLE) voxelwise for @@ -168,6 +203,61 @@ def create_bedpostx_pipeline(name="bedpostx"): ]) return bedpostx +def create_motion_correct_pipeline(name="motion_correct"): + """Creates a pipeline that corrects for motion artifact in dMRI sequences. + It takes a series of diffusion weighted images and rigidly corregisters + them to one reference image. Finally, the b-matrix is rotated accordingly, + making use of the rotation matrix obtained by FLIRT. + + Example + ------- + + >>> nipype_motioncorrect = create_motion_correct_pipeline("nipype_motioncorrect") + >>> nipype_motioncorrect.inputs.inputnode.in_file = 'diffusion.nii' + >>> nipype_motioncorrect.inputs.inputnode.in_bvec = 'diffusion.bvec' + >>> nipype_motioncorrect.inputs.inputnode.ref_num = 0 + >>> nipype_motioncorrect.run() # doctest: +SKIP + + Inputs:: + + inputnode.in_file + inputnode.ref_num + inputnode.in_bvec + + Outputs:: + + outputnode.motion_corrected + """ + + inputnode = pe.Node(interface = util.IdentityInterface(fields=["in_file", "ref_num","in_bvec"]), + name="inputnode") + + pipeline = pe.Workflow(name=name) + + split = pe.Node(fsl.Split(dimension='t'), name="split") + pick_ref = pe.Node(util.Select(), name="pick_ref") + coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1, dof=6), name = "coregistration", iterfield=["in_file"]) + rotate_bvecs = pe.Node( util.Function( input_names=["in_bvec", "in_matrix"], output_names=["out_file"], function=_rotate_bvecs ), name="rotate_b_matrix" ) + merge = pe.Node(fsl.Merge(dimension="t"), name="merge") + outputnode = pe.Node(interface = util.IdentityInterface(fields=["motion_corrected","out_bvec"]), + name="outputnode") + + pipeline.connect([ + (inputnode, split, [("in_file", "in_file")]) + ,(split, pick_ref, [("out_files", "inlist")]) + ,(inputnode, pick_ref, [("ref_num", "index")]) + ,(split, coregistration, [("out_files", "in_file")]) + ,(inputnode, rotate_bvecs, [("in_bvec","in_bvec")]) + ,(coregistration, rotate_bvecs, [("out_matrix_file","in_matrix")]) + ,(pick_ref, coregistration, [("out", "reference")]) + ,(coregistration, merge, [("out_file", "in_files")]) + ,(merge, outputnode, [("merged_file", "motion_corrected")]) + ,(rotate_bvecs, outputnode, [("out_file","out_bvec")]) + ]) + + return pipeline + + def create_eddy_correct_pipeline(name="eddy_correct"): """Creates a pipeline that replaces eddy_correct script in FSL. It takes a series of diffusion weighted images and linearly corregisters them to one @@ -197,30 +287,26 @@ def create_eddy_correct_pipeline(name="eddy_correct"): pipeline = pe.Workflow(name=name) split = pe.Node(fsl.Split(dimension='t'), name="split") - pipeline.connect([(inputnode, split, [("in_file", "in_file")])]) - pick_ref = pe.Node(util.Select(), name="pick_ref") - pipeline.connect([(split, pick_ref, [("out_files", "inlist")]), - (inputnode, pick_ref, [("ref_num", "index")])]) - - coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1), name = "coregistration", iterfield=["in_file"]) - pipeline.connect([(split, coregistration, [("out_files", "in_file")]), - (pick_ref, coregistration, [("out", "reference")])]) - + coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1 ), name = "coregistration", iterfield=["in_file"]) merge = pe.Node(fsl.Merge(dimension="t"), name="merge") - pipeline.connect([(coregistration, merge, [("out_file", "in_files")]) - ]) - outputnode = pe.Node(interface = util.IdentityInterface(fields=["eddy_corrected"]), name="outputnode") - pipeline.connect([(merge, outputnode, [("merged_file", "eddy_corrected")])]) - + pipeline.connect([ + (inputnode, split, [("in_file", "in_file")]) + ,(split, pick_ref, [("out_files", "inlist")]) + ,(inputnode, pick_ref, [("ref_num", "index")]) + ,(split, coregistration, [("out_files", "in_file")]) + ,(pick_ref, coregistration, [("out", "reference")]) + ,(coregistration, merge, [("out_file", "in_files")]) + ,(merge, outputnode, [("merged_file", "eddy_corrected")]) + ]) return pipeline -def create_epi_correct_pipeline(name="epi_correct", register_to_ref=False ): +def create_susceptibility_correct_pipeline(name="susceptibility_correct", register_to_ref=False ): """ Replaces the epidewarp.fsl script (http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl) - for EPI distortion correction with the fieldmap information in + for susceptibility distortion correction in EPI sequences with the fieldmap information in dMRI and fMRI (Jezzard et al., MRM 1995 ) using FSL's FUGUE. Example @@ -271,9 +357,6 @@ def create_epi_correct_pipeline(name="epi_correct", register_to_ref=False ): pipeline = pe.Workflow(name=name) - - matrix_file = os.path.abspath( './flirt.txt' ) - # Keep first frame from magnitude select_mag = pe.Node( interface=fsl.utils.ExtractROI(t_size=1,t_min=0), name="select_magnitude" ) @@ -372,6 +455,39 @@ def create_epi_correct_pipeline(name="epi_correct", register_to_ref=False ): return pipeline +def _rotate_bvecs( in_bvec, in_matrix ): + import os + import numpy as np + + name,fext = os.path.splitext( os.path.basename(in_bvec) ) + if fext == '.gz': name,_ = os.path.splitext(name) + out_file = os.path.abspath( './%s_rotated.bvec' % name ) + bvecs = np.loadtxt( in_bvec ) + new_bvecs = [ bvecs[:,0] ] + + for i,vol_matrix in enumerate(in_matrix[1::]): + bvec = np.matrix( bvecs[:,i] ) + rot = np.matrix(np.loadtxt( vol_matrix )[0:3,0:3]) + new_bvecs.append( (np.array( rot * bvec.T).T)[0] ) + np.savetxt( out_file, np.array(new_bvecs).T, fmt='%0.15f' ) + return out_file + +def _cat_logs( in_files ): + import shutil + import os + + name,fext = os.path.splitext( os.path.basename(in_files[0]) ) + if fext == '.gz': name,_ = os.path.splitext(name) + out_file = os.path.abspath( './%s_ecclog.log' % name ) + out_str = "" + with open( out_file, 'wb' ) as totallog: + for i,fname in enumerate(in_files): + totallog.write( "\n\npreprocessing %d\n" % i ) + with open(fname) as inlog: + for line in inlog: + totallog.write(line) + return out_file + def _compute_dwelltime( dwell_time=0.68, is_parallel=False, is_reverse_encoding=False ): if is_parallel: dwell_time*=0.5 From 25a09f388443569eda4196890d00ebbba1127eca Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 6 Mar 2013 13:37:08 +0100 Subject: [PATCH 08/11] Updated dMRI and fMRI pre-processing * ENH: modified parameter for parallel imaging, now it is possible to modify the acceleration factor * ENH: created new pipeline for full dMRI-EPI full distortions correction. Includes: motion, eddy currents and susceptibility correction --- nipype/workflows/dmri/fsl/dti.py | 131 ++++++++++++++++++++++--------- 1 file changed, 95 insertions(+), 36 deletions(-) diff --git a/nipype/workflows/dmri/fsl/dti.py b/nipype/workflows/dmri/fsl/dti.py index 2ac6a7403f..e061e2207c 100644 --- a/nipype/workflows/dmri/fsl/dti.py +++ b/nipype/workflows/dmri/fsl/dti.py @@ -12,36 +12,86 @@ def transpose(samples_over_fibres): a = a.reshape(-1,1) return a.T.tolist() -def create_dmri_preprocessing(name="dMRI_preprocessing"): +def create_dmri_preprocessing(name="dMRI_preprocessing", fieldmap_registration=False): """Creates a workflow that chains the necessary pipelines to - correct for motion, eddy currents, and susceptibility - artifacts in EPI dMRI sequences. + correct for motion, eddy currents, and susceptibility + artifacts in EPI dMRI sequences. Example ------- >>> nipype_dmri_preprocess = create_dmri_preprocessing("nipype_dmri_prep") - >>> nipype_dmri_preprocess.inputs.inputnode. = - >>> nipype_dmri_preprocess.inputs.inputnode. = - >>> nipype_dmri_preprocess.inputs.inputnode. = - >>> nipype_dmri_preprocess.inputs.inputnode. = - >>> nipype_dmri_preprocess.inputs.inputnode. = - >>> nipype_dmri_preprocess.inputs.inputnode. = - >>> nipype_dmri_preprocess.inputs.inputnode. = - >>> nipype_dmri_preprocess.inputs.inputnode. = + >>> nipype_dmri_preprocess.inputs.inputnode.in_file = 'diffusion.nii' + >>> nipype_dmri_preprocess.inputs.inputnode.in_bvec = 'diffusion.bvec' + >>> nipype_dmri_preprocess.inputs.inputnode.ref_num = 0 + >>> nipype_dmri_preprocess.inputs.inputnode.fieldmap_mag = 'magnitude.nii' + >>> nipype_dmri_preprocess.inputs.inputnode.fieldmap_pha = 'phase.nii' + >>> nipype_dmri_preprocess.inputs.inputnode.te_diff = 2.46 + >>> nipype_dmri_preprocess.inputs.inputnode.epi_echospacing = 0.77 + >>> nipype_dmri_preprocess.inputs.inputnode.epi_rev_encoding = False + >>> nipype_dmri_preprocess.inputs.inputnode.pi_accel_factor = True >>> nipype_dmri_preprocess.run() # doctest: +SKIP + Inputs:: - inputnode. + inputnode.in_file + inputnode.in_bvec + inputnode.ref_num + inputnode.fieldmap_mag + inputnode.fieldmap_pha + inputnode.te_diff + inputnode.epi_echospacing + inputnode.epi_rev_encoding + inputnode.pi_accel_factor + inputnode.vsm_sigma + Outputs:: outputnode.dmri_corrected + outputnode.bvec_rotated + + + Optional arguments:: + + fieldmap_registration - True if registration to fieldmap should be done (default False) + """ pipeline = pe.Workflow(name=name) + + inputnode = pe.Node(interface = util.IdentityInterface( + fields=['in_file', 'in_bvec','ref_num','fieldmap_mag', + 'fieldmap_pha','te_diff','epi_echospacing', + 'epi_rev_encoding','pi_accel_factor','vsm_sigma']), + name='inputnode') + + outputnode= pe.Node(interface = util.IdentityInterface( + fields=['dmri_corrected','bvec_rotated']), + name='outputnode') + + motion = create_motion_correct_pipeline() + eddy = create_eddy_correct_pipeline() + susceptibility = create_susceptibility_correct_pipeline(fieldmap_registration=fieldmap_registration) + + pipeline.connect([ + (inputnode, motion, [('in_file','inputnode.in_file'),('in_bvec','inputnode.in_bvec'),('ref_num','inputnode.ref_num')]) + ,(inputnode, eddy, [('ref_num','inputnode.ref_num')]) + ,(motion, eddy, [('outputnode.motion_corrected','inputnode.in_file')]) + ,(eddy, susceptibility, [('outputnode.eddy_corrected','inputnode.in_file')]) + ,(inputnode, susceptibility, [('ref_num','inputnode.ref_num') + ,('fieldmap_mag','inputnode.fieldmap_mag') + ,('fieldmap_pha','inputnode.fieldmap_pha') + ,('te_diff','inputnode.te_diff') + ,('epi_echospacing','inputnode.epi_echospacing') + ,('epi_rev_encoding','inputnode.epi_rev_encoding') + ,('pi_accel_factor','inputnode.pi_accel_factor') + ,('vsm_sigma','inputnode.vsm_sigma') ]) + ,(motion, outputnode, [('outputnode.out_bvec', 'bvec_rotated')]) + ,(susceptibility, outputnode, [('outputnode.epi_corrected','dmri_corrected')]) + ]) return pipeline @@ -206,7 +256,8 @@ def create_bedpostx_pipeline(name="bedpostx"): def create_motion_correct_pipeline(name="motion_correct"): """Creates a pipeline that corrects for motion artifact in dMRI sequences. It takes a series of diffusion weighted images and rigidly corregisters - them to one reference image. Finally, the b-matrix is rotated accordingly, + them to one reference image. Finally, the b-matrix is rotated accordingly + (Leemans et al. 2009 - http://www.ncbi.nlm.nih.gov/pubmed/19319973), making use of the rotation matrix obtained by FLIRT. Example @@ -227,6 +278,8 @@ def create_motion_correct_pipeline(name="motion_correct"): Outputs:: outputnode.motion_corrected + outputnode.out_bvec + """ inputnode = pe.Node(interface = util.IdentityInterface(fields=["in_file", "ref_num","in_bvec"]), @@ -260,8 +313,9 @@ def create_motion_correct_pipeline(name="motion_correct"): def create_eddy_correct_pipeline(name="eddy_correct"): """Creates a pipeline that replaces eddy_correct script in FSL. It takes a - series of diffusion weighted images and linearly corregisters them to one - reference image. + series of diffusion weighted images and linearly co-registers them to one + reference image. No rotation of the B-matrix is performed, so this pipeline + should be executed after the motion correction pipeline. Example ------- @@ -288,7 +342,7 @@ def create_eddy_correct_pipeline(name="eddy_correct"): split = pe.Node(fsl.Split(dimension='t'), name="split") pick_ref = pe.Node(util.Select(), name="pick_ref") - coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1 ), name = "coregistration", iterfield=["in_file"]) + coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1, dof=12 ), name = "coregistration", iterfield=["in_file"]) merge = pe.Node(fsl.Merge(dimension="t"), name="merge") outputnode = pe.Node(interface = util.IdentityInterface(fields=["eddy_corrected"]), name="outputnode") @@ -304,23 +358,24 @@ def create_eddy_correct_pipeline(name="eddy_correct"): ]) return pipeline -def create_susceptibility_correct_pipeline(name="susceptibility_correct", register_to_ref=False ): +def create_susceptibility_correct_pipeline(name="susceptibility_correct", fieldmap_registration=False ): """ Replaces the epidewarp.fsl script (http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl) - for susceptibility distortion correction in EPI sequences with the fieldmap information in - dMRI and fMRI (Jezzard et al., MRM 1995 ) using FSL's FUGUE. + for susceptibility distortion correction of dMRI & fMRI acquired with EPI sequences and the fieldmap + information (Jezzard et al., 1995) using FSL's FUGUE. The registration to the (warped) fieldmap + (strictly following the original script) is available using fieldmap_registration=True. Example ------- - >>> nipype_epicorrect = create_epi_correct_pipeline("nipype_epicorrect") - >>> nipype_epicorrect.inputs.inputnode.in_file = 'epi_data.nii' + >>> nipype_epicorrect = create_susceptibility_correct_pipeline("nipype_epicorrect", fieldmap_registration=False) + >>> nipype_epicorrect.inputs.inputnode.in_file = 'diffusion.nii' >>> nipype_epicorrect.inputs.inputnode.fieldmap_mag = 'magnitude.nii' >>> nipype_epicorrect.inputs.inputnode.fieldmap_pha = 'phase.nii' >>> nipype_epicorrect.inputs.inputnode.te_diff = 2.46 - >>> nipype_epicorrect.inputs.inputnode.epi_echospacing = 0.51 + >>> nipype_epicorrect.inputs.inputnode.epi_echospacing = 0.77 >>> nipype_epicorrect.inputs.inputnode.epi_rev_encoding = False - >>> nipype_epicorrect.inputs.inputnode.ref_volume = 0 - >>> nipype_epicorrect.inputs.inputnode.epi_parallel = True + >>> nipype_epicorrect.inputs.inputnode.ref_num = 0 + >>> nipype_epicorrect.inputs.inputnode.pi_accel_factor = 1.0 >>> nipype_epicorrect.run() # doctest: +SKIP Inputs:: @@ -332,14 +387,20 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", regist inputnode.epi_echospacing - The echo spacing (aka dwell time) in the EPI sequence inputnode.epi_ph_encoding_dir - The phase encoding direction in EPI acquisition (default y) inputnode.epi_rev_encoding - True if it is acquired with reverse encoding - inputnode.epi_parallel - True if EPI was acquired in a parallel imaging scheme + inputnode.pi_accel_factor - Acceleration factor used for EPI parallel imaging (GRAPPA) inputnode.vsm_sigma - Sigma value of the gaussian smoothing filter applied to the vsm (voxel shift map) - inputnode.ref_volume - The reference volume (B=0 in dMRI or a central frame in fMRI) + inputnode.ref_num - The reference volume (B=0 in dMRI or a central frame in fMRI) Outputs:: outputnode.epi_corrected + + + Optional arguments:: + + fieldmap_registration - True if registration to fieldmap should be done (default False) + """ inputnode = pe.Node(interface = util.IdentityInterface(fields=["in_file", @@ -349,9 +410,9 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", regist "epi_echospacing", "epi_ph_encoding_dir", "epi_rev_encoding", - "epi_parallel", + "pi_accel_factor", "vsm_sigma", - "ref_volume", + "ref_num", "unwarp_direction" ]), name="inputnode") @@ -365,7 +426,7 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", regist mask_mag_dil = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_dilate_mask), name='mask_dilate' ) # Compute dwell time - dwell_time = pe.Node( interface=util.Function( input_names=["dwell_time","is_parallel","is_reverse_encoding"], output_names=["dwell_time"], function=_compute_dwelltime), name='dwell_time') + dwell_time = pe.Node( interface=util.Function( input_names=["dwell_time","pi_factor","is_reverse_encoding"], output_names=["dwell_time"], function=_compute_dwelltime), name='dwell_time') # Normalize phase diff to be [-pi, pi) norm_pha = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_prepare_phasediff ), name='normalize_phasediff') @@ -392,7 +453,7 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", regist pipeline.connect([ - (inputnode, dwell_time, [('epi_echospacing','dwell_time'), ('epi_parallel','is_parallel'),('epi_rev_encoding','is_reverse_encoding') ]) + (inputnode, dwell_time, [('epi_echospacing','dwell_time'), ('pi_accel_factor','pi_factor'),('epi_rev_encoding','is_reverse_encoding') ]) ,(inputnode, select_mag, [('fieldmap_mag','in_file')] ) ,(inputnode, norm_pha, [('fieldmap_pha','in_file')] ) ,(select_mag, mask_mag, [('roi_file', 'in_file')] ) @@ -414,8 +475,7 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", regist ,(dwi_merge, outputnode, [('merged_file','epi_corrected')]) ]) - if register_to_ref: - # register to ref volume + if fieldmap_registration: """ Register magfw to example epi. There are some parameters here that may need to be tweaked. Should probably strip the mag Pre-condition: forward warp the mag in order to reg with func. What does mask do here? """ @@ -431,7 +491,7 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", regist msk_applyxfm = pe.Node( interface=fsl.ApplyXfm(), name='msk_apply_xfm') pipeline.connect([ - (inputnode, select_epi, [('in_file','in_file'), ('ref_volume','t_min')] ) + (inputnode, select_epi, [('in_file','in_file'), ('ref_num','t_min')] ) ,(select_epi, vsm_reg, [('roi_file','reference')]) ,(vsm, vsm_fwd, [('shift_out_file','shift_in_file')]) ,(mask_mag_dil, vsm_fwd, [('out_file','mask_file')] ) @@ -488,9 +548,8 @@ def _cat_logs( in_files ): totallog.write(line) return out_file -def _compute_dwelltime( dwell_time=0.68, is_parallel=False, is_reverse_encoding=False ): - if is_parallel: - dwell_time*=0.5 +def _compute_dwelltime( dwell_time=0.68, pi_factor=1.0, is_reverse_encoding=False ): + dwell_time*=(1.0/pi_factor) if is_reverse_encoding: dwell_time*=-1.0 From 4fea36f7e479d324606db6a14c75b33735e7294f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 6 Mar 2013 14:49:00 +0100 Subject: [PATCH 09/11] Updated changelog --- CHANGES | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/CHANGES b/CHANGES index 34c2dcfead..1fd57c3645 100644 --- a/CHANGES +++ b/CHANGES @@ -3,9 +3,10 @@ Next release * ENH: New interfaces: nipy.Trim * ENH: Allow control over terminal output for commandline interfaces -* ENH: New workflows: susceptibility correction for EPI imaging based on epidewarp.fsl - updated to support both dMRI and fMRI, and new parameters. -* ENH: Minor improvements to FSL's FUGUE interface +* ENH: New workflows for dMRI and fMRI pre-processing: added motion artifact correction + with rotation of the B-matrix, and susceptibility correction for EPI imaging using + fieldmaps. Updated eddy_correct pipeline to support both dMRI and fMRI, and new parameters. +* ENH: Minor improvements to FSL's FUGUE and FLIRT interfaces Release 0.7.0 (Dec 18, 2012) ============================ From 2f6338d44814bde525f4b3856e85e53478c97577 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 6 Mar 2013 15:12:21 +0100 Subject: [PATCH 10/11] Removed accidentally added files (#530) Removed files regarding MBIS tool --- nipype/interfaces/mbis/__init__.py | 2 - nipype/interfaces/mbis/postproc.py | 236 ------------------- nipype/interfaces/mbis/segmentation.py | 303 ------------------------- nipype/interfaces/setup.py | 1 - nipype/workflows/dmri/fsl/__init__.py | 2 +- 5 files changed, 1 insertion(+), 543 deletions(-) delete mode 100644 nipype/interfaces/mbis/__init__.py delete mode 100644 nipype/interfaces/mbis/postproc.py delete mode 100644 nipype/interfaces/mbis/segmentation.py diff --git a/nipype/interfaces/mbis/__init__.py b/nipype/interfaces/mbis/__init__.py deleted file mode 100644 index 931c115a0b..0000000000 --- a/nipype/interfaces/mbis/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ - -from .segmentation import (MBIS,PVMerge,MAP) diff --git a/nipype/interfaces/mbis/postproc.py b/nipype/interfaces/mbis/postproc.py deleted file mode 100644 index ef47c30098..0000000000 --- a/nipype/interfaces/mbis/postproc.py +++ /dev/null @@ -1,236 +0,0 @@ - -def loadParameters( fname ): - import csv - import numpy as np - - data = [] - with open( fname ) as csvfile: - for row in csvfile: - row = row.replace("[","").replace("]","").replace( " ", "") - data.append( [ float(val) for val in row.split(',') ] ) - - nparam = len( data[0] ) - nclasses = int( 0.5 * (nparam*4+1)**0.5 - 0.5 ) - # Load parameters - mu = [] - sigma = [] - - for row in data: - mu.append( np.matrix( [ float(val) for val in row[0:nclasses] ] ) ) - sigma.append( np.matrix( np.reshape( [ float(val) for val in row[nclasses:] ], (nclasses,nclasses) ) ) ) - - # Return arrays - return ( mu, sigma ) - -def distancePV ( sample, mask, params_tissue1, params_tissue2, distance ): - from scipy.spatial.distance import mahalanobis,euclidean - import numpy as np - - # Direction vector between pure tissues - d_vect = np.ravel(params_tissue2[0] - params_tissue1[0]).T - mu1 = np.ravel(params_tissue1[0]) - mu2 = np.ravel(params_tissue2[0]) - SI1 = params_tissue1[1].getI() - SI2 = params_tissue2[1].getI() - - if distance=='mahalanobis': - norm = np.array( [ 1/(1+ mahalanobis(pix,mu2,SI2)/ mahalanobis(pix,mu1,SI1)) for pix in sample[mask==1] ] ) - elif distance=='dummy': - norm = mask*0.5 - else: - norm = np.array( [ 1/(1+ euclidean(pix,mu2)/ euclidean(pix,mu1)) for pix in sample[mask==1] ] ) - result = np.zeros( np.shape( mask ) ) - result[mask==1] = norm - return result - -def computeMAP( tpms, normalize=True ): - import numpy as np - normalizer = np.sum( tpms, axis=0 ) - if normalize: - for tpm in tpms: - tpm[normalizer>0]/= normalizer[normalizer>0] - - out_seg = np.zeros( shape=tpms[0].shape, dtype=np.uint8 ) - out_seg[normalizer>0] = np.argmax( tpms, axis=0 )[normalizer>0] - out_seg[normalizer>0] += 1 - - return out_seg, tpms - -def fusePV2( in_files, in_maps, parameters, pt_list=( (0,1,2), (2,3), (4,) ), distance='mahalanobis', reorder=True, prefix='./' ): - import nibabel as nib - import numpy as np - import os - import sys - import collections - from scipy.spatial.distance import mahalanobis,euclidean - - nmaps = len( in_maps ) - nclasses = np.shape( parameters )[1] - assert nmaps == nclasses - - idxs = np.array(range( 0, nclasses )) - corder = range( 0, nclasses ) - npts = len( pt_list ) - - pv_task = np.array([ len(t)>1 for t in pt_list ]) # Check what tissues have pv actions - - if not pv_task.any(): - assert( nmaps == npts ) - return in_maps - - # Load probability maps - initmaps = [ nib.load(f).get_data() for f in in_maps ] - - # If reorder is True, find unordered tissue signatures - if reorder: - means = parameters[0] - firstmeans = np.array( [ np.ravel(val)[0] for val in means ] ) - idxs = np.argsort(firstmeans) - corder = np.take(corder, idxs) - - # Load images - channels = [ nib.load(c) for c in in_files ] - - # Prepare sample - data_shape = np.shape( channels[0] ) - data_sample = [ channel.get_data().reshape(-1) for channel in channels ] - data_sample = np.swapaxes( data_sample, 0, 1) - - t_flatten = y=collections.Counter( [element for tupl in pt_list for element in tupl] ) - pt_mapping = [ [ i, [] ] for i in range(0,npts) ] - #pv_mapping = [ [ np.where(idxs==i)[0] , [], i ] for i in t_flatten if t_flatten[i]>1 ] - pv_mapping = [ [ idxs[i] , [] ] for i in t_flatten if t_flatten[i]>1 ] - pt_maps = [ np.zeros( shape=data_shape, dtype=float) for i in range(0,npts)] - for pt_tuple,out_i in zip( pt_list, range(0,npts) ): - for t in pt_tuple: - in_i = idxs[t] - if t_flatten[t]==1: - pt_mapping[out_i][1].append( in_i ) - pt_maps[out_i] += initmaps[ in_i ] - else: - pvm_id = pv_mapping[:][0].index(in_i) - pv_mapping[pvm_id][1].append( out_i ) - - for m in pv_mapping: - in_id = m[0] - neighs = [ pt_mapping[val][1] for val in m[1] ] - dist = [] - mask = np.zeros(data_shape) - mask[initmaps[in_id]>0.001] = 1 - mu = [] - SI = [] - out_ids = [] - for n,l in zip(neighs, m[1][:]): - out_ids.append(l) - mu.append( parameters[0][n[-1]] ) - SI.append( parameters[1][n[-1]].getI() ) - - diff = np.zeros( shape=np.shape(mask.reshape(-1)), dtype=float) - norm = np.array([1/(1+ mahalanobis(pix,mu[1],SI[1])/ mahalanobis(pix,mu[0],SI[0])) for pix in data_sample[mask.reshape(-1)==1]]) - diff[mask.reshape(-1)==1] = norm - pt_maps[out_ids[0]]+= diff.reshape(data_shape) * initmaps[in_id] - pt_maps[out_ids[1]]+= (1-diff.reshape(data_shape)) * initmaps[in_id] - - # Normalize tpms and compute MAP - normalizer = np.sum( pt_maps, axis=0 ) - for pt_map in pt_maps: - pt_map[normalizer>0]/= normalizer[normalizer>0] - - # Generate output names - ppmnames = [ '%s_pvseg%02d.nii.gz' % ( prefix, i ) for i in range(0,len(pt_maps)) ] - - # Save - ref = nib.load( in_maps[0] ) - ref = channels[0] - for i in range(0,len(pt_maps)): - nii = nib.Nifti1Image( np.reshape( pt_maps[i], data_shape), ref.get_affine(), ref.get_header() ) - nib.save( nii, ppmnames[i] ) - - return ppmnames - - -def fusePV( in_files, in_maps, parameters, pt_list=[ 0, 2, 4 ], distance='mahalanobis', reorder=True, prefix='./' ): - import nibabel as nib - import numpy as np - import os - - nmaps = len( in_maps ) - nclasses = np.shape( parameters )[1] - assert nmaps == nclasses - - corder = range( 0, nclasses ) - npts = len( pt_list ) - - if npts == nmaps: - return in_maps - - pt_list = np.sort( pt_list ) - - # Load probability maps - initmaps = [ nib.load(f) for f in in_maps ] - - # If reorder is True, find unordered tissue signatures - if reorder: - means = parameters[0] - firstmeans = np.array( [ np.ravel(val)[0] for val in means ] ) - m = np.argsort(firstmeans) - corder = np.take(corder, m) - new_idx = np.take( corder, pt_list ) - - # Load images - channels = [ nib.load(c) for c in in_files ] - - # Prepare sample - data_shape = np.shape( channels[0] ) - data_sample = [ channel.get_data().reshape(-1) for channel in channels ] - data_sample = np.swapaxes( data_sample, 0, 1) - - - # Split between pv (partial volume) and pt (pure tissue) maps - pt_niis = [] - pv_niis = [] - pt_param = [] - - for tmap,i in zip(initmaps,range(0,nclasses)): - idx = np.where( new_idx==i )[0] - if len(idx) == 1: - pt_niis.append( tmap ) - pt_param.append( [ parameters[0][i], parameters[1][i] ] ) - else: - pv_niis.append( tmap ) - - # Compute the steps required - steps = [ val-pt_list[i-1]-1 for val,i in zip( pt_list[1:], range(1,len(pt_list[1:])+1 ) ) ] - - # Extract data and initialize normalizer - pt_maps = [ m.get_data().reshape(-1) for m in pt_niis ] - pv_maps = [ m.get_data().reshape(-1) for m in pv_niis ] - - # Process - for pt_map,i in zip( pt_maps[:-1],range(0,len(pt_maps[:-1])) ): - curr_steps = steps[i] - for step in range(1,curr_steps+1): - pv_idx = (step-1)+i - pv_map = pv_maps[pv_idx] - mask = np.zeros( np.shape( pt_map ) ) - mask[pv_map>0.001] = 1 - - if not step == curr_steps: # Direct addition of the pv map to the last pure tissue map - pt_map+= pv_map - else: # Split pv fraction proportionally to the distance to a contiguous pure tissue - dist = distancePV( data_sample, mask, pt_param[i], pt_param[i+1], distance ) - pt_map+= dist * pv_map - pt_maps[i+1]+= (1-dist) * pv_map - - # Normalize tpms and compute MAP - normalizer = np.sum( pt_maps, axis=0 ) - - # Generate output names - ppmnames = [ '%s_pvseg%02d.nii.gz' % ( prefix[:-6], i ) for i in range(0,len(pt_niis)) ] - - # Save - for i in range(0,len(pt_niis)): - nii = nib.Nifti1Image( np.reshape( pt_maps[i], data_shape) , pt_niis[0].get_affine(), pt_niis[0].get_header() ) - nib.save( nii, ppmnames[i] ) - - return ppmnames diff --git a/nipype/interfaces/mbis/segmentation.py b/nipype/interfaces/mbis/segmentation.py deleted file mode 100644 index d60d2d9c2a..0000000000 --- a/nipype/interfaces/mbis/segmentation.py +++ /dev/null @@ -1,303 +0,0 @@ -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -""" MBIS nipype interface definition - Based on FAST interface definition - - 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, os.path as op -import warnings -import numpy as np -import nibabel as nib -from nipype.interfaces.base import (TraitedSpec, File, InputMultiPath, - OutputMultiPath, Undefined, traits, - isdefined, OutputMultiPath, - CommandLineInputSpec, CommandLine, - BaseInterface, BaseInterfaceInputSpec, - traits ) -from nipype.utils.filemanip import split_filename,fname_presuffix -import postproc as pp -#import csv - -warn = warnings.warn -warnings.filterwarnings('always', category=UserWarning) - - - -class MBISInputSpec( CommandLineInputSpec ): - in_files = InputMultiPath( File(exists=True), copyfile=False, - desc='image, or multi-channel set of images, ' \ - 'to be segmented', - argstr='-C %s', position=-1, mandatory=True ) - - mask = File(exists=True, desc='binary mask file', argstr='-x %s' ) - - mask_auto = traits.Bool( desc='channels are implicitly masked', argstr='-M' ) - - out_prefix = File('outpath', desc='base name of output files', - argstr='-o %s', genfile=True) # uses in_file name as basename if none given - - initialization = traits.Enum( 'none','em','km','km+em', argstr='-I %s',desc='Initialization mode' ) - - number_classes = traits.Range(low=2, high=10, argstr='-n %d', - desc='number of tissue-type classes', value=3) - - output_steps = traits.Bool(desc='output intermediate steps', - argstr='--output-steps') - - - output_biasfield = traits.Bool(desc='output estimated bias field', - argstr='--bias-output') - - output_biascorrected = traits.Bool(desc='output restored image ' \ - '(bias-corrected image)', - argstr='--bias-corrected-output') - - output_stats = File( 'outcsvfile', desc='output file containing mixture parameters', argstr='--output-stats %s' ) - - probability_maps = traits.Bool(desc='outputs a separate binary image for each ' \ - 'tissue type', - argstr='-g' ) - - priors = InputMultiPath(File(exist=True), desc='initialize with prior images', - argstr='-P %s', minlen=3, maxlen=10) - - no_bias = traits.Bool(desc='do not remove bias field', - argstr='--bias-skip' ) - - no_normalize = traits.Bool( desc='do not normalize input channels', argstr='-N', - value=False ) - - em_iters = traits.Range(low=1, high=50, value=3, - desc='number of EM iterations', - argstr='--em-iterations %d') - mrf_iters = traits.Range(low=1, high=10, - desc='number of MRF iterations', - argstr='--mrf-iterations %d') - - mrf_lambda = traits.Range(low=0.01, high=3.0, - desc='MRF lambda parameter (segmentation spatial smoothness)', - argstr='-l %.3f') - - manual_init = File(exists=True, desc='Filename containing intensities', - argstr='-f %s') - manual_init_means = File(exists=True, desc='Filename containing intensities', - argstr='--init-means %s') - -# no_pve = traits.Bool(desc='turn off PVE (partial volume estimation)', -# argstr='--nopve') - -# use_priors = traits.Bool(desc='use priors throughout', -# argstr='-P') # must also set -a!, -# # mutually inclusive?? -# # No, conditional -# # mandatory... need to -# # figure out how to -# # handle with traits. -# verbose = traits.Bool(desc='switch on diagnostic messages', -# argstr='-v') - - - -class MBISOutputSpec( TraitedSpec ): - out_classified = File(desc='path/name of binary segmented volume file' \ - ' one val for each class _mrf') - out_parameters = traits.File(desc='csv file with tissue purameters') - bias_field = OutputMultiPath(File(desc='Estimated bias field _bias')) - probability_maps = OutputMultiPath(File(desc='filenames, one for each class, for each ' \ - 'input, mrf_x')) - restored_image = OutputMultiPath(File(desc='restored images (one for each input image) ' \ - 'named according to the input images _corrected_chXX')) - - normalized_inputs = OutputMultiPath( File( desc='normalized input channels' ) ) - - -class MBIS(CommandLine): - """ Use MBIS for segmentation and bias correction. - - Examples - -------- - >>> from nipype.interfaces import mbis - >>> from nipype.testing import example_data - - Assign options through the ``inputs`` attribute: - - >>> mbisr = mbis.MBIS() - >>> mbisr.inputs.in_files = example_data('structural.nii') - >>> out = mbisr.run() #doctest: +SKIP - - """ - _cmd = 'brain_seg' - input_spec = MBISInputSpec - output_spec = MBISOutputSpec - - def _list_outputs(self): - outputs = self.output_spec().get() - if not isdefined(self.inputs.number_classes): - nclasses = 3 - else: - nclasses = self.inputs.number_classes - - outprefix = self._getdefaultprefix() - - outputs['out_classified'] = '%s_mrf.nii.gz' % outprefix - - if self.inputs.probability_maps: - outputs['probability_maps'] = [] - for i in range(1,nclasses+1): - outputs['probability_maps'].append('%s_mrf_seg%02d.nii.gz' % (outprefix,i) ) - - if (not (isdefined(self.inputs.no_bias))) and self.inputs.output_biascorrected: - outputs['restored_image'] = [] - for val, f in enumerate(self.inputs.in_files): - # image numbering is 1-based - outputs['restored_image'].append('%s_corrected_ch%02d.nii.gz' % (outprefix,val+1) ) - - if (not (isdefined(self.inputs.no_bias))) and self.inputs.output_biasfield: - outputs['bias_field'] = [] - for val, f in enumerate(self.inputs.in_files): - # image numbering is 1-based - outputs['bias_field'].append('%s_bias_ch%02d.nii.gz' % (outprefix,val) ) - - if ((not self.inputs.no_normalize) and self.inputs.output_steps ): - outputs['normalized_inputs' ] = [] - for val, f in enumerate(self.inputs.in_files): - outputs['normalized_inputs'].append('%s_normin_%0d.nii.gz' % (outprefix,val) ) - - if isdefined(self.inputs.output_stats): - fname= outprefix + '_stats_final' + self.inputs.output_stats - outputs['out_parameters'] = fname -# outputs['out_parameters'] = np.loadtxt( fname, delimiter='[],' ) -# with open( self.inputs.output_stats ) as csvfile: -# dataReader = csv.reader( csvfile ) -# outputs['out_parameters'] = np.array( [ [ row ] for row in dataReader ] ) -# csvfile.close() - - return outputs - - def _getdefaultprefix( self, name='mbis' ): - if not isdefined(self.inputs.out_prefix): - return os.path.abspath( os.path.join( os.getcwd(), name ) ) - else: - return self.inputs.out_prefix - - def _gen_filename(self,name): - if name == 'out_prefix': - return self._getdefaultprefix() - return None - - - def _gen_fname(self, prefix, suffix=None, ext='.nii.gz', cwd=None): - """Generate a filename based on the given parameters. - - The filename will take the form: preffix. - - Parameters - ---------- - prefix : str - Filename to base the new filename on. - - suffix : str - Suffix to add to the `basename`. (defaults is '' ) - - ext : str - Desired extension (default is nii.gz) - - Returns - ------- - fname : str - New filename based on given parameters. - - """ - if ((prefix == '') or (prefix is None) ): - prefix = './' - if suffix is None: - suffix = '' - if cwd is None: - cwd = os.getcwd() - - suffix = ''.join((suffix,ext)) - fname = fname_presuffix(prefix, suffix=suffix, use_ext=False, newpath=cwd ) - return fname - - -class PVMergeInputSpec(BaseInterfaceInputSpec): - in_files = InputMultiPath( File(exists=True), copyfile=False, desc='image, or multi-channel set of images,' \ - 'used for segmentation', mandatory=True ) - in_maps = InputMultiPath( File(exists=True), copyfile=False, desc='Resulting tissue probability maps', mandatory=True ) - in_labelling = File( exists=True, desc='Hard labelling resulting from segmentation', mandatory=True) - parameters = File( exists=True, desc='CSV file containing parameters for all the classes', mandatory=True) - pure_tissues = traits.Tuple( value=( (0,), (1,), (2,) ), minlen=1, desc='Identifiers of the pure tissue classes' ) - dist = traits.String( value="euclidean", desc="Distance definition to be used" ) - reorder = traits.Bool( value=True, desc='Reorder classes if the classification is not ordered by first contrast means' ) - -class PVMergeOutputSpec(TraitedSpec): - out_files = OutputMultiPath( File( desc='filenames, one for each pure tissue according to prefix_msegXX.nii.gz' ) ) - -class PVMerge(BaseInterface): - input_spec = PVMergeInputSpec - output_spec = PVMergeOutputSpec - _ppmnames = [] - - def _run_interface(self, runtime): - params = pp.loadParameters(self.inputs.parameters) - path, base, _ = split_filename( self.inputs.in_maps[0] ) - - self._ppmnames = pp.fusePV2( self.inputs.in_files, - self.inputs.in_maps, - params, - self.inputs.pure_tissues, - self.inputs.dist, - self.inputs.reorder, - os.path.abspath( os.path.join( os.getcwd(), base ) ) ) - - return runtime - - def _list_outputs(self): - outputs = self._outputs().get() - outputs["out_files"] = self._ppmnames - return outputs - -class MAPInputSpec(BaseInterfaceInputSpec): - in_files = InputMultiPath( File(exists=True), copyfile=False, desc='Input tissue probability maps', mandatory=True ) - normalize = traits.Bool( desc='apply normalization to input tpms', value=True ) - prefix = traits.String( desc='prefix path and basename for output', value='' ) - -class MAPOutputSpec(TraitedSpec): - out_tpms = OutputMultiPath( File( desc='Normalized tpms' ) ) - out_file = File( desc='output labelling after MAP' ) - -class MAP(BaseInterface): - input_spec = MAPInputSpec - output_spec = MAPOutputSpec - - def _run_interface(self, runtime): - import nibabel as nib - import numpy as np - path, base, _ = split_filename( self.inputs.in_files[0] ) - prefix=self.inputs.prefix - if not isdefined(self.inputs.prefix): - prefix = os.getcwd() + "/"+ base - self._out_fname = '%s_map.nii.gz' % prefix - assert( len(self.inputs.in_files) > 1 ) - tpm_data = np.array( [ nib.load(tpm).get_data() for tpm in self.inputs.in_files ] ) - (seg_file,tpms_norm) = pp.computeMAP( tpm_data, self.inputs.normalize ) - ref = nib.load(self.inputs.in_files[0] ) - nib.save( nib.Nifti1Image( seg_file, ref.get_affine(), ref.get_header() ), self._out_fname ) - if self.inputs.normalize: - self._tpmnames = [ '%s_tpm%02d.nii.gz' % (prefix,i) for i in range(1,len(self.inputs.in_files)+1) ] - for tpm,name in zip(tpms_norm,self._tpmnames): - nib.save( nib.Nifti1Image( tpm, ref.get_affine(), ref.get_header() ), name ) - return runtime - - def _list_outputs(self): - outputs = self._outputs().get() - if self.inputs.normalize: - outputs["out_tpms"] = self._tpmnames - outputs["out_file"] = self._out_fname - return outputs diff --git a/nipype/interfaces/setup.py b/nipype/interfaces/setup.py index 430947e3de..b4fcb6dc36 100644 --- a/nipype/interfaces/setup.py +++ b/nipype/interfaces/setup.py @@ -16,7 +16,6 @@ def configuration(parent_package='', top_path=None): config.add_subpackage('dipy') config.add_subpackage('freesurfer') config.add_subpackage('fsl') - config.add_subpackage('mbis') config.add_subpackage('mne') config.add_subpackage('mrtrix') config.add_subpackage('nipy') diff --git a/nipype/workflows/dmri/fsl/__init__.py b/nipype/workflows/dmri/fsl/__init__.py index db5679e7da..916ba0d7d8 100644 --- a/nipype/workflows/dmri/fsl/__init__.py +++ b/nipype/workflows/dmri/fsl/__init__.py @@ -1,4 +1,4 @@ -from dti import create_bedpostx_pipeline, create_eddy_correct_pipeline, create_susceptibility_correct_pipeline +from dti import create_bedpostx_pipeline, create_eddy_correct_pipeline, create_susceptibility_correct_pipeline, create_dmri_preprocessing from tbss import (create_tbss_1_preproc, create_tbss_2_reg, create_tbss_3_postreg, create_tbss_4_prestats, create_tbss_all, create_tbss_non_FA) From 9325696b75409723cfcb48d509eb992a0a649928 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 6 Mar 2013 16:17:09 +0100 Subject: [PATCH 11/11] Fixed some details for the PR #530 * Improved documentation, and added warning about bvec rotation * Improved FSL FLIRT interface to support spline interpolation * Added spline interpolation wherever resampling is required * Syntax checked with autopep8 --- nipype/interfaces/fsl/preprocess.py | 353 ++++++++++-------- nipype/workflows/dmri/fsl/dti.py | 536 +++++++++++++++------------- 2 files changed, 474 insertions(+), 415 deletions(-) diff --git a/nipype/interfaces/fsl/preprocess.py b/nipype/interfaces/fsl/preprocess.py index 362818010e..c6c03ac06a 100644 --- a/nipype/interfaces/fsl/preprocess.py +++ b/nipype/interfaces/fsl/preprocess.py @@ -11,7 +11,8 @@ >>> os.chdir(datadir) """ -import os, os.path as op +import os +import os.path as op import warnings import numpy as np @@ -33,10 +34,10 @@ class BETInputSpec(FSLCommandInputSpec): # We use position args here as list indices - so a negative number # will put something on the end in_file = File(exists=True, - desc='input file to skull strip', - argstr='%s', position=0, mandatory=True) + desc='input file to skull strip', + argstr='%s', position=0, mandatory=True) out_file = File(desc='name of output skull stripped image', - argstr='%s', position=1, genfile=True, hash_files=False) + argstr='%s', position=1, genfile=True, hash_files=False) outline = traits.Bool(desc='create surface outline image', argstr='-o') mask = traits.Bool(desc='create binary mask image', @@ -44,19 +45,19 @@ class BETInputSpec(FSLCommandInputSpec): skull = traits.Bool(desc='create skull image', argstr='-s') no_output = traits.Bool(argstr='-n', - desc="Don't generate segmented output") + desc="Don't generate segmented output") frac = traits.Float(desc='fractional intensity threshold', argstr='-f %.2f') vertical_gradient = traits.Float(argstr='-g %.2f', - desc='vertical gradient in fractional intensity ' \ - 'threshold (-1, 1)') + desc='vertical gradient in fractional intensity ' + 'threshold (-1, 1)') radius = traits.Int(argstr='-r %d', units='mm', desc="head radius") center = traits.List(traits.Int, desc='center of gravity in voxels', argstr='-c %s', minlen=0, maxlen=3, units='voxels') threshold = traits.Bool(argstr='-t', - desc="apply thresholding to segmented brain image and mask") + desc="apply thresholding to segmented brain image and mask") mesh = traits.Bool(argstr='-e', desc="generate a vtk mesh brain surface") # the remaining 'options' are more like modes (mutually exclusive) that @@ -66,20 +67,20 @@ class BETInputSpec(FSLCommandInputSpec): # supported _xor_inputs = ('functional', 'reduce_bias', 'robust', 'padding', 'remove_eyes', 'surfaces', 't2_guided') - robust = traits.Bool(desc='robust brain centre estimation ' \ + robust = traits.Bool(desc='robust brain centre estimation ' '(iterates BET several times)', - argstr='-R', xor=_xor_inputs) - padding = traits.Bool(desc='improve BET if FOV is very small in Z ' \ + argstr='-R', xor=_xor_inputs) + padding = traits.Bool(desc='improve BET if FOV is very small in Z ' '(by temporarily padding end slices)', - argstr='-Z', xor=_xor_inputs) - remove_eyes = traits.Bool(desc='eye & optic nerve cleanup (can be ' \ + argstr='-Z', xor=_xor_inputs) + remove_eyes = traits.Bool(desc='eye & optic nerve cleanup (can be ' 'useful in SIENA)', - argstr='-S', xor=_xor_inputs) - surfaces = traits.Bool(desc='run bet2 and then betsurf to get additional ' \ - 'skull and scalp surfaces (includes ' \ + argstr='-S', xor=_xor_inputs) + surfaces = traits.Bool(desc='run bet2 and then betsurf to get additional ' + 'skull and scalp surfaces (includes ' 'registrations)', argstr='-A', xor=_xor_inputs) - t2_guided = File(desc='as with creating surfaces, when also feeding in ' \ + t2_guided = File(desc='as with creating surfaces, when also feeding in ' 'non-brain-extracted T2 (includes registrations)', argstr='-A2 %s', xor=_xor_inputs) functional = traits.Bool(argstr='-F', xor=_xor_inputs, @@ -156,31 +157,33 @@ def _list_outputs(self): if ((isdefined(self.inputs.mesh) and self.inputs.mesh) or (isdefined(self.inputs.surfaces) and self.inputs.surfaces)): outputs['meshfile'] = self._gen_fname(outputs['out_file'], - suffix='_mesh.vtk', - change_ext=False) + suffix='_mesh.vtk', + change_ext=False) if (isdefined(self.inputs.mask) and self.inputs.mask) or \ - (isdefined(self.inputs.reduce_bias) and \ - self.inputs.reduce_bias): + (isdefined(self.inputs.reduce_bias) and + self.inputs.reduce_bias): outputs['mask_file'] = self._gen_fname(outputs['out_file'], - suffix='_mask') + suffix='_mask') if isdefined(self.inputs.outline) and self.inputs.outline: outputs['outline_file'] = self._gen_fname(outputs['out_file'], - suffix='_overlay') + suffix='_overlay') if isdefined(self.inputs.surfaces) and self.inputs.surfaces: outputs['inskull_mask_file'] = self._gen_fname(outputs['out_file'], - suffix='_inskull_mask') + suffix='_inskull_mask') outputs['inskull_mesh_file'] = self._gen_fname(outputs['out_file'], - suffix='_inskull_mesh') - outputs['outskull_mask_file'] = self._gen_fname(outputs['out_file'], - suffix='_outskull_mask') - outputs['outskull_mesh_file'] = self._gen_fname(outputs['out_file'], - suffix='_outskull_mesh') + suffix='_inskull_mesh') + outputs[ + 'outskull_mask_file'] = self._gen_fname(outputs['out_file'], + suffix='_outskull_mask') + outputs[ + 'outskull_mesh_file'] = self._gen_fname(outputs['out_file'], + suffix='_outskull_mesh') outputs['outskin_mask_file'] = self._gen_fname(outputs['out_file'], - suffix='_outskin_mask') + suffix='_outskin_mask') outputs['outskin_mesh_file'] = self._gen_fname(outputs['out_file'], - suffix='_outskin_mesh') + suffix='_outskin_mesh') outputs['skull_mask_file'] = self._gen_fname(outputs['out_file'], - suffix='_skull_mask') + suffix='_skull_mask') if isdefined(self.inputs.no_output) and self.inputs.no_output: outputs['out_file'] = Undefined return outputs @@ -194,45 +197,46 @@ def _gen_filename(self, name): class FASTInputSpec(FSLCommandInputSpec): """ Defines inputs (trait classes) for FAST """ in_files = InputMultiPath(File(exists=True), copyfile=False, - desc='image, or multi-channel set of images, ' \ + desc='image, or multi-channel set of images, ' 'to be segmented', - argstr='%s', position=-1, mandatory=True) + argstr='%s', position=-1, mandatory=True) out_basename = File(desc='base name of output files', argstr='-o %s') # uses in_file name as basename if none given number_classes = traits.Range(low=1, high=10, argstr='-n %d', desc='number of tissue-type classes') output_biasfield = traits.Bool(desc='output estimated bias field', argstr='-b') - output_biascorrected = traits.Bool(desc='output restored image ' \ - '(bias-corrected image)', + output_biascorrected = traits.Bool(desc='output restored image ' + '(bias-corrected image)', argstr='-B') - img_type = traits.Enum((1, 2, 3), desc='int specifying type of image: ' \ - '(1 = T1, 2 = T2, 3 = PD)', + img_type = traits.Enum((1, 2, 3), desc='int specifying type of image: ' + '(1 = T1, 2 = T2, 3 = PD)', argstr='-t %d') bias_iters = traits.Range(low=1, high=10, argstr='-I %d', - desc='number of main-loop iterations during ' \ - 'bias-field removal') + desc='number of main-loop iterations during ' + 'bias-field removal') bias_lowpass = traits.Range(low=4, high=40, - desc='bias field smoothing extent (FWHM) ' \ - 'in mm', + desc='bias field smoothing extent (FWHM) ' + 'in mm', argstr='-l %d', units='mm') init_seg_smooth = traits.Range(low=0.0001, high=0.1, - desc='initial segmentation spatial ' \ - 'smoothness (during bias field ' \ - 'estimation)', + desc='initial segmentation spatial ' + 'smoothness (during bias field ' + 'estimation)', argstr='-f %.3f') - segments = traits.Bool(desc='outputs a separate binary image for each ' \ - 'tissue type', + segments = traits.Bool(desc='outputs a separate binary image for each ' + 'tissue type', argstr='-g') - init_transform = File(exists=True, desc=' initialise'\ - ' using priors', + init_transform = File(exists=True, desc=' initialise' + ' using priors', argstr='-a %s') - other_priors = InputMultiPath(File(exist=True), desc='alternative prior images', + other_priors = InputMultiPath( + File(exist=True), desc='alternative prior images', argstr='-A %s', minlen=3, maxlen=3) no_pve = traits.Bool(desc='turn off PVE (partial volume estimation)', - argstr='--nopve') + argstr='--nopve') no_bias = traits.Bool(desc='do not remove bias field', - argstr='-N') + argstr='-N') use_priors = traits.Bool(desc='use priors throughout', argstr='-P') # must also set -a!, # mutually inclusive?? @@ -241,15 +245,15 @@ class FASTInputSpec(FSLCommandInputSpec): # figure out how to # handle with traits. segment_iters = traits.Range(low=1, high=50, - desc='number of segmentation-initialisation'\ - ' iterations', + desc='number of segmentation-initialisation' + ' iterations', argstr='-W %d') mixel_smooth = traits.Range(low=0.0, high=1.0, desc='spatial smoothness for mixeltype', argstr='-R %.2f') iters_afterbias = traits.Range(low=1, hight=20, - desc='number of main-loop iterations ' \ - 'after bias-field removal', + desc='number of main-loop iterations ' + 'after bias-field removal', argstr='-O %d') hyper = traits.Range(low=0.0, high=1.0, desc='segmentation spatial smoothness', @@ -257,7 +261,7 @@ class FASTInputSpec(FSLCommandInputSpec): verbose = traits.Bool(desc='switch on diagnostic messages', argstr='-v') manual_seg = File(exists=True, desc='Filename containing intensities', - argstr='-s %s') + argstr='-s %s') probability_maps = traits.Bool(desc='outputs individual probability maps', argstr='-p') @@ -265,22 +269,22 @@ class FASTInputSpec(FSLCommandInputSpec): class FASTOutputSpec(TraitedSpec): """Specify possible outputs from FAST""" tissue_class_map = File(exists=True, - desc='path/name of binary segmented volume file' \ + desc='path/name of binary segmented volume file' ' one val for each class _seg') - tissue_class_files = OutputMultiPath(File(desc='path/name of binary segmented volumes ' \ - 'one file for each class _seg_x')) - restored_image = OutputMultiPath(File(desc='restored images (one for each input image) ' \ - 'named according to the input images _restore')) + tissue_class_files = OutputMultiPath(File(desc='path/name of binary segmented volumes ' + 'one file for each class _seg_x')) + restored_image = OutputMultiPath(File(desc='restored images (one for each input image) ' + 'named according to the input images _restore')) mixeltype = File(desc="path/name of mixeltype volume file _mixeltype") partial_volume_map = File(desc="path/name of partial volume file _pveseg") - partial_volume_files = OutputMultiPath(File(desc='path/name of partial volumes files ' \ - 'one for each class, _pve_x')) + partial_volume_files = OutputMultiPath(File(desc='path/name of partial volumes files ' + 'one for each class, _pve_x')) bias_field = OutputMultiPath(File(desc='Estimated bias field _bias')) - probability_maps = OutputMultiPath(File(desc='filenames, one for each class, for each ' \ - 'input, prob_x')) + probability_maps = OutputMultiPath(File(desc='filenames, one for each class, for each ' + 'input, prob_x')) class FAST(FSLCommand): @@ -331,9 +335,9 @@ def _list_outputs(self): suffix='_seg') if self.inputs.segments: outputs['tissue_class_files'] = [] - for i in range(nclasses): + for i in range(nclasses): outputs['tissue_class_files'].append( - self._gen_fname(basefile, suffix='_seg_%d' % i)) + self._gen_fname(basefile, suffix='_seg_%d' % i)) if isdefined(self.inputs.output_biascorrected): outputs['restored_image'] = [] if len(self.inputs.in_files) > 1: @@ -342,19 +346,21 @@ def _list_outputs(self): for val, f in enumerate(self.inputs.in_files): # image numbering is 1-based outputs['restored_image'].append( - self._gen_fname(basefile, suffix='_restore_%d' % (val + 1))) + self._gen_fname(basefile, suffix='_restore_%d' % (val + 1))) else: # single image segmentation has unnumbered output image outputs['restored_image'].append( - self._gen_fname(basefile, suffix='_restore')) + self._gen_fname(basefile, suffix='_restore')) outputs['mixeltype'] = self._gen_fname(basefile, suffix='_mixeltype') if not self.inputs.no_pve: - outputs['partial_volume_map'] = self._gen_fname(basefile, suffix='_pveseg') + outputs['partial_volume_map'] = self._gen_fname( + basefile, suffix='_pveseg') outputs['partial_volume_files'] = [] for i in range(nclasses): - outputs['partial_volume_files'].append(self._gen_fname(basefile, - suffix='_pve_%d' % i)) + outputs[ + 'partial_volume_files'].append(self._gen_fname(basefile, + suffix='_pve_%d' % i)) if self.inputs.output_biasfield: outputs['bias_field'] = [] if len(self.inputs.in_files) > 1: @@ -363,17 +369,17 @@ def _list_outputs(self): for val, f in enumerate(self.inputs.in_files): # image numbering is 1-based outputs['bias_field'].append( - self._gen_fname(basefile, suffix='_bias_%d' % (val + 1))) + self._gen_fname(basefile, suffix='_bias_%d' % (val + 1))) else: # single image segmentation has unnumbered output image outputs['bias_field'].append( - self._gen_fname(basefile, suffix='_bias')) + self._gen_fname(basefile, suffix='_bias')) if self.inputs.probability_maps: outputs['probability_maps'] = [] for i in range(nclasses): outputs['probability_maps'].append( - self._gen_fname(basefile, suffix='_prob_%d' % i)) + self._gen_fname(basefile, suffix='_prob_%d' % i)) return outputs @@ -387,16 +393,16 @@ class FLIRTInputSpec(FSLCommandInputSpec): reference = File(exists=True, argstr='-ref %s', mandatory=True, position=1, desc='reference file') out_file = File(argstr='-out %s', desc='registered output file', - genfile=True, position=2, hash_files=False) + genfile=True, position=2, hash_files=False) out_matrix_file = File(argstr='-omat %s', - desc='output affine matrix in 4x4 asciii format', - genfile=True, position=3, hash_files=False) + desc='output affine matrix in 4x4 asciii format', + genfile=True, position=3, hash_files=False) - out_log = File( desc='output log' ) + out_log = File(desc='output log') in_matrix_file = File(argstr='-init %s', desc='input 4x4 affine matrix') apply_xfm = traits.Bool(argstr='-applyxfm', requires=['in_matrix_file'], - desc='apply transformation supplied by in_matrix_file') + desc='apply transformation supplied by in_matrix_file') datatype = traits.Enum('char', 'short', 'int', 'float', 'double', argstr='-datatype %s', desc='force output data type') @@ -407,17 +413,17 @@ class FLIRTInputSpec(FSLCommandInputSpec): # XXX What is the difference between 'cost' and 'searchcost'? Are # these both necessary or do they map to the same variable. cost_func = traits.Enum('mutualinfo', 'corratio', 'normcorr', 'normmi', - 'leastsq', 'labeldiff', - argstr='-searchcost %s', - desc='cost function') + 'leastsq', 'labeldiff', + argstr='-searchcost %s', + desc='cost function') uses_qform = traits.Bool(argstr='-usesqform', - desc='initialize using sform or qform') + desc='initialize using sform or qform') display_init = traits.Bool(argstr='-displayinit', - desc='display initial matrix') + desc='display initial matrix') angle_rep = traits.Enum('quaternion', 'euler', argstr='-anglerep %s', desc='representation of rotation angles') - interp = traits.Enum('trilinear', 'nearestneighbour', 'sinc', + interp = traits.Enum('trilinear', 'nearestneighbour', 'sinc', 'spline', argstr='-interp %s', desc='final interpolation method used in reslicing') sinc_width = traits.Int(argstr='-sincwidth %d', units='voxels', @@ -435,7 +441,7 @@ class FLIRTInputSpec(FSLCommandInputSpec): min_sampling = traits.Float(argstr='-minsampling %f', units='mm', desc='set minimum voxel dimension for sampling') padding_size = traits.Int(argstr='-paddingsize %d', units='voxels', - desc='for applyxfm: interpolates outside image '\ + desc='for applyxfm: interpolates outside image ' 'by size') searchr_x = traits.List(traits.Int, minlen=2, maxlen=2, units='degrees', argstr='-searchrx %s', @@ -474,9 +480,9 @@ class FLIRTOutputSpec(TraitedSpec): out_file = File(exists=True, desc='path/name of registered file (if generated)') out_matrix_file = File(exists=True, - desc='path/name of calculated affine transform ' \ + desc='path/name of calculated affine transform ' '(if generated)') - out_log = File(desc='path/name of output log (if generated)' ) + out_log = File(desc='path/name of output log (if generated)') class FLIRT(FSLCommand): @@ -520,36 +526,37 @@ def _list_outputs(self): outputs['out_matrix_file'] = self._gen_fname(self.inputs.in_file, suffix='_flirt.mat', change_ext=False) - outputs['out_matrix_file'] = os.path.abspath(outputs['out_matrix_file']) + outputs['out_matrix_file'] = os.path.abspath( + outputs['out_matrix_file']) - if isdefined( self.inputs.save_log ) and self.inputs.save_log: + if isdefined(self.inputs.save_log) and self.inputs.save_log: outputs['out_log'] = self.inputs.out_log # Generate an out_log file if one is not provided - if not isdefined( outputs['out_log']): - outputs['out_log'] = self._gen_fname( self.inputs.in_file, + if not isdefined(outputs['out_log']): + outputs['out_log'] = self._gen_fname(self.inputs.in_file, suffix='_flirt.log', change_ext=False) - outputs['out_log'] = os.path.abspath( outputs['out_log'] ) + outputs['out_log'] = os.path.abspath(outputs['out_log']) return outputs def aggregate_outputs(self, runtime=None, needed_outputs=None): - outputs = super(FLIRT,self).aggregate_outputs(runtime=runtime, needed_outputs=needed_outputs) - if isdefined( self.inputs.save_log ) and self.inputs.save_log: + outputs = super(FLIRT, self).aggregate_outputs( + runtime=runtime, needed_outputs=needed_outputs) + if isdefined(self.inputs.save_log) and self.inputs.save_log: with open(outputs.out_log, "a") as text_file: - text_file.write(runtime.stdout + '\n' ) + text_file.write(runtime.stdout + '\n') return outputs def _parse_inputs(self, skip=None): skip = [] - if isdefined( self.inputs.save_log ) and self.inputs.save_log: - if not isdefined( self.inputs.verbose ) or self.inputs.verbose==0: - self.inputs.verbose=1 - - skip.append( 'save_log' ) + if isdefined(self.inputs.save_log) and self.inputs.save_log: + if not isdefined(self.inputs.verbose) or self.inputs.verbose == 0: + self.inputs.verbose = 1 - return super(FLIRT,self)._parse_inputs(skip=skip) + skip.append('save_log') + return super(FLIRT, self)._parse_inputs(skip=skip) def _gen_filename(self, name): if name in ('out_file', 'out_matrix_file'): @@ -557,8 +564,10 @@ def _gen_filename(self, name): else: return None + class ApplyXfmInputSpec(FLIRTInputSpec): - apply_xfm = traits.Bool(True, argstr='-applyxfm', requires=['in_matrix_file'], + apply_xfm = traits.Bool( + True, argstr='-applyxfm', requires=['in_matrix_file'], desc='apply transformation supplied by in_matrix_file', usedefault=True) @@ -591,27 +600,40 @@ class MCFLIRTInputSpec(FSLCommandInputSpec): desc="timeseries to motion-correct") out_file = File(argstr='-out %s', genfile=True, desc="file to write", hash_files=False) - cost = traits.Enum('mutualinfo', 'woods', 'corratio', 'normcorr', 'normmi', 'leastsquares', + cost = traits.Enum( + 'mutualinfo', 'woods', 'corratio', 'normcorr', 'normmi', 'leastsquares', argstr='-cost %s', desc="cost function to optimize") bins = traits.Int(argstr='-bins %d', desc="number of histogram bins") - dof = traits.Int(argstr='-dof %d', desc="degrees of freedom for the transformation") + dof = traits.Int( + argstr='-dof %d', desc="degrees of freedom for the transformation") ref_vol = traits.Int(argstr='-refvol %d', desc="volume to align frames to") - scaling = traits.Float(argstr='-scaling %.2f', desc="scaling factor to use") - smooth = traits.Float(argstr='-smooth %.2f', desc="smoothing factor for the cost function") - rotation = traits.Int(argstr='-rotation %d', desc="scaling factor for rotation tolerances") + scaling = traits.Float( + argstr='-scaling %.2f', desc="scaling factor to use") + smooth = traits.Float( + argstr='-smooth %.2f', desc="smoothing factor for the cost function") + rotation = traits.Int( + argstr='-rotation %d', desc="scaling factor for rotation tolerances") stages = traits.Int(argstr='-stages %d', desc="stages (if 4, perform final search with sinc interpolation") - init = File(exists=True, argstr='-init %s', desc="inital transformation matrix") + init = File(exists=True, argstr='-init %s', + desc="inital transformation matrix") interpolation = traits.Enum("spline", "nn", "sinc", argstr="-%s_final", desc="interpolation method for transformation") - use_gradient = traits.Bool(argstr='-gdt', desc="run search on gradient images") - use_contour = traits.Bool(argstr='-edge', desc="run search on contour images") + use_gradient = traits.Bool( + argstr='-gdt', desc="run search on gradient images") + use_contour = traits.Bool( + argstr='-edge', desc="run search on contour images") mean_vol = traits.Bool(argstr='-meanvol', desc="register to mean volume") - stats_imgs = traits.Bool(argstr='-stats', desc="produce variance and std. dev. images") - save_mats = traits.Bool(argstr='-mats', desc="save transformation matrices") - save_plots = traits.Bool(argstr='-plots', desc="save transformation parameters") - save_rms = traits.Bool(argstr='-rmsabs -rmsrel', desc="save rms displacement parameters") - ref_file = File(exists=True, argstr='-reffile %s', desc="target image for motion correction") + stats_imgs = traits.Bool( + argstr='-stats', desc="produce variance and std. dev. images") + save_mats = traits.Bool( + argstr='-mats', desc="save transformation matrices") + save_plots = traits.Bool( + argstr='-plots', desc="save transformation parameters") + save_rms = traits.Bool( + argstr='-rmsabs -rmsrel', desc="save rms displacement parameters") + ref_file = File(exists=True, argstr='-reffile %s', + desc="target image for motion correction") class MCFLIRTOutputSpec(TraitedSpec): @@ -620,7 +642,8 @@ class MCFLIRTOutputSpec(TraitedSpec): std_img = File(exists=True, desc="standard deviation image") mean_img = File(exists=True, desc="mean timeseries image") par_file = File(exists=True, desc="text-file with motion parameters") - mat_file = OutputMultiPath(File(exists=True), desc="transformation matrices") + mat_file = OutputMultiPath(File( + exists=True), desc="transformation matrices") rms_files = OutputMultiPath(File(exists=True), desc="absolute and relative displacement parameters") @@ -658,11 +681,11 @@ def _list_outputs(self): outputs['out_file'] = self._gen_outfilename() if isdefined(self.inputs.stats_imgs) and self.inputs.stats_imgs: - outputs['variance_img'] = self._gen_fname(outputs['out_file'] + \ + outputs['variance_img'] = self._gen_fname(outputs['out_file'] + '_variance.ext', cwd=cwd) - outputs['std_img'] = self._gen_fname(outputs['out_file'] + \ + outputs['std_img'] = self._gen_fname(outputs['out_file'] + '_sigma.ext', cwd=cwd) - outputs['mean_img'] = self._gen_fname(outputs['out_file'] + \ + outputs['mean_img'] = self._gen_fname(outputs['out_file'] + '_meanvol.ext', cwd=cwd) if isdefined(self.inputs.save_mats) and self.inputs.save_mats: _, filename = os.path.split(outputs['out_file']) @@ -706,7 +729,7 @@ class FNIRTInputSpec(FSLCommandInputSpec): inwarp_file = File(exists=True, argstr='--inwarp=%s', desc='name of file containing initial non-linear warps') in_intensitymap_file = File(exists=True, argstr='--intin=%s', - desc='name of file/files containing initial intensity maping'\ + desc='name of file/files containing initial intensity maping' 'usually generated by previos fnirt run') fieldcoeff_file = traits.Either(traits.Bool, File, argstr='--cout=%s', desc='name of output file with field coefficients or true') @@ -717,37 +740,41 @@ class FNIRTInputSpec(FSLCommandInputSpec): desc='name of output file with field or true', hash_files=False) jacobian_file = traits.Either(traits.Bool, File, argstr='--jout=%s', - desc='name of file for writing out the Jacobian'\ + desc='name of file for writing out the Jacobian' 'of the field (for diagnostic or VBM purposes)', hash_files=False) modulatedref_file = traits.Either(traits.Bool, File, argstr='--refout=%s', - desc='name of file for writing out intensity modulated'\ + desc='name of file for writing out intensity modulated' '--ref (for diagnostic purposes)', hash_files=False) out_intensitymap_file = traits.Either(traits.Bool, File, argstr='--intout=%s', - desc='name of files for writing information pertaining '\ + desc='name of files for writing information pertaining ' 'to intensity mapping', hash_files=False) log_file = File(argstr='--logout=%s', desc='Name of log-file', genfile=True, hash_files=False) - config_file = traits.Either(traits.Enum("T1_2_MNI152_2mm", "FA_2_FMRIB58_1mm"), File(exists=True), argstr='--config=%s', + config_file = traits.Either( + traits.Enum("T1_2_MNI152_2mm", "FA_2_FMRIB58_1mm"), File(exists=True), argstr='--config=%s', desc='Name of config file specifying command line arguments') refmask_file = File(exists=True, argstr='--refmask=%s', desc='name of file with mask in reference space') inmask_file = File(exists=True, argstr='--inmask=%s', desc='name of file with mask in input image space') - skip_refmask = traits.Bool(argstr='--applyrefmask=0', xor=['apply_refmask'], + skip_refmask = traits.Bool( + argstr='--applyrefmask=0', xor=['apply_refmask'], desc='Skip specified refmask if set, default false') skip_inmask = traits.Bool(argstr='--applyinmask=0', xor=['apply_inmask'], desc='skip specified inmask if set, default false') - apply_refmask = traits.List(traits.Enum(0, 1), argstr='--applyrefmask=%s', xor=['skip_refmask'], + apply_refmask = traits.List( + traits.Enum(0, 1), argstr='--applyrefmask=%s', xor=['skip_refmask'], desc='list of iterations to use reference mask on (1 to use, 0 to skip)', sep=",") - apply_inmask = traits.List(traits.Enum(0, 1), argstr='--applyinmask=%s', xor=['skip_inmask'], + apply_inmask = traits.List( + traits.Enum(0, 1), argstr='--applyinmask=%s', xor=['skip_inmask'], desc='list of iterations to use input mask on (1 to use, 0 to skip)', sep=",") skip_implicit_ref_masking = traits.Bool(argstr='--imprefm 0', - desc='skip implicit masking based on value'\ + desc='skip implicit masking based on value' 'in --ref image. Default = 0') skip_implicit_in_masking = traits.Bool(argstr='--impinm 0', - desc='skip implicit masking based on value'\ + desc='skip implicit masking based on value' 'in --in image. Default = 0') refmask_val = traits.Float(argstr='--imprefval=%f', desc='Value to mask out in --ref image. Default =0.0') @@ -762,7 +789,7 @@ class FNIRTInputSpec(FSLCommandInputSpec): sep=",") warp_resolution = traits.Tuple(traits.Int, traits.Int, traits.Int, argstr='--warpres=%d,%d,%d', - desc='(approximate) resolution (in mm) of warp basis '\ + desc='(approximate) resolution (in mm) of warp basis ' 'in x-, y- and z-direction, default 10, 10, 10') spline_order = traits.Int(argstr='--splineorder=%d', desc='Order of spline, 2->Qadratic spline, 3->Cubic spline. Default=3') @@ -774,7 +801,7 @@ class FNIRTInputSpec(FSLCommandInputSpec): argstr='--regmod=%s', desc='Model for regularisation of warp-field [membrane_energy bending_energy], default bending_energy') regularization_lambda = traits.List(traits.Float, argstr='--lambda=%s', - desc='Weight of regularisation, default depending on --ssqlambda and --regmod '\ + desc='Weight of regularisation, default depending on --ssqlambda and --regmod ' 'switches. See user documetation.', sep=",") skip_lambda_ssq = traits.Bool(argstr='--ssqlambda 0', desc='If true, lambda is not weighted by current ssq, default false') @@ -791,13 +818,15 @@ class FNIRTInputSpec(FSLCommandInputSpec): desc='Order of poynomial for mapping intensities, default 5') biasfield_resolution = traits.Tuple(traits.Int, traits.Int, traits.Int, argstr='--biasres=%d,%d,%d', - desc='Resolution (in mm) of bias-field modelling '\ + desc='Resolution (in mm) of bias-field modelling ' 'local intensities, default 50, 50, 50') bias_regularization_lambda = traits.Float(argstr='--biaslambda=%f', desc='Weight of regularisation for bias-field, default 10000') - skip_intensity_mapping = traits.Bool(argstr='--estint=0', xor=['apply_intensity_mapping'], + skip_intensity_mapping = traits.Bool( + argstr='--estint=0', xor=['apply_intensity_mapping'], desc='Skip estimate intensity-mapping default false') - apply_intensity_mapping = traits.List(traits.Enum(0, 1), argstr='--estint=%s', xor=['skip_intensity_mapping'], + apply_intensity_mapping = traits.List( + traits.Enum(0, 1), argstr='--estint=%s', xor=['skip_intensity_mapping'], desc='List of subsampling levels to apply intensity mapping for (0 to skip, 1 to apply)', sep=",") hessian_precision = traits.Enum('double', 'float', argstr='--numprec=%s', desc='Precision for representing Hessian, double or float. Default double') @@ -809,7 +838,7 @@ class FNIRTOutputSpec(TraitedSpec): field_file = File(desc='file with warp field') jacobian_file = File(desc='file containing Jacobian of the field') modulatedref_file = File(desc='file containing intensity modulated --ref') - out_intensitymap_file = File(\ + out_intensitymap_file = File( desc='file containing info pertaining to intensity mapping') log_file = File(desc='Name of log-file') @@ -937,7 +966,8 @@ class ApplyWarpInputSpec(FSLCommandInputSpec): desc='filename for post-transform (affine matrix)') mask_file = File(exists=True, argstr='--mask=%s', desc='filename for mask image (in reference space)') - interp = traits.Enum('nn', 'trilinear', 'sinc', 'spline', argstr='--interp=%s', + interp = traits.Enum( + 'nn', 'trilinear', 'sinc', 'spline', argstr='--interp=%s', desc='interpolation method') @@ -1008,7 +1038,8 @@ class SliceTimerInputSpec(FSLCommandInputSpec): class SliceTimerOutputSpec(TraitedSpec): - slice_time_corrected_file = File(exists=True, desc='slice time corrected file') + slice_time_corrected_file = File( + exists=True, desc='slice time corrected file') class SliceTimer(FSLCommand): @@ -1060,7 +1091,8 @@ class SUSANInputSpec(FSLCommandInputSpec): desc='within-plane (2) or fully 3D (3)') use_median = traits.Enum(1, 0, argstr='%d', position=5, usedefault=True, desc='whether to use a local median filter in the cases where single-point noise is detected') - usans = traits.List(traits.Tuple(File(exists=True), traits.Float), maxlen=2, + usans = traits.List( + traits.Tuple(File(exists=True), traits.Float), maxlen=2, argstr='', position=6, default=[], usedefault=True, desc='determines whether the smoothing area (USAN) is to be ' 'found from secondary images (0, 1 or 2). A negative ' @@ -1128,7 +1160,7 @@ class FUGUEInputSpec(FSLCommandInputSpec): unwarped_file = File(argstr='--unwarp=%s', genfile=True, desc='apply unwarping and save as filename', hash_files=False) - save_warped = traits.Bool( desc='apply forward warp and save' ) + save_warped = traits.Bool(desc='apply forward warp and save') warped_file = File(argstr='--warp=%s', genfile=True, desc='apply forward warp and save as filename', hash_files=False) @@ -1146,9 +1178,9 @@ class FUGUEInputSpec(FSLCommandInputSpec): fmap_in_file = File(exists=True, argstr='--loadfmap=%s', desc='filename for loading fieldmap (rad/s)') - save_shift = traits.Bool( desc='output pixel shift volume' ) + save_shift = traits.Bool(desc='output pixel shift volume') - shift_out_file = traits.File( argstr='--saveshift=%s', genfile=True, + shift_out_file = traits.File(argstr='--saveshift=%s', genfile=True, desc='filename for saving pixel shift volume', hash_files=False) shift_in_file = File(exists=True, argstr='--loadshift=%s', @@ -1194,13 +1226,14 @@ class FUGUEInputSpec(FSLCommandInputSpec): argstr='--unmaskshift=%s', requires=['shift_out_file'], desc='saves the unmasked shiftmap when using --saveshift', hash_files=False) - nokspace = traits.Bool(argstr='--nokspace', desc='do not use k-space forward warping') + nokspace = traits.Bool( + argstr='--nokspace', desc='do not use k-space forward warping') class FUGUEOutputSpec(TraitedSpec): unwarped_file = File(exists=True, desc='unwarped file') - shift_out_file = File( desc='voxel shift map file' ) - warped_file = File( desc='warped file' ) + shift_out_file = File(desc='voxel shift map file') + warped_file = File(desc='warped file') class FUGUE(FSLCommand): @@ -1219,18 +1252,19 @@ class FUGUE(FSLCommand): def __init__(self, **kwargs): super(FUGUE, self).__init__(**kwargs) - warn('This interface has not been fully tested. Please report any failures.') + warn( + 'This interface has not been fully tested. Please report any failures.') def _parse_inputs(self, skip=None): skip = [] - if not isdefined( self.inputs.save_shift ) or not self.inputs.save_shift: + if not isdefined(self.inputs.save_shift) or not self.inputs.save_shift: skip.append('shift_out_file') - if not isdefined( self.inputs.save_warped ) or not self.inputs.save_warped: + if not isdefined(self.inputs.save_warped) or not self.inputs.save_warped: skip.append('warped_file') else: skip.append('unwarped_file') - return super(FUGUE,self)._parse_inputs(skip=skip) + return super(FUGUE, self)._parse_inputs(skip=skip) def _list_outputs(self): outputs = self._outputs().get() @@ -1240,19 +1274,21 @@ def _list_outputs(self): suffix='_unwarped') outputs['unwarped_file'] = os.path.abspath(out_file) - if isdefined( self.inputs.save_shift ) and self.inputs.save_shift: + if isdefined(self.inputs.save_shift) and self.inputs.save_shift: shift_out = self.inputs.shift_out_file if not isdefined(shift_out): - shift_out = self._gen_fname( self.inputs.in_file, suffix='_shift' ) + shift_out = self._gen_fname( + self.inputs.in_file, suffix='_shift') - outputs['shift_out_file'] = os.path.abspath( shift_out ) + outputs['shift_out_file'] = os.path.abspath(shift_out) - if isdefined( self.inputs.save_warped ) and self.inputs.save_warped: + if isdefined(self.inputs.save_warped) and self.inputs.save_warped: warped_out = self.inputs.warped_file if not isdefined(warped_out): - warped_out = self._gen_fname( self.inputs.in_file, suffix='_fwdwarp' ) + warped_out = self._gen_fname( + self.inputs.in_file, suffix='_fwdwarp') - outputs['warped_file'] = os.path.abspath( warped_out ) + outputs['warped_file'] = os.path.abspath(warped_out) del outputs['unwarped_file'] return outputs @@ -1272,7 +1308,8 @@ def _gen_filename(self, name): class PRELUDEInputSpec(FSLCommandInputSpec): complex_phase_file = File(exists=True, argstr='--complex=%s', - mandatory=True, xor=['magnitude_file', 'phase_file'], + mandatory=True, xor=[ + 'magnitude_file', 'phase_file'], desc='complex phase input volume') magnitude_file = File(exists=True, argstr='--abs=%s', mandatory=True, diff --git a/nipype/workflows/dmri/fsl/dti.py b/nipype/workflows/dmri/fsl/dti.py index e061e2207c..3b95bb5356 100644 --- a/nipype/workflows/dmri/fsl/dti.py +++ b/nipype/workflows/dmri/fsl/dti.py @@ -5,18 +5,27 @@ import nipype.interfaces.fsl as fsl import os + def transpose(samples_over_fibres): import numpy as np a = np.array(samples_over_fibres) - if len(a.shape)==1: - a = a.reshape(-1,1) + if len(a.shape) == 1: + a = a.reshape(-1, 1) return a.T.tolist() + def create_dmri_preprocessing(name="dMRI_preprocessing", fieldmap_registration=False): """Creates a workflow that chains the necessary pipelines to correct for motion, eddy currents, and susceptibility artifacts in EPI dMRI sequences. + .. warning:: + + IMPORTANT NOTICE: this workflow rotates the b-vectors, so please be adviced + that not all the dicom converters ensure the consistency between the resulting + nifti orientation and the b matrix table (e.g. dcm2nii checks it). + + Example ------- @@ -35,16 +44,16 @@ def create_dmri_preprocessing(name="dMRI_preprocessing", fieldmap_registration=F Inputs:: - inputnode.in_file - inputnode.in_bvec - inputnode.ref_num - inputnode.fieldmap_mag - inputnode.fieldmap_pha - inputnode.te_diff - inputnode.epi_echospacing - inputnode.epi_rev_encoding - inputnode.pi_accel_factor - inputnode.vsm_sigma + inputnode.in_file - The diffusion data + inputnode.in_bvec - The b-matrix file, in FSL format and consistent with the in_file orientation + inputnode.ref_num - The reference volume (a b=0 volume in dMRI) + inputnode.fieldmap_mag - The magnitude of the fieldmap + inputnode.fieldmap_pha - The phase difference of the fieldmap + inputnode.te_diff - TE increment used (in msec.) on the fieldmap acquisition (generally 2.46ms for 3T scanners) + inputnode.epi_echospacing - The EPI EchoSpacing parameter (in msec.) + inputnode.epi_rev_encoding - True if reverse encoding was used (generally False) + inputnode.pi_accel_factor - Parallel imaging factor (aka GRAPPA acceleration factor) + inputnode.vsm_sigma - Sigma (in mm.) of the gaussian kernel used for in-slice smoothing of the deformation field (voxel shift map, vsm) Outputs:: @@ -55,47 +64,35 @@ def create_dmri_preprocessing(name="dMRI_preprocessing", fieldmap_registration=F Optional arguments:: - fieldmap_registration - True if registration to fieldmap should be done (default False) + fieldmap_registration - True if registration to fieldmap should be performed (default False) """ pipeline = pe.Workflow(name=name) - - inputnode = pe.Node(interface = util.IdentityInterface( - fields=['in_file', 'in_bvec','ref_num','fieldmap_mag', - 'fieldmap_pha','te_diff','epi_echospacing', - 'epi_rev_encoding','pi_accel_factor','vsm_sigma']), - name='inputnode') - outputnode= pe.Node(interface = util.IdentityInterface( - fields=['dmri_corrected','bvec_rotated']), - name='outputnode') + inputnode = pe.Node(interface=util.IdentityInterface( + fields=['in_file', 'in_bvec', 'ref_num', 'fieldmap_mag', + 'fieldmap_pha', 'te_diff', 'epi_echospacing', + 'epi_rev_encoding', 'pi_accel_factor', 'vsm_sigma']), + name='inputnode') + + outputnode = pe.Node(interface=util.IdentityInterface( + fields=['dmri_corrected', 'bvec_rotated']), + name='outputnode') motion = create_motion_correct_pipeline() eddy = create_eddy_correct_pipeline() - susceptibility = create_susceptibility_correct_pipeline(fieldmap_registration=fieldmap_registration) + susceptibility = create_susceptibility_correct_pipeline( + fieldmap_registration=fieldmap_registration) pipeline.connect([ - (inputnode, motion, [('in_file','inputnode.in_file'),('in_bvec','inputnode.in_bvec'),('ref_num','inputnode.ref_num')]) - ,(inputnode, eddy, [('ref_num','inputnode.ref_num')]) - ,(motion, eddy, [('outputnode.motion_corrected','inputnode.in_file')]) - ,(eddy, susceptibility, [('outputnode.eddy_corrected','inputnode.in_file')]) - ,(inputnode, susceptibility, [('ref_num','inputnode.ref_num') - ,('fieldmap_mag','inputnode.fieldmap_mag') - ,('fieldmap_pha','inputnode.fieldmap_pha') - ,('te_diff','inputnode.te_diff') - ,('epi_echospacing','inputnode.epi_echospacing') - ,('epi_rev_encoding','inputnode.epi_rev_encoding') - ,('pi_accel_factor','inputnode.pi_accel_factor') - ,('vsm_sigma','inputnode.vsm_sigma') ]) - ,(motion, outputnode, [('outputnode.out_bvec', 'bvec_rotated')]) - ,(susceptibility, outputnode, [('outputnode.epi_corrected','dmri_corrected')]) - ]) + (inputnode, motion, [('in_file', 'inputnode.in_file'), ( + 'in_bvec', 'inputnode.in_bvec'), ('ref_num', 'inputnode.ref_num')]), (inputnode, eddy, [('ref_num', 'inputnode.ref_num')]), (motion, eddy, [('outputnode.motion_corrected', 'inputnode.in_file')]), (eddy, susceptibility, [('outputnode.eddy_corrected', 'inputnode.in_file')]), (inputnode, susceptibility, [('ref_num', 'inputnode.ref_num'), ('fieldmap_mag', 'inputnode.fieldmap_mag'), ('fieldmap_pha', 'inputnode.fieldmap_pha'), ('te_diff', 'inputnode.te_diff'), ('epi_echospacing', 'inputnode.epi_echospacing'), ('epi_rev_encoding', 'inputnode.epi_rev_encoding'), ('pi_accel_factor', 'inputnode.pi_accel_factor'), ('vsm_sigma', 'inputnode.vsm_sigma')]), (motion, outputnode, [('outputnode.out_bvec', 'bvec_rotated')]), (susceptibility, outputnode, [('outputnode.epi_corrected', 'dmri_corrected')]) + ]) return pipeline - def create_bedpostx_pipeline(name="bedpostx"): """Creates a pipeline that does the same as bedpostx script from FSL - @@ -135,14 +132,14 @@ def create_bedpostx_pipeline(name="bedpostx"): """ + inputnode = pe.Node( + interface=util.IdentityInterface(fields=["dwi", "mask"]), + name="inputnode") - inputnode = pe.Node(interface = util.IdentityInterface(fields=["dwi", "mask"]), - name="inputnode") - - mask_dwi = pe.Node(interface = fsl.ImageMaths(op_string = "-mas"), + mask_dwi = pe.Node(interface=fsl.ImageMaths(op_string="-mas"), name="mask_dwi") - slice_dwi = pe.Node(interface = fsl.Split(dimension="z"), name="slice_dwi") - slice_mask = pe.Node(interface = fsl.Split(dimension="z"), + slice_dwi = pe.Node(interface=fsl.Split(dimension="z"), name="slice_dwi") + slice_mask = pe.Node(interface=fsl.Split(dimension="z"), name="slice_mask") preproc = pe.Workflow(name="preproc") @@ -156,7 +153,6 @@ def create_bedpostx_pipeline(name="bedpostx"): xfibres = pe.MapNode(interface=fsl.XFibres(), name="xfibres", iterfield=['dwi', 'mask']) - # Normal set of parameters xfibres.inputs.n_fibres = 2 xfibres.inputs.fudge = 1 @@ -167,12 +163,12 @@ def create_bedpostx_pipeline(name="bedpostx"): xfibres.inputs.non_linear = True xfibres.inputs.update_proposal_every = 24 - inputnode = pe.Node(interface = util.IdentityInterface(fields=["thsamples", - "phsamples", - "fsamples", - "dyads", - "mean_dsamples", - "mask"]), + inputnode = pe.Node(interface=util.IdentityInterface(fields=["thsamples", + "phsamples", + "fsamples", + "dyads", + "mean_dsamples", + "mask"]), name="inputnode") merge_thsamples = pe.MapNode(fsl.Merge(dimension="z"), @@ -182,7 +178,6 @@ def create_bedpostx_pipeline(name="bedpostx"): merge_fsamples = pe.MapNode(fsl.Merge(dimension="z"), name="merge_fsamples", iterfield=['in_files']) - merge_mean_dsamples = pe.Node(fsl.Merge(dimension="z"), name="merge_mean_dsamples") @@ -197,24 +192,33 @@ def create_bedpostx_pipeline(name="bedpostx"): postproc = pe.Workflow(name="postproc") - postproc.connect([(inputnode, merge_thsamples, [(('thsamples',transpose), 'in_files')]), - (inputnode, merge_phsamples, [(('phsamples',transpose), 'in_files')]), - (inputnode, merge_fsamples, [(('fsamples',transpose), 'in_files')]), - (inputnode, merge_mean_dsamples, [('mean_dsamples', 'in_files')]), - - (merge_thsamples, mean_thsamples, [('merged_file', 'in_file')]), - (merge_phsamples, mean_phsamples, [('merged_file', 'in_file')]), - (merge_fsamples, mean_fsamples, [('merged_file', 'in_file')]), - (merge_thsamples, make_dyads, [('merged_file', 'theta_vol')]), - (merge_phsamples, make_dyads, [('merged_file', 'phi_vol')]), - (inputnode, make_dyads, [('mask', 'mask')]), - ]) - - inputnode = pe.Node(interface = util.IdentityInterface(fields=["dwi", - "mask", - "bvecs", - "bvals"]), - name="inputnode") + postproc.connect( + [(inputnode, merge_thsamples, [(('thsamples', transpose), 'in_files')]), + (inputnode, merge_phsamples, [(( + 'phsamples', transpose), 'in_files')]), + (inputnode, merge_fsamples, [(( + 'fsamples', transpose), 'in_files')]), + (inputnode, merge_mean_dsamples, [ + ('mean_dsamples', 'in_files')]), + + (merge_thsamples, mean_thsamples, [ + ('merged_file', 'in_file')]), + (merge_phsamples, mean_phsamples, [ + ('merged_file', 'in_file')]), + (merge_fsamples, mean_fsamples, [ + ('merged_file', 'in_file')]), + (merge_thsamples, make_dyads, [ + ('merged_file', 'theta_vol')]), + (merge_phsamples, make_dyads, [ + ('merged_file', 'phi_vol')]), + (inputnode, make_dyads, [('mask', 'mask')]), + ]) + + inputnode = pe.Node(interface=util.IdentityInterface(fields=["dwi", + "mask", + "bvecs", + "bvals"]), + name="inputnode") bedpostx = pe.Workflow(name=name) bedpostx.connect([(inputnode, preproc, [('mask', 'inputnode.mask')]), @@ -226,40 +230,57 @@ def create_bedpostx_pipeline(name="bedpostx"): (inputnode, xfibres, [('bvecs', 'bvecs')]), (inputnode, postproc, [('mask', 'inputnode.mask')]), - (xfibres, postproc, [('thsamples','inputnode.thsamples'), - ('phsamples', 'inputnode.phsamples'), - ('fsamples', 'inputnode.fsamples'), - ('dyads', 'inputnode.dyads'), - ('mean_dsamples', 'inputnode.mean_dsamples')]), - ]) - - outputnode = pe.Node(interface = util.IdentityInterface(fields=["thsamples", - "phsamples", - "fsamples", - "mean_thsamples", - "mean_phsamples", - "mean_fsamples", - "dyads", - "dyads_dispersion"]), - name="outputnode") - bedpostx.connect([(postproc, outputnode, [("merge_thsamples.merged_file", "thsamples"), - ("merge_phsamples.merged_file", "phsamples"), - ("merge_fsamples.merged_file", "fsamples"), - ("mean_thsamples.out_file", "mean_thsamples"), - ("mean_phsamples.out_file", "mean_phsamples"), - ("mean_fsamples.out_file", "mean_fsamples"), + (xfibres, postproc, [ + ('thsamples', 'inputnode.thsamples'), + ('phsamples', + 'inputnode.phsamples'), + ('fsamples', 'inputnode.fsamples'), + ('dyads', 'inputnode.dyads'), + ('mean_dsamples', 'inputnode.mean_dsamples')]), + ]) + + outputnode = pe.Node( + interface=util.IdentityInterface(fields=["thsamples", + "phsamples", + "fsamples", + "mean_thsamples", + "mean_phsamples", + "mean_fsamples", + "dyads", + "dyads_dispersion"]), + name="outputnode") + bedpostx.connect( + [(postproc, outputnode, [("merge_thsamples.merged_file", "thsamples"), + ("merge_phsamples.merged_file", + "phsamples"), + ("merge_fsamples.merged_file", + "fsamples"), + ("mean_thsamples.out_file", + "mean_thsamples"), + ("mean_phsamples.out_file", + "mean_phsamples"), + ("mean_fsamples.out_file", + "mean_fsamples"), ("make_dyads.dyads", "dyads"), ("make_dyads.dispersion", "dyads_dispersion")]) ]) return bedpostx + def create_motion_correct_pipeline(name="motion_correct"): """Creates a pipeline that corrects for motion artifact in dMRI sequences. - It takes a series of diffusion weighted images and rigidly corregisters + It takes a series of diffusion weighted images and rigidly corregisters them to one reference image. Finally, the b-matrix is rotated accordingly (Leemans et al. 2009 - http://www.ncbi.nlm.nih.gov/pubmed/19319973), making use of the rotation matrix obtained by FLIRT. + .. warning:: + + IMPORTANT NOTICE: this workflow rotates the b-vectors, so please be adviced + that not all the dicom converters ensure the consistency between the resulting + nifti orientation and the b matrix table (e.g. dcm2nii checks it). + + Example ------- @@ -282,30 +303,28 @@ def create_motion_correct_pipeline(name="motion_correct"): """ - inputnode = pe.Node(interface = util.IdentityInterface(fields=["in_file", "ref_num","in_bvec"]), + inputnode = pe.Node( + interface=util.IdentityInterface( + fields=["in_file", "ref_num", "in_bvec"]), name="inputnode") pipeline = pe.Workflow(name=name) split = pe.Node(fsl.Split(dimension='t'), name="split") pick_ref = pe.Node(util.Select(), name="pick_ref") - coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1, dof=6), name = "coregistration", iterfield=["in_file"]) - rotate_bvecs = pe.Node( util.Function( input_names=["in_bvec", "in_matrix"], output_names=["out_file"], function=_rotate_bvecs ), name="rotate_b_matrix" ) + coregistration = pe.MapNode(fsl.FLIRT(no_search=True, interp='spline', + padding_size=1, dof=6), name="coregistration", iterfield=["in_file"]) + rotate_bvecs = pe.Node(util.Function(input_names=["in_bvec", "in_matrix"], output_names=[ + "out_file"], function=_rotate_bvecs), name="rotate_b_matrix") merge = pe.Node(fsl.Merge(dimension="t"), name="merge") - outputnode = pe.Node(interface = util.IdentityInterface(fields=["motion_corrected","out_bvec"]), + outputnode = pe.Node( + interface=util.IdentityInterface( + fields=["motion_corrected", "out_bvec"]), name="outputnode") pipeline.connect([ - (inputnode, split, [("in_file", "in_file")]) - ,(split, pick_ref, [("out_files", "inlist")]) - ,(inputnode, pick_ref, [("ref_num", "index")]) - ,(split, coregistration, [("out_files", "in_file")]) - ,(inputnode, rotate_bvecs, [("in_bvec","in_bvec")]) - ,(coregistration, rotate_bvecs, [("out_matrix_file","in_matrix")]) - ,(pick_ref, coregistration, [("out", "reference")]) - ,(coregistration, merge, [("out_file", "in_files")]) - ,(merge, outputnode, [("merged_file", "motion_corrected")]) - ,(rotate_bvecs, outputnode, [("out_file","out_bvec")]) + (inputnode, split, [("in_file", "in_file")]), ( + split, pick_ref, [("out_files", "inlist")]), (inputnode, pick_ref, [("ref_num", "index")]), (split, coregistration, [("out_files", "in_file")]), (inputnode, rotate_bvecs, [("in_bvec", "in_bvec")]), (coregistration, rotate_bvecs, [("out_matrix_file", "in_matrix")]), (pick_ref, coregistration, [("out", "reference")]), (coregistration, merge, [("out_file", "in_files")]), (merge, outputnode, [("merged_file", "motion_corrected")]), (rotate_bvecs, outputnode, [("out_file", "out_bvec")]) ]) return pipeline @@ -335,30 +354,29 @@ def create_eddy_correct_pipeline(name="eddy_correct"): outputnode.eddy_corrected """ - inputnode = pe.Node(interface = util.IdentityInterface(fields=["in_file", "ref_num"]), + inputnode = pe.Node( + interface=util.IdentityInterface(fields=["in_file", "ref_num"]), name="inputnode") pipeline = pe.Workflow(name=name) split = pe.Node(fsl.Split(dimension='t'), name="split") pick_ref = pe.Node(util.Select(), name="pick_ref") - coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1, dof=12 ), name = "coregistration", iterfield=["in_file"]) + coregistration = pe.MapNode(fsl.FLIRT(no_search=True, padding_size=1, + dof=12, interp='spline'), name="coregistration", iterfield=["in_file"]) merge = pe.Node(fsl.Merge(dimension="t"), name="merge") - outputnode = pe.Node(interface = util.IdentityInterface(fields=["eddy_corrected"]), + outputnode = pe.Node( + interface=util.IdentityInterface(fields=["eddy_corrected"]), name="outputnode") pipeline.connect([ - (inputnode, split, [("in_file", "in_file")]) - ,(split, pick_ref, [("out_files", "inlist")]) - ,(inputnode, pick_ref, [("ref_num", "index")]) - ,(split, coregistration, [("out_files", "in_file")]) - ,(pick_ref, coregistration, [("out", "reference")]) - ,(coregistration, merge, [("out_file", "in_files")]) - ,(merge, outputnode, [("merged_file", "eddy_corrected")]) + (inputnode, split, [("in_file", "in_file")]), ( + split, pick_ref, [("out_files", "inlist")]), (inputnode, pick_ref, [("ref_num", "index")]), (split, coregistration, [("out_files", "in_file")]), (pick_ref, coregistration, [("out", "reference")]), (coregistration, merge, [("out_file", "in_files")]), (merge, outputnode, [("merged_file", "eddy_corrected")]) ]) return pipeline -def create_susceptibility_correct_pipeline(name="susceptibility_correct", fieldmap_registration=False ): + +def create_susceptibility_correct_pipeline(name="susceptibility_correct", fieldmap_registration=False): """ Replaces the epidewarp.fsl script (http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl) for susceptibility distortion correction of dMRI & fMRI acquired with EPI sequences and the fieldmap information (Jezzard et al., 1995) using FSL's FUGUE. The registration to the (warped) fieldmap @@ -400,16 +418,16 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", fieldm Optional arguments:: fieldmap_registration - True if registration to fieldmap should be done (default False) - + """ - inputnode = pe.Node(interface = util.IdentityInterface(fields=["in_file", - "fieldmap_mag", - "fieldmap_pha", - "te_diff", + inputnode = pe.Node(interface=util.IdentityInterface(fields=["in_file", + "fieldmap_mag", + "fieldmap_pha", + "te_diff", "epi_echospacing", - "epi_ph_encoding_dir", - "epi_rev_encoding", + "epi_ph_encoding_dir", + "epi_rev_encoding", "pi_accel_factor", "vsm_sigma", "ref_num", @@ -419,217 +437,221 @@ def create_susceptibility_correct_pipeline(name="susceptibility_correct", fieldm pipeline = pe.Workflow(name=name) # Keep first frame from magnitude - select_mag = pe.Node( interface=fsl.utils.ExtractROI(t_size=1,t_min=0), name="select_magnitude" ) - + select_mag = pe.Node(interface=fsl.utils.ExtractROI( + t_size=1, t_min=0), name="select_magnitude") + # mask_brain - mask_mag = pe.Node( interface=fsl.BET(mask=True),name='mask_magnitude' ) - mask_mag_dil = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_dilate_mask), name='mask_dilate' ) + mask_mag = pe.Node(interface=fsl.BET(mask=True), name='mask_magnitude') + mask_mag_dil = pe.Node(interface=util.Function(input_names=[ + "in_file"], output_names=["out_file"], function=_dilate_mask), name='mask_dilate') # Compute dwell time - dwell_time = pe.Node( interface=util.Function( input_names=["dwell_time","pi_factor","is_reverse_encoding"], output_names=["dwell_time"], function=_compute_dwelltime), name='dwell_time') + dwell_time = pe.Node(interface=util.Function(input_names=["dwell_time", "pi_factor", "is_reverse_encoding"], output_names=[ + "dwell_time"], function=_compute_dwelltime), name='dwell_time') # Normalize phase diff to be [-pi, pi) - norm_pha = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_prepare_phasediff ), name='normalize_phasediff') + norm_pha = pe.Node(interface=util.Function(input_names=["in_file"], output_names=[ + "out_file"], function=_prepare_phasediff), name='normalize_phasediff') # Execute FSL PRELUDE: prelude -p %s -a %s -o %s -f -v -m %s - prelude = pe.Node( interface=fsl.PRELUDE(process3d=True), name='phase_unwrap' ) - fill_phase = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_fill_phase ), name='fill_phasediff' ) + prelude = pe.Node(interface=fsl.PRELUDE( + process3d=True), name='phase_unwrap') + fill_phase = pe.Node(interface=util.Function(input_names=["in_file"], output_names=[ + "out_file"], function=_fill_phase), name='fill_phasediff') # to assure that vsm is same dimension as mag. The input only affects the output dimension. # The content of the input has no effect on the vsm. The de-warped mag volume is # meaningless and will be thrown away - # fugue -i %s -u %s -p %s --dwell=%s --asym=%s --mask=%s --saveshift=%s % ( mag_name, magdw_name, ph_name, esp, tediff, mask_name, vsmmag_name) - vsm = pe.Node( interface=fsl.FUGUE(save_shift=True), name="generate_vsm" ) - vsm_mean = pe.Node( interface=util.Function( input_names=["in_file","mask_file","in_unwarped"], output_names=["out_file"], function=_vsm_remove_mean), name="vsm_mean_shift" ) - + # fugue -i %s -u %s -p %s --dwell=%s --asym=%s --mask=%s --saveshift=%s % + # ( mag_name, magdw_name, ph_name, esp, tediff, mask_name, vsmmag_name) + vsm = pe.Node(interface=fsl.FUGUE(save_shift=True), name="generate_vsm") + vsm_mean = pe.Node(interface=util.Function(input_names=["in_file", "mask_file", "in_unwarped"], output_names=[ + "out_file"], function=_vsm_remove_mean), name="vsm_mean_shift") + # fugue_epi - dwi_split = pe.Node( interface=util.Function( input_names=['in_file'], output_names=['out_files'], function=_split_dwi ), name='dwi_split' ) + dwi_split = pe.Node(interface=util.Function(input_names=[ + 'in_file'], output_names=['out_files'], function=_split_dwi), name='dwi_split') # 'fugue -i %s -u %s --loadshift=%s --mask=%s' % ( vol_name, out_vol_name, vsm_name, mask_name ) - dwi_applyxfm = pe.MapNode( interface=fsl.FUGUE(icorr=True,save_shift=False), iterfield=['in_file'], name='dwi_fugue' ) + dwi_applyxfm = pe.MapNode(interface=fsl.FUGUE( + icorr=True, save_shift=False), iterfield=['in_file'], name='dwi_fugue') # Merge back all volumes - dwi_merge = pe.Node( interface=fsl.utils.Merge(dimension='t' ), name='dwi_merge') + dwi_merge = pe.Node(interface=fsl.utils.Merge( + dimension='t'), name='dwi_merge') - outputnode = pe.Node(interface = util.IdentityInterface(fields=["epi_corrected"]), + outputnode = pe.Node( + interface=util.IdentityInterface(fields=["epi_corrected"]), name="outputnode") - - pipeline.connect([ - (inputnode, dwell_time, [('epi_echospacing','dwell_time'), ('pi_accel_factor','pi_factor'),('epi_rev_encoding','is_reverse_encoding') ]) - ,(inputnode, select_mag, [('fieldmap_mag','in_file')] ) - ,(inputnode, norm_pha, [('fieldmap_pha','in_file')] ) - ,(select_mag, mask_mag, [('roi_file', 'in_file')] ) - ,(mask_mag, mask_mag_dil, [('mask_file','in_file')] ) - ,(select_mag, prelude, [('roi_file', 'magnitude_file') ]) - ,(norm_pha, prelude, [('out_file', 'phase_file' ) ]) - ,(mask_mag_dil, prelude, [('out_file', 'mask_file' ) ]) - ,(prelude, fill_phase, [('unwrapped_phase_file','in_file')] ) - ,(inputnode, vsm, [('fieldmap_mag', 'in_file')]) - ,(fill_phase, vsm, [('out_file','phasemap_file')]) - ,(inputnode, vsm, [(('te_diff',_ms2sec),'asym_se_time'),('vsm_sigma','smooth2d')]) - ,(dwell_time, vsm, [(('dwell_time',_ms2sec),'dwell_time') ]) - ,(mask_mag_dil, vsm, [('out_file','mask_file')] ) - ,(mask_mag_dil, vsm_mean, [('out_file','mask_file')] ) - ,(vsm, vsm_mean, [('unwarped_file','in_unwarped'),('shift_out_file','in_file')]) - ,(inputnode, dwi_split, [('in_file','in_file')]) - ,(dwi_split, dwi_applyxfm, [('out_files','in_file')]) - ,(dwi_applyxfm, dwi_merge, [('unwarped_file','in_files')]) - ,(dwi_merge, outputnode, [('merged_file','epi_corrected')]) + pipeline.connect([ + (inputnode, dwell_time, [('epi_echospacing', 'dwell_time'), ( + 'pi_accel_factor', 'pi_factor'), ('epi_rev_encoding', 'is_reverse_encoding')]), (inputnode, select_mag, [('fieldmap_mag', 'in_file')]), (inputnode, norm_pha, [('fieldmap_pha', 'in_file')]), (select_mag, mask_mag, [('roi_file', 'in_file')]), (mask_mag, mask_mag_dil, [('mask_file', 'in_file')]), (select_mag, prelude, [('roi_file', 'magnitude_file')]), (norm_pha, prelude, [('out_file', 'phase_file')]), (mask_mag_dil, prelude, [('out_file', 'mask_file')]), (prelude, fill_phase, [('unwrapped_phase_file', 'in_file')]), (inputnode, vsm, [('fieldmap_mag', 'in_file')]), (fill_phase, vsm, [('out_file', 'phasemap_file')]), (inputnode, vsm, [(('te_diff', _ms2sec), 'asym_se_time'), ('vsm_sigma', 'smooth2d')]), (dwell_time, vsm, [(('dwell_time', _ms2sec), 'dwell_time')]), (mask_mag_dil, vsm, [('out_file', 'mask_file')]), (mask_mag_dil, vsm_mean, [('out_file', 'mask_file')]), (vsm, vsm_mean, [('unwarped_file', 'in_unwarped'), ('shift_out_file', 'in_file')]), (inputnode, dwi_split, [('in_file', 'in_file')]), (dwi_split, dwi_applyxfm, [('out_files', 'in_file')]), (dwi_applyxfm, dwi_merge, [('unwarped_file', 'in_files')]), (dwi_merge, outputnode, [('merged_file', 'epi_corrected')]) ]) if fieldmap_registration: """ Register magfw to example epi. There are some parameters here that may need to be tweaked. Should probably strip the mag Pre-condition: forward warp the mag in order to reg with func. What does mask do here? """ - # Select reference volume from EPI (B0 in dMRI and a middle frame in fMRI) - select_epi = pe.Node( interface=fsl.utils.ExtractROI(t_size=1), name="select_epi" ) - - # fugue -i %s -w %s --loadshift=%s --mask=%s % ( mag_name, magfw_name, vsmmag_name, mask_name ), log ) # Forward Map - vsm_fwd = pe.Node( interface=fsl.FUGUE(save_warped=True), name="vsm_fwd") - vsm_reg = pe.Node( interface=fsl.FLIRT( bins=256, cost='corratio', dof=6, interp='trilinear', searchr_x=[-10,10], searchr_y=[-10,10], searchr_z=[-10,10] ), name="vsm_registration" ) + # Select reference volume from EPI (B0 in dMRI and a middle frame in + # fMRI) + select_epi = pe.Node(interface=fsl.utils.ExtractROI( + t_size=1), name="select_epi") + + # fugue -i %s -w %s --loadshift=%s --mask=%s % ( mag_name, magfw_name, + # vsmmag_name, mask_name ), log ) # Forward Map + vsm_fwd = pe.Node(interface=fsl.FUGUE( + save_warped=True), name="vsm_fwd") + vsm_reg = pe.Node(interface=fsl.FLIRT(bins=256, cost='corratio', dof=6, interp='spline', searchr_x=[ + -10, 10], searchr_y=[-10, 10], searchr_z=[-10, 10]), name="vsm_registration") # 'flirt -in %s -ref %s -out %s -init %s -applyxfm' % ( vsmmag_name, ref_epi, vsmmag_name, magfw_mat_out ) - vsm_applyxfm = pe.Node( interface=fsl.ApplyXfm(), name='vsm_apply_xfm') + vsm_applyxfm = pe.Node(interface=fsl.ApplyXfm( + interp='spline'), name='vsm_apply_xfm') # 'flirt -in %s -ref %s -out %s -init %s -applyxfm' % ( mask_name, ref_epi, mask_name, magfw_mat_out ) - msk_applyxfm = pe.Node( interface=fsl.ApplyXfm(), name='msk_apply_xfm') + msk_applyxfm = pe.Node(interface=fsl.ApplyXfm( + interp='nearestneighbour'), name='msk_apply_xfm') pipeline.connect([ - (inputnode, select_epi, [('in_file','in_file'), ('ref_num','t_min')] ) - ,(select_epi, vsm_reg, [('roi_file','reference')]) - ,(vsm, vsm_fwd, [('shift_out_file','shift_in_file')]) - ,(mask_mag_dil, vsm_fwd, [('out_file','mask_file')] ) - ,(inputnode, vsm_fwd, [('fieldmap_mag', 'in_file')]) - ,(vsm_fwd, vsm_reg, [('warped_file','in_file')]) - ,(vsm_reg, msk_applyxfm, [('out_matrix_file','in_matrix_file')]) - ,(select_epi, msk_applyxfm, [('roi_file','reference')]) - ,(mask_mag_dil, msk_applyxfm, [('out_file','in_file') ]) - ,(vsm_reg, vsm_applyxfm, [('out_matrix_file','in_matrix_file')]) - ,(select_epi, vsm_applyxfm, [('roi_file','reference')]) - ,(vsm_mean, vsm_applyxfm, [('out_file','in_file')]) - ,(msk_applyxfm, dwi_applyxfm, [('out_file','mask_file')]) - ,(vsm_applyxfm, dwi_applyxfm, [('out_file','shift_in_file')]) + (inputnode, select_epi, [( + 'in_file', 'in_file'), ('ref_num', 't_min')]), (select_epi, vsm_reg, [('roi_file', 'reference')]), (vsm, vsm_fwd, [('shift_out_file', 'shift_in_file')]), (mask_mag_dil, vsm_fwd, [('out_file', 'mask_file')]), (inputnode, vsm_fwd, [('fieldmap_mag', 'in_file')]), (vsm_fwd, vsm_reg, [('warped_file', 'in_file')]), (vsm_reg, msk_applyxfm, [('out_matrix_file', 'in_matrix_file')]), (select_epi, msk_applyxfm, [('roi_file', 'reference')]), (mask_mag_dil, msk_applyxfm, [('out_file', 'in_file')]), (vsm_reg, vsm_applyxfm, [('out_matrix_file', 'in_matrix_file')]), (select_epi, vsm_applyxfm, [('roi_file', 'reference')]), (vsm_mean, vsm_applyxfm, [('out_file', 'in_file')]), (msk_applyxfm, dwi_applyxfm, [('out_file', 'mask_file')]), (vsm_applyxfm, dwi_applyxfm, [('out_file', 'shift_in_file')]) ]) else: pipeline.connect([ - (mask_mag_dil, dwi_applyxfm, [('out_file','mask_file')]) - ,(vsm_mean, dwi_applyxfm, [('out_file','shift_in_file')]) + (mask_mag_dil, dwi_applyxfm, [('out_file', 'mask_file')]), ( + vsm_mean, dwi_applyxfm, [('out_file', 'shift_in_file')]) ]) - return pipeline -def _rotate_bvecs( in_bvec, in_matrix ): + +def _rotate_bvecs(in_bvec, in_matrix): import os import numpy as np - name,fext = os.path.splitext( os.path.basename(in_bvec) ) - if fext == '.gz': name,_ = os.path.splitext(name) - out_file = os.path.abspath( './%s_rotated.bvec' % name ) - bvecs = np.loadtxt( in_bvec ) - new_bvecs = [ bvecs[:,0] ] - - for i,vol_matrix in enumerate(in_matrix[1::]): - bvec = np.matrix( bvecs[:,i] ) - rot = np.matrix(np.loadtxt( vol_matrix )[0:3,0:3]) - new_bvecs.append( (np.array( rot * bvec.T).T)[0] ) - np.savetxt( out_file, np.array(new_bvecs).T, fmt='%0.15f' ) + name, fext = os.path.splitext(os.path.basename(in_bvec)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('./%s_rotated.bvec' % name) + bvecs = np.loadtxt(in_bvec) + new_bvecs = [bvecs[:, 0]] + + for i, vol_matrix in enumerate(in_matrix[1::]): + bvec = np.matrix(bvecs[:, i]) + rot = np.matrix(np.loadtxt(vol_matrix)[0:3, 0:3]) + new_bvecs.append((np.array(rot * bvec.T).T)[0]) + np.savetxt(out_file, np.array(new_bvecs).T, fmt='%0.15f') return out_file -def _cat_logs( in_files ): + +def _cat_logs(in_files): import shutil import os - name,fext = os.path.splitext( os.path.basename(in_files[0]) ) - if fext == '.gz': name,_ = os.path.splitext(name) - out_file = os.path.abspath( './%s_ecclog.log' % name ) + name, fext = os.path.splitext(os.path.basename(in_files[0])) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('./%s_ecclog.log' % name) out_str = "" - with open( out_file, 'wb' ) as totallog: - for i,fname in enumerate(in_files): - totallog.write( "\n\npreprocessing %d\n" % i ) + with open(out_file, 'wb') as totallog: + for i, fname in enumerate(in_files): + totallog.write("\n\npreprocessing %d\n" % i) with open(fname) as inlog: for line in inlog: totallog.write(line) return out_file -def _compute_dwelltime( dwell_time=0.68, pi_factor=1.0, is_reverse_encoding=False ): - dwell_time*=(1.0/pi_factor) + +def _compute_dwelltime(dwell_time=0.68, pi_factor=1.0, is_reverse_encoding=False): + dwell_time *= (1.0/pi_factor) if is_reverse_encoding: - dwell_time*=-1.0 + dwell_time *= -1.0 return dwell_time -def _prepare_phasediff( in_file ): + +def _prepare_phasediff(in_file): import nibabel as nib import os import numpy as np - img = nib.load( in_file ) - max_diff = np.max( img.get_data().reshape(-1) ) - min_diff = np.min( img.get_data().reshape(-1) ) - A = ( 2.0 * np.pi )/( max_diff-min_diff ) - B = np.pi - ( A * max_diff ) + img = nib.load(in_file) + max_diff = np.max(img.get_data().reshape(-1)) + min_diff = np.min(img.get_data().reshape(-1)) + A = (2.0 * np.pi)/(max_diff-min_diff) + B = np.pi - (A * max_diff) diff_norm = img.get_data() * A + B - name,fext = os.path.splitext( os.path.basename(in_file) ) - if fext == '.gz': name,_ = os.path.splitext(name) - out_file = os.path.abspath( './%s_2pi.nii.gz' % name ) - nib.save( nib.Nifti1Image( diff_norm, img.get_affine(), img.get_header() ), out_file ) + name, fext = os.path.splitext(os.path.basename(in_file)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('./%s_2pi.nii.gz' % name) + nib.save(nib.Nifti1Image( + diff_norm, img.get_affine(), img.get_header()), out_file) return out_file -def _dilate_mask ( in_file, iterations=4 ): + +def _dilate_mask(in_file, iterations=4): import nibabel as nib import scipy.ndimage as ndimage import os - img = nib.load( in_file ) - img._data = ndimage.binary_dilation( img.get_data(), iterations=iterations ) - - name,fext = os.path.splitext( os.path.basename(in_file) ) - if fext == '.gz': name,_ = os.path.splitext(name) - out_file = os.path.abspath( './%s_dil.nii.gz' % name ) - nib.save( img, out_file ) + img = nib.load(in_file) + img._data = ndimage.binary_dilation(img.get_data(), iterations=iterations) + + name, fext = os.path.splitext(os.path.basename(in_file)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('./%s_dil.nii.gz' % name) + nib.save(img, out_file) return out_file - -def _fill_phase( in_file ): + +def _fill_phase(in_file): import nibabel as nib import os import numpy as np - img = nib.load( in_file ) - dumb_img = nib.Nifti1Image( np.zeros( img.get_shape() ), img.get_affine(), img.get_header() ) - out_nii = nib.funcs.concat_images( ( img, dumb_img ) ) - name,fext = os.path.splitext( os.path.basename(in_file) ) - if fext == '.gz': name,_ = os.path.splitext(name) - out_file = os.path.abspath( './%s_phase_unwrapped.nii.gz' % name ) - nib.save( out_nii, out_file ) + img = nib.load(in_file) + dumb_img = nib.Nifti1Image(np.zeros( + img.get_shape()), img.get_affine(), img.get_header()) + out_nii = nib.funcs.concat_images((img, dumb_img)) + name, fext = os.path.splitext(os.path.basename(in_file)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('./%s_phase_unwrapped.nii.gz' % name) + nib.save(out_nii, out_file) return out_file -def _vsm_remove_mean( in_file, mask_file, in_unwarped ): + +def _vsm_remove_mean(in_file, mask_file, in_unwarped): import nibabel as nib import os import numpy as np import numpy.ma as ma - img = nib.load( in_file ) - msk = nib.load( mask_file ).get_data() + img = nib.load(in_file) + msk = nib.load(mask_file).get_data() img_data = img.get_data() - img_data[ msk==0 ] = 0 - vsmmag_masked = ma.masked_values( img_data.reshape(-1), 0.0 ) + img_data[msk == 0] = 0 + vsmmag_masked = ma.masked_values(img_data.reshape(-1), 0.0) vsmmag_masked = vsmmag_masked - vsmmag_masked.mean() - img._data = vsmmag_masked.reshape( img.get_shape() ) - name,fext = os.path.splitext( os.path.basename(in_file) ) - if fext == '.gz': name,_ = os.path.splitext(name) - out_file = os.path.abspath( './%s_demeaned.nii.gz' % name ) - nib.save( img, out_file ) + img._data = vsmmag_masked.reshape(img.get_shape()) + name, fext = os.path.splitext(os.path.basename(in_file)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('./%s_demeaned.nii.gz' % name) + nib.save(img, out_file) return out_file -def _ms2sec( val ): + +def _ms2sec(val): return val*1e-3; -def _split_dwi( in_file ): + +def _split_dwi(in_file): import nibabel as nib import os out_files = [] - frames = nib.funcs.four_to_three( nib.load( in_file ) ) - name,fext = os.path.splitext( os.path.basename(in_file) ) - if fext == '.gz': name,_ = os.path.splitext(name) - for i,frame in enumerate(frames): - out_file = os.path.abspath( './%s_%03d.nii.gz' % (name,i) ) - nib.save( frame, out_file ) - out_files.append( out_file ) + frames = nib.funcs.four_to_three(nib.load(in_file)) + name, fext = os.path.splitext(os.path.basename(in_file)) + if fext == '.gz': + name, _ = os.path.splitext(name) + for i, frame in enumerate(frames): + out_file = os.path.abspath('./%s_%03d.nii.gz' % (name, i)) + nib.save(frame, out_file) + out_files.append(out_file) return out_files