diff --git a/examples/rsfmri_preprocessing.py b/examples/rsfmri_preprocessing.py new file mode 100644 index 0000000000..1cc8360ab6 --- /dev/null +++ b/examples/rsfmri_preprocessing.py @@ -0,0 +1,760 @@ +#!/usr/bin/env python +""" +A preprocessing workflow for Siemens resting state data. + +This workflow makes use of: + +- AFNI +- ANTS +- C3D_Affine_Tool +- DicomStack +- FreeSurfer +- FSL +- NiPy + +For example: + +python rsfmri_preprocessing.py -d /data/12345-34-1.dcm -f /data/Resting.nii + -s subj001 -n 2 --despike -o output + -p PBS --plugin_args "dict(qsub_args='-q many')" +""" + +import os + +from nipype import (ants, afni, fsl, freesurfer, nipy, Function, DataSink) +from nipype import Workflow, Node, MapNode + +from nipype.algorithms.rapidart import ArtifactDetect +from nipype.algorithms.misc import TSNR +from nipype.interfaces.fsl.utils import EPIDeWarp +from nipype.interfaces.io import FreeSurferSource +from nipype.interfaces.c3 import C3dAffineTool +from nipype.interfaces.utility import Merge, IdentityInterface +from nipype.utils.filemanip import filename_to_list + +import numpy as np +import scipy as sp +import nibabel as nb +from dcmstack.extract import default_extractor +from dicom import read_file + +imports = ['import os', + 'import nibabel as nb', + 'import numpy as np', + 'import scipy as sp', + 'from nipype.utils.filemanip import filename_to_list' + ] + + +def get_info(dicom_files): + """Given a Siemens dicom file return metadata + + Returns + ------- + RepetitionTime + Slice Acquisition Times + Spacing between slices + """ + meta = default_extractor(read_file(filename_to_list(dicom_files)[0], + stop_before_pixels=True, + force=True)) + return (meta['RepetitionTime']/1000., meta['CsaImage.MosaicRefAcqTimes'], + meta['SpacingBetweenSlices']) + + +def median(in_files): + """Computes an average of the median of each realigned timeseries + + Parameters + ---------- + + in_files: one or more realigned Nifti 4D time series + + Returns + ------- + + out_file: a 3D Nifti file + """ + + average = None + for idx, filename in enumerate(filename_to_list(in_files)): + img = nb.load(filename) + data = np.median(img.get_data(), axis=3) + if not average: + average = data + else: + average = average + data + median_img = nb.Nifti1Image(average/float(idx + 1), + img.get_affine(), img.get_header()) + filename = os.path.join(os.getcwd(), 'median.nii.gz') + median_img.to_filename(filename) + return filename + + +def get_aparc_aseg(files): + """Return the aparc+aseg.mgz file""" + for name in files: + if 'aparc+aseg.mgz' in name: + return name + raise ValueError('aparc+aseg.mgz not found') + + +def motion_regressors(motion_params, order=2, derivatives=2): + """Compute motion regressors upto given order and derivative + + motion + d(motion)/dt + d2(motion)/dt2 (linear + quadratic) + """ + out_files = [] + for idx, filename in enumerate(filename_to_list(motion_params)): + params = np.genfromtxt(filename) + out_params = params + for d in range(1, derivatives + 1): + cparams = np.vstack((np.repeat(params[0, :][None, :], d, axis=0), + params)) + out_params = np.hstack((out_params, np.diff(cparams, d, axis=0))) + out_params2 = out_params + for i in range(2, order + 1): + out_params2 = np.hstack((out_params2, np.power(out_params, i))) + filename = os.path.join(os.getcwd(), "motion_regressor%02d.txt" % idx) + np.savetxt(filename, out_params2, fmt="%.10f") + out_files.append(filename) + return out_files + + +def build_filter1(motion_params, comp_norm, outliers): + """Builds a regressor set comprisong motion parameters, composite norm and + outliers + + The outliers are added as a single time point column for each outlier + + + Parameters + ---------- + + motion_params: a text file containing motion parameters and its derivatives + comp_norm: a text file containing the composite norm + outliers: a text file containing 0-based outlier indices + + Returns + ------- + components_file: a text file containing all the regressors + + """ + out_files = [] + for idx, filename in enumerate(filename_to_list(motion_params)): + params = np.genfromtxt(filename) + norm_val = np.genfromtxt(filename_to_list(comp_norm)[idx]) + out_params = np.hstack((params, norm_val[:, None])) + try: + outlier_val = np.genfromtxt(filename_to_list(outliers)[idx]) + except IOError: + outlier_val = np.empty((0)) + for index in np.atleast_1d(outlier_val): + outlier_vector = np.zeros((out_params.shape[0], 1)) + outlier_vector[index] = 1 + out_params = np.hstack((out_params, outlier_vector)) + filename = os.path.join(os.getcwd(), "filter_regressor%02d.txt" % idx) + np.savetxt(filename, out_params, fmt="%.10f") + out_files.append(filename) + return out_files + + +def extract_noise_components(realigned_file, mask_file, num_components=6): + """Derive components most reflective of physiological noise + + Parameters + ---------- + realigned_file: a 4D Nifti file containing realigned volumes + mask_file: a 3D Nifti file containing white matter + ventricular masks + num_components: number of components to use for noise decomposition + + Returns + ------- + components_file: a text file containing the noise components + """ + imgseries = nb.load(realigned_file) + noise_mask = nb.load(mask_file) + voxel_timecourses = imgseries.get_data()[np.nonzero(noise_mask.get_data())] + voxel_timecourses = voxel_timecourses.byteswap().newbyteorder() + voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0 + _, _, v = sp.linalg.svd(voxel_timecourses, full_matrices=False) + components_file = os.path.join(os.getcwd(), 'noise_components.txt') + np.savetxt(components_file, v[:num_components, :].T) + return components_file + + +def extract_subrois(timeseries_file, label_file, indices): + """Extract voxel time courses for each subcortical roi index + + Parameters + ---------- + + timeseries_file: a 4D Nifti file + label_file: a 3D file containing rois in the same space/size of the 4D file + indices: a list of indices for ROIs to extract. Currently a dictionary + mapping freesurfer indices to CMA/Label Fusion indices are being used + + Returns + ------- + out_file: a text file containing time courses for each voxel of each roi + The first four columns are: freesurfer index, i, j, k positions in the + label file + """ + img = nb.load(timeseries_file) + data = img.get_data() + roiimg = nb.load(label_file) + rois = roiimg.get_data() + out_ts_file = os.path.join(os.getcwd(), 'subcortical_timeseries.txt') + with open(out_ts_file, 'wt') as fp: + for fsindex, cmaindex in sorted(indices.items()): + ijk = np.nonzero(rois == cmaindex) + ts = data[ijk] + for i0, row in enumerate(ts): + fp.write('%d,%d,%d,%d,' % (fsindex, ijk[0][i0], + ijk[1][i0], ijk[2][i0]) + + ','.join(['%.10f' % val for val in row]) + '\n') + return out_ts_file + + +def combine_hemi(left, right): + """Combine left and right hemisphere time series into a single text file + """ + lh_data = nb.load(left).get_data() + rh_data = nb.load(right).get_data() + + indices = np.vstack((1000000 + np.arange(0, lh_data.shape[0])[:, None], + 2000000 + np.arange(0, rh_data.shape[0])[:, None])) + all_data = np.hstack((indices, np.vstack((lh_data.squeeze(), + rh_data.squeeze())))) + filename = 'combined_surf.txt' + np.savetxt(filename, all_data, + fmt=','.join(['%d'] + ['%.10f'] * (all_data.shape[1] - 1))) + return os.path.abspath(filename) + + +""" +Creates the main preprocessing workflow +""" + + +def create_workflow(files, + subject_id, + n_vol=0, + despike=True, + TR=None, + slice_times=None, + slice_thickness=None, + fieldmap_images=[], + norm_threshold=1, + num_components=6, + vol_fwhm=None, + surf_fwhm=None, + lowpass_freq=-1, + highpass_freq=-1, + sink_directory=os.getcwd(), + FM_TEdiff=2.46, + FM_sigma=2, + FM_echo_spacing=.7, + target_subject=['fsaverage3', 'fsaverage4'], + name='resting'): + + wf = Workflow(name=name) + + # Skip starting volumes + remove_vol = MapNode(fsl.ExtractROI(t_min=n_vol, t_size=-1), + iterfield=['in_file'], + name="remove_volumes") + remove_vol.inputs.in_file = files + + # Run AFNI's despike. This is always run, however, whether this is fed to + # realign depends on the input configuration + despiker = MapNode(afni.Despike(outputtype='NIFTI_GZ'), + iterfield=['in_file'], + name='despike') + #despiker.plugin_args = {'qsub_args': '-l nodes=1:ppn='} + + wf.connect(remove_vol, 'roi_file', despiker, 'in_file') + + # Run Nipy joint slice timing and realignment algorithm + realign = Node(nipy.SpaceTimeRealigner(), name='realign') + realign.inputs.tr = TR + realign.inputs.slice_times = slice_times + realign.inputs.slice_info = 2 + + if despike: + wf.connect(despiker, 'out_file', realign, 'in_file') + else: + wf.connect(remove_vol, 'roi_file', realign, 'in_file') + + # Comute TSNR on realigned data regressing polynomials upto order 2 + tsnr = MapNode(TSNR(regress_poly=2), iterfield=['in_file'], name='tsnr') + wf.connect(realign, 'out_file', tsnr, 'in_file') + + # Compute the median image across runs + calc_median = Node(Function(input_names=['in_files'], + output_names=['median_file'], + function=median, + imports=imports), + name='median') + wf.connect(tsnr, 'detrended_file', calc_median, 'in_files') + + # Coregister the median to the surface + register = Node(freesurfer.BBRegister(), + name='bbregister') + register.inputs.subject_id = subject_id + register.inputs.init = 'fsl' + register.inputs.contrast_type = 't2' + register.inputs.out_fsl_file = True + register.inputs.epi_mask = True + + # Compute fieldmaps and unwarp using them + if fieldmap_images: + fieldmap = Node(interface=EPIDeWarp(), name='fieldmap_unwarp') + fieldmap.inputs.tediff = FM_TEdiff + fieldmap.inputs.esp = FM_echo_spacing + fieldmap.inputs.sigma = FM_sigma + fieldmap.inputs.mag_file = fieldmap_images[0] + fieldmap.inputs.dph_file = fieldmap_images[1] + wf.connect(calc_median, 'median_file', fieldmap, 'exf_file') + + dewarper = MapNode(interface=fsl.FUGUE(), iterfield=['in_file'], + name='dewarper') + wf.connect(tsnr, 'detrended_file', dewarper, 'in_file') + wf.connect(fieldmap, 'exf_mask', dewarper, 'mask_file') + wf.connect(fieldmap, 'vsm_file', dewarper, 'shift_in_file') + wf.connect(fieldmap, 'exfdw', register, 'source_file') + else: + wf.connect(calc_median, 'median_file', register, 'source_file') + + # Get the subject's freesurfer source directory + fssource = Node(FreeSurferSource(), + name='fssource') + fssource.inputs.subject_id = subject_id + fssource.inputs.subjects_dir = os.environ['SUBJECTS_DIR'] + + # Extract wm+csf, brain masks by eroding freesurfer lables and then + # transform the masks into the space of the median + wmcsf = Node(freesurfer.Binarize(), name='wmcsfmask') + mask = wmcsf.clone('anatmask') + wmcsftransform = Node(freesurfer.ApplyVolTransform(inverse=True, + interp='nearest'), + name='wmcsftransform') + wmcsftransform.inputs.subjects_dir = os.environ['SUBJECTS_DIR'] + wmcsf.inputs.wm_ven_csf = True + wmcsf.inputs.match = [4, 5, 14, 15, 24, 31, 43, 44, 63] + wmcsf.inputs.binary_file = 'wmcsf.nii.gz' + wmcsf.inputs.erode = int(np.ceil(slice_thickness)) + wf.connect(fssource, ('aparc_aseg', get_aparc_aseg), wmcsf, 'in_file') + if fieldmap_images: + wf.connect(fieldmap, 'exf_mask', wmcsftransform, 'source_file') + else: + wf.connect(calc_median, 'median_file', wmcsftransform, 'source_file') + wf.connect(register, 'out_reg_file', wmcsftransform, 'reg_file') + wf.connect(wmcsf, 'binary_file', wmcsftransform, 'target_file') + + mask.inputs.binary_file = 'mask.nii.gz' + mask.inputs.dilate = int(np.ceil(slice_thickness)) + 1 + mask.inputs.erode = int(np.ceil(slice_thickness)) + mask.inputs.min = 0.5 + wf.connect(fssource, ('aparc_aseg', get_aparc_aseg), mask, 'in_file') + masktransform = wmcsftransform.clone("masktransform") + if fieldmap_images: + wf.connect(fieldmap, 'exf_mask', masktransform, 'source_file') + else: + wf.connect(calc_median, 'median_file', masktransform, 'source_file') + wf.connect(register, 'out_reg_file', masktransform, 'reg_file') + wf.connect(mask, 'binary_file', masktransform, 'target_file') + + # Compute Art outliers + art = Node(interface=ArtifactDetect(use_differences=[True, False], + use_norm=True, + norm_threshold=norm_threshold, + zintensity_threshold=3, + parameter_source='NiPy', + bound_by_brainmask=True, + save_plot=False, + mask_type='file'), + name="art") + if fieldmap_images: + wf.connect(dewarper, 'unwarped_file', art, 'realigned_files') + else: + wf.connect(tsnr, 'detrended_file', art, 'realigned_files') + wf.connect(realign, 'par_file', + art, 'realignment_parameters') + wf.connect(masktransform, 'transformed_file', art, 'mask_file') + + # Compute motion regressors + motreg = Node(Function(input_names=['motion_params', 'order', + 'derivatives'], + output_names=['out_files'], + function=motion_regressors, + imports=imports), + name='getmotionregress') + wf.connect(realign, 'par_file', motreg, 'motion_params') + + # Create a filter to remove motion and art confounds + createfilter1 = Node(Function(input_names=['motion_params', 'comp_norm', + 'outliers'], + output_names=['out_files'], + function=build_filter1, + imports=imports), + name='makemotionbasedfilter') + wf.connect(motreg, 'out_files', createfilter1, 'motion_params') + wf.connect(art, 'norm_files', createfilter1, 'comp_norm') + wf.connect(art, 'outlier_files', createfilter1, 'outliers') + + # Filter the motion and art confounds + filter1 = MapNode(fsl.GLM(out_res_name='timeseries.nii.gz', + demean=True), + iterfield=['in_file', 'design'], + name='filtermotion') + if fieldmap_images: + wf.connect(dewarper, 'unwarped_file', filter1, 'in_file') + else: + wf.connect(tsnr, 'detrended_file', filter1, 'in_file') + wf.connect(createfilter1, 'out_files', filter1, 'design') + wf.connect(masktransform, 'transformed_file', filter1, 'mask') + + # Create a filter to remove noise components based on white matter and CSF + createfilter2 = MapNode(Function(input_names=['realigned_file', 'mask_file', + 'num_components'], + output_names=['out_files'], + function=extract_noise_components, + imports=imports), + iterfield=['realigned_file'], + name='makecompcorrfilter') + createfilter2.inputs.num_components = num_components + wf.connect(filter1, 'out_res', createfilter2, 'realigned_file') + wf.connect(masktransform, 'transformed_file', createfilter2, 'mask_file') + + # Filter noise components + filter2 = MapNode(fsl.GLM(out_res_name='timeseries_cleaned.nii.gz', + demean=True), + iterfield=['in_file', 'design'], + name='filtercompcorr') + wf.connect(filter1, 'out_res', filter2, 'in_file') + wf.connect(createfilter2, 'out_files', filter2, 'design') + wf.connect(masktransform, 'transformed_file', filter2, 'mask') + + # Smoothing using surface and volume smoothing + smooth = MapNode(freesurfer.Smooth(), + iterfield=['in_file'], + name='smooth') + smooth.inputs.proj_frac_avg = (0.1, 0.9, 0.1) + if surf_fwhm is None: + surf_fwhm = 5 * slice_thickness + smooth.inputs.surface_fwhm = surf_fwhm + if vol_fwhm is None: + vol_fwhm = 2 * slice_thickness + smooth.inputs.vol_fwhm = vol_fwhm + wf.connect(filter2, 'out_res', smooth, 'in_file') + wf.connect(register, 'out_reg_file', smooth, 'reg_file') + + # Bandpass filter the data + bandpass = MapNode(fsl.TemporalFilter(), + iterfield=['in_file'], + name='bandpassfilter') + if highpass_freq < 0: + bandpass.inputs.highpass_sigma = -1 + else: + bandpass.inputs.highpass_sigma = 1. / (2 * TR * highpass_freq) + if lowpass_freq < 0: + bandpass.inputs.lowpass_sigma = -1 + else: + bandpass.inputs.lowpass_sigma = 1. / (2 * TR * lowpass_freq) + wf.connect(smooth, 'smoothed_file', bandpass, 'in_file') + + # Convert aparc to subject functional space + aparctransform = wmcsftransform.clone("aparctransform") + if fieldmap_images: + wf.connect(fieldmap, 'exf_mask', aparctransform, 'source_file') + else: + wf.connect(calc_median, 'median_file', aparctransform, 'source_file') + wf.connect(register, 'out_reg_file', aparctransform, 'reg_file') + wf.connect(fssource, ('aparc_aseg', get_aparc_aseg), + aparctransform, 'target_file') + + # Sample the average time series in aparc ROIs + sampleaparc = MapNode(freesurfer.SegStats(avgwf_txt_file=True, + default_color_table=True), + iterfield=['in_file'], + name='aparc_ts') + sampleaparc.inputs.segment_id = ([8] + range(10, 14) + [17, 18, 26, 47] + + range(49, 55) + [58] + range(1001, 1036) + + range(2001, 2036)) + + wf.connect(aparctransform, 'transformed_file', + sampleaparc, 'segmentation_file') + wf.connect(bandpass, 'out_file', sampleaparc, 'in_file') + + # Sample the time series onto the surface of the target surface. Performs + # sampling into left and right hemisphere + target = Node(IdentityInterface(fields=['target_subject']), name='target') + target.iterables = ('target_subject', filename_to_list(target_subject)) + + samplerlh = MapNode(freesurfer.SampleToSurface(), + iterfield=['source_file'], + name='sampler_lh') + samplerlh.inputs.sampling_method = "average" + samplerlh.inputs.sampling_range = (0.1, 0.9, 0.1) + samplerlh.inputs.sampling_units = "frac" + samplerlh.inputs.interp_method = "trilinear" + #samplerlh.inputs.cortex_mask = True + samplerlh.inputs.out_type = 'niigz' + samplerlh.inputs.subjects_dir = os.environ['SUBJECTS_DIR'] + + samplerrh = samplerlh.clone('sampler_rh') + + samplerlh.inputs.hemi = 'lh' + wf.connect(bandpass, 'out_file', samplerlh, 'source_file') + wf.connect(register, 'out_reg_file', samplerlh, 'reg_file') + wf.connect(target, 'target_subject', samplerlh, 'target_subject') + + samplerrh.set_input('hemi', 'rh') + wf.connect(bandpass, 'out_file', samplerrh, 'source_file') + wf.connect(register, 'out_reg_file', samplerrh, 'reg_file') + wf.connect(target, 'target_subject', samplerrh, 'target_subject') + + # Combine left and right hemisphere to text file + combiner = MapNode(Function(input_names=['left', 'right'], + output_names=['out_file'], + function=combine_hemi, + imports=imports), + iterfield=['left', 'right'], + name="combiner") + wf.connect(samplerlh, 'out_file', combiner, 'left') + wf.connect(samplerrh, 'out_file', combiner, 'right') + + # Compute registration between the subject's structural and MNI template + # This is currently set to perform a very quick registration. However, the + # registration can be made significantly more accurate for cortical + # structures by increasing the number of iterations + # All parameters are set using the example from: + # + reg = Node(ants.Registration(), name='antsRegister') + reg.inputs.output_transform_prefix = "output_" + reg.inputs.transforms = ['Translation', 'Rigid', 'Affine', 'SyN'] + reg.inputs.transform_parameters = [(0.1,), (0.1,), (0.1,), (0.2, 3.0, 0.0)] + # reg.inputs.number_of_iterations = ([[10000, 111110, 11110]]*3 + + # [[100, 50, 30]]) + reg.inputs.number_of_iterations = [[100, 100, 100]] * 3 + [[100, 20, 10]] + reg.inputs.dimension = 3 + reg.inputs.write_composite_transform = True + reg.inputs.collapse_output_transforms = False + reg.inputs.metric = ['Mattes'] * 3 + [['Mattes', 'CC']] + reg.inputs.metric_weight = [1] * 3 + [[0.5, 0.5]] + reg.inputs.radius_or_number_of_bins = [32] * 3 + [[32, 4]] + reg.inputs.sampling_strategy = ['Regular'] * 3 + [[None, None]] + reg.inputs.sampling_percentage = [0.3] * 3 + [[None, None]] + reg.inputs.convergence_threshold = [1.e-8] * 3 + [-0.01] + reg.inputs.convergence_window_size = [20] * 3 + [5] + reg.inputs.smoothing_sigmas = [[4, 2, 1]] * 3 + [[1, 0.5, 0]] + reg.inputs.sigma_units = ['vox'] * 4 + reg.inputs.shrink_factors = [[6, 4, 2]] + [[3, 2, 1]]*2 + [[4, 2, 1]] + reg.inputs.use_estimate_learning_rate_once = [True] * 4 + reg.inputs.use_histogram_matching = [False] * 3 + [True] + reg.inputs.output_warped_image = 'output_warped_image.nii.gz' + reg.inputs.fixed_image = \ + os.path.abspath('OASIS-TRT-20_template_to_MNI152_2mm.nii.gz') + reg.inputs.num_threads = 4 + reg.inputs.terminal_output = 'file' + reg.plugin_args = {'qsub_args': '-l nodes=1:ppn=4'} + + # Convert T1.mgz to nifti for using with ANTS + convert = Node(freesurfer.MRIConvert(out_type='niigz'), name='convert2nii') + wf.connect(fssource, 'T1', convert, 'in_file') + + # Mask the T1.mgz file with the brain mask computed earlier + maskT1 = Node(fsl.BinaryMaths(operation='mul'), name='maskT1') + wf.connect(mask, 'binary_file', maskT1, 'operand_file') + wf.connect(convert, 'out_file', maskT1, 'in_file') + wf.connect(maskT1, 'out_file', reg, 'moving_image') + + # Convert the BBRegister transformation to ANTS ITK format + convert2itk = MapNode(C3dAffineTool(), + iterfield=['transform_file', 'source_file'], + name='convert2itk') + convert2itk.inputs.fsl2ras = True + convert2itk.inputs.itk_transform = True + wf.connect(register, 'out_fsl_file', convert2itk, 'transform_file') + if fieldmap_images: + wf.connect(fieldmap, 'exf_mask', convert2itk, 'source_file') + else: + wf.connect(calc_median, 'median_file', convert2itk, 'source_file') + wf.connect(convert, 'out_file', convert2itk, 'reference_file') + + # Concatenate the affine and ants transforms into a list + pickfirst = lambda x: x[0] + merge = MapNode(Merge(2), iterfield=['in2'], name='mergexfm') + wf.connect(convert2itk, 'itk_transform', merge, 'in2') + wf.connect(reg, ('composite_transform', pickfirst), merge, 'in1') + + # Apply the combined transform to the time series file + sample2mni = MapNode(ants.ApplyTransforms(), + iterfield=['input_image', 'transforms'], + name='sample2mni') + sample2mni.inputs.input_image_type = 3 + sample2mni.inputs.interpolation = 'BSpline' + sample2mni.inputs.invert_transform_flags = [False, False] + sample2mni.inputs.reference_image = \ + os.path.abspath('OASIS-TRT-20_template_to_MNI152_2mm.nii.gz') + sample2mni.inputs.terminal_output = 'file' + wf.connect(bandpass, 'out_file', sample2mni, 'input_image') + wf.connect(merge, 'out', sample2mni, 'transforms') + + # Sample the time series file for each subcortical roi + ts2txt = MapNode(Function(input_names=['timeseries_file', 'label_file', + 'indices'], + output_names=['out_file'], + function=extract_subrois, + imports=imports), + iterfield=['timeseries_file'], + name='getsubcortts') + ts2txt.inputs.indices = dict(zip([8] + range(10, 14) + [17, 18, 26, 47] + + range(49, 55) + [58], + [39, 60, 37, 58, 56, 48, 32, 30, + 38, 59, 36, 57, 55, 47, 31, 23])) + ts2txt.inputs.label_file = \ + os.path.abspath(('OASIS-TRT-20_DKT31_CMA_jointfusion_labels_in_MNI152' + '_2mm.nii.gz')) + wf.connect(sample2mni, 'output_image', ts2txt, 'timeseries_file') + + # Save the relevant data into an output directory + datasink = Node(interface=DataSink(), name="datasink") + datasink.inputs.base_directory = sink_directory + datasink.inputs.container = subject_id + datasink.inputs.substitutions = [('_target_subject_', '')] + datasink.inputs.regexp_substitutions = (r'(/_.*(\d+/))', r'/run\2') + wf.connect(despiker, 'out_file', datasink, 'resting.qa.despike') + wf.connect(realign, 'par_file', datasink, 'resting.qa.motion') + wf.connect(tsnr, 'tsnr_file', datasink, 'resting.qa.tsnr') + wf.connect(tsnr, 'mean_file', datasink, 'resting.qa.tsnr.@mean') + wf.connect(tsnr, 'stddev_file', datasink, 'resting.qa.@tsnr_stddev') + if fieldmap_images: + wf.connect(fieldmap, 'exf_mask', datasink, 'resting.reference') + else: + wf.connect(calc_median, 'median_file', datasink, 'resting.reference') + wf.connect(art, 'norm_files', datasink, 'resting.qa.art.@norm') + wf.connect(art, 'intensity_files', datasink, 'resting.qa.art.@intensity') + wf.connect(art, 'outlier_files', datasink, 'resting.qa.art.@outlier_files') + wf.connect(mask, 'binary_file', datasink, 'resting.mask') + wf.connect(masktransform, 'transformed_file', + datasink, 'resting.mask.@transformed_file') + wf.connect(register, 'out_reg_file', datasink, 'resting.registration.bbreg') + wf.connect(reg, ('composite_transform', pickfirst), + datasink, 'resting.registration.ants') + wf.connect(register, 'min_cost_file', + datasink, 'resting.qa.bbreg.@mincost') + wf.connect(smooth, 'smoothed_file', datasink, 'resting.timeseries.fullpass') + wf.connect(bandpass, 'out_file', datasink, 'resting.timeseries.bandpassed') + wf.connect(sample2mni, 'output_image', datasink, 'resting.timeseries.mni') + wf.connect(createfilter1, 'out_files', + datasink, 'resting.regress.@regressors') + wf.connect(createfilter2, 'out_files', + datasink, 'resting.regress.@compcorr') + wf.connect(sampleaparc, 'summary_file', + datasink, 'resting.parcellations.aparc') + wf.connect(sampleaparc, 'avgwf_txt_file', + datasink, 'resting.parcellations.aparc.@avgwf') + wf.connect(ts2txt, 'out_file', + datasink, 'resting.parcellations.grayo.@subcortical') + datasink2 = Node(interface=DataSink(), name="datasink2") + datasink2.inputs.base_directory = sink_directory + datasink2.inputs.container = subject_id + datasink2.inputs.substitutions = [('_target_subject_', '')] + datasink2.inputs.regexp_substitutions = (r'(/_.*(\d+/))', r'/run\2') + wf.connect(combiner, 'out_file', + datasink2, 'resting.parcellations.grayo.@surface') + return wf + + +""" +Creates the full workflow including getting information from dicom files +""" + + +def create_resting_workflow(args): + TR = args.TR + slice_times = args.slice_times + slice_thickness = None + if args.dicom_file: + TR, slice_times, slice_thickness = get_info(args.dicom_file) + slice_times = (np.array(slice_times)/1000.).tolist() + + if slice_thickness is None: + from nibabel import load + img = load(args.files[0]) + slice_thickness = max(img.get_header().get_zooms()[:3]) + + kwargs = dict(files=[os.path.abspath(filename) for + filename in args.files], + subject_id=args.subject_id, + n_vol=args.n_vol, + despike=args.despike, + TR=TR, + slice_times=slice_times, + slice_thickness=slice_thickness, + lowpass_freq=args.lowpass_freq, + highpass_freq=args.highpass_freq, + sink_directory=os.path.abspath(args.sink)) + if args.field_maps: + kwargs.update(**dict(fieldmap_images=args.field_maps, + FM_TEdiff=args.TE_diff, + FM_echo_spacing=args.echo_spacing, + FM_sigma=args.sigma)) + wf = create_workflow(**kwargs) + return wf + +if __name__ == "__main__": + from argparse import ArgumentParser + parser = ArgumentParser(description=__doc__) + parser.add_argument("-d", "--dicom_file", dest="dicom_file", + help="an example dicom file from the resting series") + parser.add_argument("-f", "--files", dest="files", nargs="+", + help="4d nifti files for resting state", + required=True) + parser.add_argument("-s", "--subject_id", dest="subject_id", + help="FreeSurfer subject id", required=True) + parser.add_argument("-n", "--n_vol", dest="n_vol", default=0, type=int, + help="Volumes to skip at the beginning") + parser.add_argument("--despike", dest="despike", default=False, + action="store_true", help="Use despiked data") + parser.add_argument("--TR", dest="TR", default=None, + help="TR if dicom not provided in seconds") + parser.add_argument("--slice_times", dest="slice_times", nargs="+", + type=float, help="Slice times in seconds") + parser.add_argument("-l", "--lowpass_freq", dest="lowpass_freq", + default=-1, help="Low pass frequency (Hz)") + parser.add_argument("-u", "--highpass_freq", dest="highpass_freq", + default=-1, help="High pass frequency (Hz)") + parser.add_argument("-o", "--output_dir", dest="sink", + help="Output directory base") + parser.add_argument("-w", "--work_dir", dest="work_dir", + help="Output directory base") + parser.add_argument("-p", "--plugin", dest="plugin", + default='Linear', + help="Plugin to use") + parser.add_argument("--plugin_args", dest="plugin_args", + help="Plugin arguments") + parser.add_argument("--field_maps", dest="field_maps", nargs="+", + help="field map niftis") + parser.add_argument("--fm_echospacing", dest="echo_spacing", type=float, + help="field map echo spacing") + parser.add_argument("--fm_TE_diff", dest='TE_diff', type=float, + help="field map echo time difference") + parser.add_argument("--fm_sigma", dest='sigma', type=float, + help="field map sigma value") + args = parser.parse_args() + + wf = create_resting_workflow(args) + + if args.work_dir: + work_dir = os.path.abspath(args.work_dir) + else: + work_dir = os.getcwd() + + wf.base_dir = work_dir + if args.plugin_args: + wf.run(args.plugin, plugin_args=eval(args.plugin_args)) + else: + wf.run(args.plugin) diff --git a/nipype/interfaces/freesurfer/model.py b/nipype/interfaces/freesurfer/model.py index 9f35ae8d0e..19963f3832 100644 --- a/nipype/interfaces/freesurfer/model.py +++ b/nipype/interfaces/freesurfer/model.py @@ -342,9 +342,9 @@ def __init__(self, **kwargs): class BinarizeInputSpec(FSTraitedSpec): in_file = File(exists=True, argstr='--i %s', mandatory=True, copyfile=False, desc='input volume') - min = traits.Float(argstr='--min %f', + min = traits.Float(argstr='--min %f', xor=['wm_ven_csf'], desc='min thresh') - max = traits.Float(argstr='--max %f', + max = traits.Float(argstr='--max %f', xor=['wm_ven_csf'], desc='max thresh') rmin = traits.Float(argstr='--rmin %f', desc='compute min based on rmin*globalmean') @@ -356,7 +356,7 @@ class BinarizeInputSpec(FSTraitedSpec): desc='set match vals to 2 and 41 (aseg for cerebral WM)') ventricles = traits.Bool(argstr='--ventricles', desc='set match vals those for aseg ventricles+choroid (not 4th)') - wm_ven_csf = traits.Bool(argstr='--wm+vcsf', + wm_ven_csf = traits.Bool(argstr='--wm+vcsf', xor=['min', 'max'], desc='WM and ventricular CSF, including choroid (not 4th)') binary_file = File(argstr='--o %s', genfile=True, desc='binary output volume') @@ -430,7 +430,7 @@ def _list_outputs(self): outfile = fname_presuffix(self.inputs.in_file, newpath=os.getcwd(), suffix='_thresh') - outputs['binary_file'] = outfile + outputs['binary_file'] = os.path.abspath(outfile) value = self.inputs.count_file if isdefined(value): if isinstance(value, bool): diff --git a/nipype/interfaces/fsl/model.py b/nipype/interfaces/fsl/model.py index 273be0c662..b883f8d4be 100755 --- a/nipype/interfaces/fsl/model.py +++ b/nipype/interfaces/fsl/model.py @@ -1804,7 +1804,6 @@ def _list_outputs(self): return outputs def _gen_filename(self, name): - if name in ('out_file'): - return fname_presuffix(self.inputs.in_file, - suffix='_glm.txt', use_ext=False) + if name in ['out_file']: + return self._gen_fname(self.inputs.in_file, suffix='_glm') return None diff --git a/nipype/interfaces/io.py b/nipype/interfaces/io.py index 07a53b831b..e70c1349c6 100644 --- a/nipype/interfaces/io.py +++ b/nipype/interfaces/io.py @@ -214,7 +214,7 @@ class DataSink(IOBase): input_spec = DataSinkInputSpec output_spec = DataSinkOutputSpec - def __init__(self, infields=None, **kwargs): + def __init__(self, infields=None, force_run=True, **kwargs): """ Parameters ---------- @@ -232,6 +232,8 @@ def __init__(self, infields=None, **kwargs): self.inputs._outputs[key] = Undefined undefined_traits[key] = Undefined self.inputs.trait_set(trait_change_notify=False, **undefined_traits) + if force_run: + self._always_run = True def _get_dst(self, src): ## If path is directory with trailing os.path.sep, diff --git a/nipype/interfaces/nipy/__init__.py b/nipype/interfaces/nipy/__init__.py index ec541b443f..d78ae64447 100644 --- a/nipype/interfaces/nipy/__init__.py +++ b/nipype/interfaces/nipy/__init__.py @@ -1,3 +1,3 @@ from .model import FitGLM, EstimateContrast -from .preprocess import ComputeMask, FmriRealign4d +from .preprocess import ComputeMask, FmriRealign4d, SpaceTimeRealigner from .utils import Similarity diff --git a/nipype/interfaces/nipy/preprocess.py b/nipype/interfaces/nipy/preprocess.py index 46f514f9da..b0c12ea9e0 100644 --- a/nipype/interfaces/nipy/preprocess.py +++ b/nipype/interfaces/nipy/preprocess.py @@ -24,6 +24,8 @@ from nipy.labs.mask import compute_mask from nipy.algorithms.registration import FmriRealign4d as FR4d from nipy import save_image, load_image + from nipy.algorithms.registration import SpaceTimeRealign + from nipy.algorithms.registration.groupwise_registration import SpaceRealign from ..base import (TraitedSpec, BaseInterface, traits, BaseInterfaceInputSpec, isdefined, File, @@ -82,7 +84,7 @@ class FmriRealign4dInputSpec(BaseInterfaceInputSpec): desc="File to realign") tr = traits.Float(desc="TR in seconds", mandatory=True) - slice_order = traits.List(traits.Int(), maxver=0.3, + slice_order = traits.List(traits.Int(), desc=('0 based slice order. This would be equivalent to entering' 'np.argsort(spm_slice_order) for this field. This effects' 'interleaved acquisition. This field will be deprecated in' @@ -117,6 +119,7 @@ class FmriRealign4dOutputSpec(TraitedSpec): desc="Motion parameter files") +@np.deprecate_with_doc('Please use SpaceTimeRealign instead') class FmriRealign4d(BaseInterface): """Simultaneous motion and slice timing correction algorithm @@ -196,6 +199,131 @@ def _list_outputs(self): return outputs +class SpaceTimeRealignerInputSpec(BaseInterfaceInputSpec): + + in_file = InputMultiPath(exists=True, + mandatory=True, + desc="File to realign") + tr = traits.Float(desc="TR in seconds", requires=['slice_times']) + slice_times = traits.Either(traits.List(traits.Float()), + traits.Enum('asc_alt_2', 'asc_alt_2_1', + 'asc_alt_half', 'asc_alt_siemens', + 'ascending', 'desc_alt_2', + 'desc_alt_half', 'descending'), + desc=('Actual slice acquisition times.')) + slice_info = traits.Either(traits.Int, + traits.List(min_len=2, max_len=2), + desc=('Single integer or length 2 sequence ' + 'If int, the axis in `images` that is the ' + 'slice axis. In a 4D image, this will ' + 'often be axis = 2. If a 2 sequence, then' + ' elements are ``(slice_axis, ' + 'slice_direction)``, where ``slice_axis`` ' + 'is the slice axis in the image as above, ' + 'and ``slice_direction`` is 1 if the ' + 'slices were acquired slice 0 first, slice' + ' -1 last, or -1 if acquired slice -1 ' + 'first, slice 0 last. If `slice_info` is ' + 'an int, assume ' + '``slice_direction`` == 1.'), + requires=['slice_times'], + ) + + +class SpaceTimeRealignerOutputSpec(TraitedSpec): + out_file = OutputMultiPath(File(exists=True), + desc="Realigned files") + par_file = OutputMultiPath(File(exists=True), + desc=("Motion parameter files. Angles are not " + "euler angles")) + + +class SpaceTimeRealigner(BaseInterface): + """Simultaneous motion and slice timing correction algorithm + + If slice_times is not specified, this algorithm performs spatial motion + correction + + This interface wraps nipy's SpaceTimeRealign algorithm [1]_ or simply the + SpatialRealign algorithm when timing info is not provided. + + Examples + -------- + >>> from nipype.interfaces.nipy import SpaceTimeRealigner + >>> #Run spatial realignment only + >>> realigner = SpaceTimeRealigner() + >>> realigner.inputs.in_file = ['functional.nii'] + >>> res = realigner.run() # doctest: +SKIP + + >>> realigner = SpaceTimeRealigner() + >>> realigner.inputs.in_file = ['functional.nii'] + >>> realigner.inputs.tr = 2 + >>> realigner.inputs.slice_times = range(0, 3, 67) + >>> realigner.inputs.slice_info = 2 + >>> res = realigner.run() # doctest: +SKIP + + + References + ---------- + .. [1] Roche A. A four-dimensional registration algorithm with \ + application to joint correction of motion and slice timing \ + in fMRI. IEEE Trans Med Imaging. 2011 Aug;30(8):1546-54. DOI_. + + .. _DOI: http://dx.doi.org/10.1109/TMI.2011.2131152 + + """ + + input_spec = SpaceTimeRealignerInputSpec + output_spec = SpaceTimeRealignerOutputSpec + keywords = ['slice timing', 'motion correction'] + + def _run_interface(self, runtime): + + all_ims = [load_image(fname) for fname in self.inputs.in_file] + + if not isdefined(self.inputs.slice_times): + R = SpaceRealign(all_ims) + else: + R = SpaceTimeRealign(all_ims, + tr=self.inputs.tr, + slice_times=self.inputs.slice_times, + slice_info=self.inputs.slice_info, + ) + + R.estimate(refscan=None) + + corr_run = R.resample() + self._out_file_path = [] + self._par_file_path = [] + + for j, corr in enumerate(corr_run): + self._out_file_path.append(os.path.abspath('corr_%s.nii.gz' % + (split_filename(self.inputs.in_file[j])[1]))) + save_image(corr, self._out_file_path[j]) + + self._par_file_path.append(os.path.abspath('%s.par' % + (os.path.split(self.inputs.in_file[j])[1]))) + mfile = open(self._par_file_path[j], 'w') + motion = R._transforms[j] + # nipy does not encode euler angles. return in original form of + # translation followed by rotation vector see: + # http://en.wikipedia.org/wiki/Rodrigues'_rotation_formula + for i, mo in enumerate(motion): + params = ['%.10f' % item for item in np.hstack((mo.translation, + mo.rotation))] + string = ' '.join(params) + '\n' + mfile.write(string) + mfile.close() + + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + outputs['out_file'] = self._out_file_path + outputs['par_file'] = self._par_file_path + return outputs + + class TrimInputSpec(BaseInterfaceInputSpec): in_file = File( exists=True, mandatory=True, diff --git a/nipype/pipeline/engine.py b/nipype/pipeline/engine.py index 81610c9404..54c0569ba9 100644 --- a/nipype/pipeline/engine.py +++ b/nipype/pipeline/engine.py @@ -1211,6 +1211,7 @@ def hash_exists(self, updatehash=False): outdir = self.output_dir() hashfiles = glob(os.path.join(outdir, '_0x*.json')) if len(hashfiles) > 1: + logger.info(hashfiles) logger.info('Removing multiple hashfiles and forcing node to rerun') for hashfile in hashfiles: os.unlink(hashfile) diff --git a/nipype/pipeline/plugins/multiproc.py b/nipype/pipeline/plugins/multiproc.py index d9718e3ccc..fe24878750 100644 --- a/nipype/pipeline/plugins/multiproc.py +++ b/nipype/pipeline/plugins/multiproc.py @@ -75,8 +75,14 @@ def _get_result(self, taskid): def _submit_job(self, node, updatehash=False): self._taskid += 1 + try: + if node.inputs.terminal_output == 'stream': + node.inputs.terminal_output = 'file' + except: + pass self._taskresult[self._taskid] = self.pool.apply_async(run_node, - (node, updatehash,)) + (node, + updatehash,)) return self._taskid def _report_crash(self, node, result=None):