diff --git a/nianalysis/interfaces/custom/dicom.py b/nianalysis/interfaces/custom/dicom.py index 853f01ec..42c6e140 100644 --- a/nianalysis/interfaces/custom/dicom.py +++ b/nianalysis/interfaces/custom/dicom.py @@ -252,8 +252,12 @@ def _run_interface(self, runtime): pet_start_time = hd.AcquisitionTime except AttributeError: pet_start_time = None - with open(pet_image, 'r') as f: + with open(pet_image, 'rb') as f: for line in f: + try: + line = line[:-1].decode('utf-8') + except UnicodeDecodeError: + continue if 'image duration' in line: pet_duration = line.strip() pet_duration = int(pet_duration.split(':=')[-1]) diff --git a/nianalysis/interfaces/custom/fmri.py b/nianalysis/interfaces/custom/fmri.py index 14a1e7cf..c38deb8b 100644 --- a/nianalysis/interfaces/custom/fmri.py +++ b/nianalysis/interfaces/custom/fmri.py @@ -1,8 +1,12 @@ from nipype.interfaces.base import ( - BaseInterface, BaseInterfaceInputSpec, TraitedSpec, Directory, File) + BaseInterface, BaseInterfaceInputSpec, TraitedSpec, Directory, File, + traits) import os import shutil +import pydicom +import numpy as np +import glob class PrepareFIXInputSpec(BaseInterfaceInputSpec): @@ -59,7 +63,7 @@ def _run_interface(self, runtime): shutil.copy2(epi_mean, 'melodic_ica/mean_func.nii.gz') shutil.copytree(melodic_dir, 'melodic_ica/filtered_func_data.ica') shutil.copy2(filtered_epi, 'melodic_ica/filtered_func_data.nii.gz') - + with open('hand_label_file.txt', 'w') as f: f.write('not_provided') @@ -72,3 +76,43 @@ def _list_outputs(self): outputs["hand_label_file"] = os.getcwd()+'/hand_label_file.txt' return outputs + + +class FieldMapTimeInfoInputSpec(BaseInterfaceInputSpec): + + fm_mag = Directory() + + +class FieldMapTimeInfoOutputSpec(TraitedSpec): + + delta_te = traits.Float() + + +class FieldMapTimeInfo(BaseInterface): + + input_spec = FieldMapTimeInfoInputSpec + output_spec = FieldMapTimeInfoOutputSpec + + def _run_interface(self, runtime): + + fm_mag = sorted(glob.glob(self.inputs.fm_mag+'/*')) + tes = [pydicom.read_file(x).EchoTime for x in fm_mag] + tes = list(set(tes)) + if len(tes) != 2: + print('Something went wrong when trying to estimate ' + 'the delta TE between the two echos field map ' + 'images. delta_TE will be set equal to 2.46, ' + 'but please check if this is correct.') + self.delta_te = 2.46 + else: + self.delta_te = float(np.abs(tes[1]-tes[0])) + print('delta_TE: {}'.format(self.delta_te)) + + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + + outputs["delta_te"] = self.delta_te + + return outputs diff --git a/nianalysis/interfaces/custom/motion_correction.py b/nianalysis/interfaces/custom/motion_correction.py index ec47ca5d..3a4aee92 100644 --- a/nianalysis/interfaces/custom/motion_correction.py +++ b/nianalysis/interfaces/custom/motion_correction.py @@ -11,13 +11,15 @@ import scipy.ndimage.measurements as snm import datetime as dt try: - import matplotlib + import matplotlib + matplotlib.use('Agg') import matplotlib.pyplot as plot except ImportError: pass from nipype.interfaces import fsl import pydicom import math +import subprocess as sp class MotionMatCalculationInputSpec(BaseInterfaceInputSpec): @@ -501,9 +503,10 @@ class MeanDisplacementCalculation(BaseInterface): def _run_interface(self, runtime): - list_inputs = list(zip(self.inputs.motion_mats, self.inputs.start_times, - self.inputs.real_durations, self.inputs.trs, - self.inputs.input_names)) + list_inputs = list(zip( + self.inputs.motion_mats, self.inputs.start_times, + self.inputs.real_durations, self.inputs.trs, + self.inputs.input_names)) ref = nib.load(self.inputs.reference) ref_data = ref.get_data() # centre of gravity @@ -516,7 +519,7 @@ def _run_interface(self, runtime): .total_seconds(), x[2], x[3], x[4]) for x in list_inputs] study_len = int((list_inputs[-1][1]+float(list_inputs[-1][2]))*1000) mean_displacement_rc = np.zeros(study_len)-1 - motion_par_rc = np.zeros((6, study_len)) + motion_par_rc = np.zeros((6, study_len))-1 mean_displacement = [] motion_par = [] idt_mat = np.eye(4) @@ -604,6 +607,9 @@ def _run_interface(self, runtime): if (mean_displacement_rc[i] == -1 and mean_displacement_rc[i-1] != -1): mean_displacement_rc[i] = mean_displacement_rc[i-1] + for i in range(motion_par_rc.shape[1]): + if (motion_par_rc[0, i] == -1 and motion_par_rc[0, i-1] != -1): + motion_par_rc[:, i] = motion_par_rc[:, i-1] to_save = [mean_displacement, mean_displacement_consecutive, mean_displacement_rc, motion_par_rc, start_times, @@ -914,6 +920,8 @@ class PlotMeanDisplacementRCInputSpec(BaseInterfaceInputSpec): mean_disp_rc = File(exists=True, desc='Text file containing the mean ' 'displacement real clock.') + motion_par_rc = File(exists=True, desc='Text file containing the motion ' + 'parameters real clock.') frame_start_times = File(exists=True, desc='Frame start times as detected' 'by the motion framing pipeline') false_indexes = File(exists=True, desc='Time indexes were the scanner was ' @@ -925,6 +933,8 @@ class PlotMeanDisplacementRCInputSpec(BaseInterfaceInputSpec): class PlotMeanDisplacementRCOutputSpec(TraitedSpec): mean_disp_plot = File(exists=True, desc='Mean displacement plot.') + rot_plot = File(exists=True, desc='Rotation parameters plot.') + trans_plot = File(exists=True, desc='Translation parameters plot.') class PlotMeanDisplacementRC(BaseInterface): @@ -935,9 +945,13 @@ class PlotMeanDisplacementRC(BaseInterface): def _run_interface(self, runtime): mean_disp_rc = np.loadtxt(self.inputs.mean_disp_rc) - frame_start_times = np.loadtxt(self.inputs.frame_start_times) false_indexes = np.loadtxt(self.inputs.false_indexes, dtype=int) - framing = self.inputs.framing + + if isdefined(self.inputs.motion_par_rc): + motion_par_rc = np.loadtxt(self.inputs.motion_par_rc) + plot_mp = True + else: + plot_mp = False plot_offset = True dates = np.arange(0, len(mean_disp_rc), 1) indxs = np.zeros(len(mean_disp_rc), int)+1 @@ -953,28 +967,60 @@ def _run_interface(self, runtime): if len(start_true_period) == len(end_true_period)-1: end_true_period.remove(end_true_period[-1]) elif len(start_true_period) != len(end_true_period): - print ('Something went wrong in the indentification of the MR ' + print ('Something went wrong in the identification of the MR ' 'idling time. It will not be plotted.') plot_offset = False - fig, ax = plot.subplots() - fig.set_size_inches(21, 9) + self.gen_plot(dates, mean_disp_rc, plot_offset, start_true_period, + end_true_period) + if plot_mp: + for i in range(2): + mp = motion_par_rc[i*3:(i+1)*3, :] + self.gen_plot( + dates, mp, plot_offset, start_true_period, end_true_period, + plot_mp=plot_mp, mp_ind=i) + + return runtime + + def gen_plot(self, dates, to_plot, plot_offset, start_true_period, + end_true_period, plot_mp=False, mp_ind=None): + + frame_start_times = np.loadtxt(self.inputs.frame_start_times) + framing = self.inputs.framing font = {'weight': 'bold', 'size': 30} matplotlib.rc('font', **font) + fig, ax = plot.subplots() + fig.set_size_inches(21, 9) ax.set_xlim(0, dates[-1]) - ax.set_ylim(-0.3, np.max(mean_disp_rc) + 1) + if plot_mp: + col = ['b', 'g', 'r'] if plot_offset: for i in range(0, len(start_true_period)): - ax.plot(dates[start_true_period[i]-1:end_true_period[i]+1], - mean_disp_rc[start_true_period[i]-1: - end_true_period[i]+1], - c='b', linewidth=2) + if plot_mp: + for ii in range(3): + ax.plot( + dates[start_true_period[i]-1:end_true_period[i]+1], + to_plot[ii, start_true_period[i]-1: + end_true_period[i]+1], + c=col[ii], linewidth=2) + else: + ax.plot(dates[start_true_period[i]-1:end_true_period[i]+1], + to_plot[start_true_period[i]-1: + end_true_period[i]+1], c='b', linewidth=2) for i in range(0, len(end_true_period)-1): - ax.plot( - dates[end_true_period[i]-1:start_true_period[i+1]+1], - mean_disp_rc[end_true_period[i]-1: - start_true_period[i+1]+1], - c='b', linewidth=2, ls='--', dashes=(2, 3)) + if plot_mp: + for ii in range(3): + ax.plot( + dates[end_true_period[i]-1: + start_true_period[i+1]+1], + to_plot[ii, end_true_period[i]-1: + start_true_period[i+1]+1], + c=col[ii], linewidth=2, ls='--', dashes=(2, 3)) + else: + ax.plot( + dates[end_true_period[i]-1:start_true_period[i+1]+1], + to_plot[end_true_period[i]-1:start_true_period[i+1]+1], + c='b', linewidth=2, ls='--', dashes=(2, 3)) if framing: cl = 'yellow' @@ -1006,19 +1052,35 @@ def _run_interface(self, runtime): cl = 'yellow' indx = np.arange(0, len(dates), 300000) - my_thick = [str(i) for i in np.arange(0, len(dates)/60000, 5)] + my_thick = [str(i) for i in np.arange(0, len(dates)/60000, 5, + dtype=int)] plot.xticks(dates[indx], my_thick) - - plot.savefig('mean_displacement_real_clock.png') + plot.xlabel('Time [min]', fontsize=25) + if mp_ind == 0: + ax.set_ylim(np.min(to_plot)-0.1, np.max(to_plot)+0.1) + plot.legend(['Rotation X', 'Rotation Y', 'Rotation Z'], loc=0) + plot.ylabel('Rotation [rad]', fontsize=25) + plot.savefig('Rotation_real_clock.png') + elif mp_ind == 1: + ax.set_ylim(np.min(to_plot)-0.5, np.max(to_plot)+0.5) + plot.legend( + ['Translation X', 'Translation Y', 'Translation Z'], loc=0) + plot.ylabel('Translation [mm]', fontsize=25) + plot.savefig('Translation_real_clock.png') + else: + plot.ylabel('Mean displacement [mm]', fontsize=25) + plot.savefig('mean_displacement_real_clock.png') plot.close() - return runtime - def _list_outputs(self): outputs = self._outputs().get() outputs["mean_disp_plot"] = ( os.getcwd()+'/mean_displacement_real_clock.png') + outputs["rot_plot"] = ( + os.getcwd()+'/Rotation_real_clock.png') + outputs["trans_plot"] = ( + os.getcwd()+'/Translation_real_clock.png') return outputs @@ -1449,3 +1511,46 @@ def _list_outputs(self): outputs["average_bin_mats"] = os.getcwd()+'/average_bin_mats' return outputs + + +class ReorientUmapInputSpec(BaseInterfaceInputSpec): + + umap = Directory(exists=True) + niftis = traits.List() + + +class ReorientUmapOutputSpec(TraitedSpec): + + reoriented_umaps = traits.List() + + +class ReorientUmap(BaseInterface): + """Had to write this pipeline because I found out that sometimes the + reoriented umaps (in nifti format) had different orientation with respect + to the original umap (dicom). This pipeline uses mrconvert to extract the + strides from the dicom umap and to apply them to each of the nifti umaps. + """ + input_spec = ReorientUmapInputSpec + output_spec = ReorientUmapOutputSpec + + def _run_interface(self, runtime): + + niftis = self.inputs.niftis + umap = self.inputs.umap + + for nifti in niftis: + _, base, ext = split_filename(nifti) + outname = base+'_reorient'+ext + cmd = ('mrconvert -strides {0} {1} {2}' + .format(umap, nifti, outname)) + sp.check_output(cmd, shell=True) + + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + + outputs["reoriented_umaps"] = sorted(glob.glob( + os.getcwd()+'/*_reorient.*')) + + return outputs diff --git a/nianalysis/interfaces/custom/pet.py b/nianalysis/interfaces/custom/pet.py index 3d381713..f1022873 100644 --- a/nianalysis/interfaces/custom/pet.py +++ b/nianalysis/interfaces/custom/pet.py @@ -650,7 +650,7 @@ def _run_interface(self, runtime): def get_qform(self, image): cmd = 'mrinfo {}'.format(image) - hd = sp.check_output(cmd, shell=True) + hd = (sp.check_output(cmd, shell=True)).decode('utf-8') i = [n for n, el in enumerate(hd.split('\n')) if 'Transform' in el][0] mat = [] for j in range(3): @@ -755,7 +755,8 @@ def applyxfm(self, in_file, ref, mat, outname): def extract_qform(self, image): cmd = 'fslhd {}'.format(image) - image_info = sp.check_output(cmd, shell=True).strip().split('\n') + image_info = (sp.check_output(cmd, shell=True)).decode('utf-8') + image_info = image_info.strip().split('\n') qform = np.eye(4) for i, line in enumerate(image_info): if 'qto_xyz:1' in line: diff --git a/nianalysis/interfaces/mrtrix/preproc.py b/nianalysis/interfaces/mrtrix/preproc.py index 6504890c..442a293f 100644 --- a/nianalysis/interfaces/mrtrix/preproc.py +++ b/nianalysis/interfaces/mrtrix/preproc.py @@ -44,7 +44,7 @@ class DWIPreprocInputSpec(MRTrix3BaseInputSpec): "these two images must have the same phase"), position=3) eddy_parameters = traits.Str( - argstr='-eddy_parameters "%s"', desc='parameters to be passed to eddy') + argstr='-eddy_options "%s"', desc='parameters to be passed to eddy') no_clean_up = traits.Bool(True, argstr='-nocleanup', desc='Do not delete the temporary folder') temp_dir = Directory(genfile=True, argstr='-tempdir %s', @@ -220,6 +220,7 @@ class DWIDenoiseInputSpec(MRTrix3BaseInputSpec): "*nix) are found in the provided output file then the CWD (when " "the workflow is run, i.e. the working directory) will be " "prepended to the output path.")) + out_file_ext = traits.Str(desc='Extention of the output file.') noise = File( genfile=True, argstr="-noise %s", desc=("The estimated spatially-varying noise level")) @@ -261,8 +262,13 @@ def _gen_outfname(self): if isdefined(self.inputs.out_file): out_name = self.inputs.out_file else: - base, ext = split_extension( - os.path.basename(self.inputs.in_file)) + if isdefined(self.inputs.out_file_ext): + ext = self.inputs.out_file_ext + base, _ = split_extension(os.path.basename( + self.inputs.in_file)) + else: + base, ext = split_extension( + os.path.basename(self.inputs.in_file)) out_name = os.path.join(os.getcwd(), "{}_conv{}".format(base, ext)) return out_name @@ -271,8 +277,13 @@ def _gen_noisefname(self): if isdefined(self.inputs.out_file): out_name = self.inputs.out_file else: - base, ext = split_extension( - os.path.basename(self.inputs.in_file)) + if isdefined(self.inputs.out_file_ext): + ext = self.inputs.out_file_ext + base, _ = split_extension(os.path.basename( + self.inputs.in_file)) + else: + base, ext = split_extension( + os.path.basename(self.inputs.in_file)) out_name = os.path.join(os.getcwd(), "{}_noise{}".format(base, ext)) return out_name diff --git a/nianalysis/interfaces/mrtrix/utils.py b/nianalysis/interfaces/mrtrix/utils.py index bd388efa..db4db820 100644 --- a/nianalysis/interfaces/mrtrix/utils.py +++ b/nianalysis/interfaces/mrtrix/utils.py @@ -375,6 +375,7 @@ class MRCalcInputSpec(CommandLineInputSpec): out_file = File(genfile=True, argstr='%s', position=-1, desc="Extracted DW or b-zero images") + out_ext = traits.Str(desc='Extention of the output file.') operation = traits.Enum( 'abs', 'neg', 'sqrt', 'exp', 'log', 'log10', 'cos', 'sin', 'tan', @@ -418,8 +419,11 @@ def _gen_outfilename(self): if isdefined(self.inputs.out_file): filename = self.inputs.out_file else: - _, ext = split_extension( - os.path.basename(self.inputs.operands[0])) + if isdefined(self.inputs.out_ext): + ext = self.inputs.out_ext + else: + _, ext = split_extension( + os.path.basename(self.inputs.operands[0])) filename = os.getcwd() for op in self.inputs.operands: try: diff --git a/nianalysis/motion_correction_utils.py b/nianalysis/motion_correction_utils.py index f591ce57..323f0933 100644 --- a/nianalysis/motion_correction_utils.py +++ b/nianalysis/motion_correction_utils.py @@ -341,10 +341,12 @@ def dwi_type_assignment(input_dir, dmri_images): for dwi in dmri_images: cmd = 'mrinfo {0}'.format(input_dir+'/'+dwi) - info = sp.check_output(cmd, shell=True).strip().split('\n') + info = (sp.check_output(cmd, shell=True)).decode('utf-8') + info = info.strip().split('\n') for line in info: if 'Dimensions:' in line: dim = line.split('Dimensions:')[-1].strip().split('x') + break hd_extraction = DicomHeaderInfoExtraction() hd_extraction.inputs.dicom_folder = input_dir+'/'+dwi dcm_info = hd_extraction.run() @@ -376,9 +378,9 @@ def dwi_type_assignment(input_dir, dmri_images): for j in range(len(b0)): ped_b0 = b0[j][2] if ped_b0 == ped_main: - if main_dwi[i][1] == b0[j][1]: + if main_dwi[i][1] == b0[j][1] and (j == i or j == i+1): dmris.append([b0[j][0], '1']) - elif main_dwi[i][1] != b0[j][1]: + elif main_dwi[i][1] != b0[j][1] and (j == i or j == i+1): dmris.append([b0[j][0], '-1']) else: unused_b0.append(b0[j][0]) @@ -424,10 +426,12 @@ def check_image_start_time(input_dir, scans): toremove = [] for scan in scans: try: +# scan_name = scan.split('-')[1] + scan_number = scan.split('-')[0].zfill(3) hd_extraction = DicomHeaderInfoExtraction() hd_extraction.inputs.dicom_folder = input_dir+'/'+scan dcm_info = hd_extraction.run() - start_times.append([dcm_info.outputs.start_time, scan]) + start_times.append([dcm_info.outputs.start_time, scan_number, scan]) except: print(('This folder {} seems to not contain DICOM files. It will ' 'be ingnored.'.format(scan))) @@ -436,8 +440,8 @@ def check_image_start_time(input_dir, scans): diff = ((dt.datetime.strptime(start_times[i][0], '%H%M%S.%f') - dt.datetime.strptime(start_times[i-1][0], '%H%M%S.%f')) .total_seconds()) - if diff < 10: - toremove.append(start_times[i][1]) + if diff < 5: + toremove.append(start_times[i][-1]) return toremove diff --git a/nianalysis/study/mri/base.py b/nianalysis/study/mri/base.py index ad028db3..aaa4fcf7 100644 --- a/nianalysis/study/mri/base.py +++ b/nianalysis/study/mri/base.py @@ -22,7 +22,7 @@ import logging from nianalysis.interfaces.ants import AntsRegSyn from nipype.interfaces.ants.resampling import ApplyTransforms -from arcana.parameter import ParameterSpec +from arcana.parameter import ParameterSpec, SwitchSpec from nianalysis.interfaces.custom.motion_correction import ( MotionMatCalculation) @@ -84,8 +84,6 @@ class MRIStudy(Study, metaclass=StudyMetaClass): ParameterSpec('bet_f_threshold', 0.5), ParameterSpec('bet_reduce_bias', False), ParameterSpec('bet_g_threshold', 0.0), - ParameterSpec('bet_method', 'fsl_bet', - choices=('fsl_bet', 'optibet')), ParameterSpec('MNI_template', os.path.join(atlas_path, 'MNI152_T1_2mm.nii.gz')), ParameterSpec('MNI_template_brain', @@ -93,16 +91,12 @@ class MRIStudy(Study, metaclass=StudyMetaClass): ParameterSpec('MNI_template_mask', os.path.join( atlas_path, 'MNI152_T1_2mm_brain_mask.nii.gz')), ParameterSpec('optibet_gen_report', False), - ParameterSpec('atlas_coreg_tool', 'ants', - choices=('fnirt', 'ants')), ParameterSpec('fnirt_atlas', 'MNI152'), ParameterSpec('fnirt_resolution', '2mm'), ParameterSpec('fnirt_intensity_model', 'global_non_linear_with_bias'), ParameterSpec('fnirt_subsampling', [4, 4, 2, 2, 1, 1]), ParameterSpec('preproc_new_dims', ('RL', 'AP', 'IS')), ParameterSpec('preproc_resolution', None, dtype=list), - ParameterSpec('linear_reg_method', 'flirt', - choices=('flirt', 'spm', 'ants')), ParameterSpec('flirt_degrees_of_freedom', 6, desc=( "Number of degrees of freedom used in the registration. " "Default is 6 -> affine transformation.")), @@ -114,6 +108,14 @@ class MRIStudy(Study, metaclass=StudyMetaClass): "Whether to use the QS form supplied in the input image " "header (the image coordinates of the FOV supplied by the " "scanner"))] + + add_switch_specs = [ + SwitchSpec('linear_reg_method', 'flirt', + choices=('flirt', 'spm', 'ants')), + SwitchSpec('atlas_coreg_tool', 'ants', + choices=('fnirt', 'ants')), + SwitchSpec('bet_method', 'fsl_bet', + choices=('fsl_bet', 'optibet'))] @property def coreg_brain_spec(self): @@ -328,7 +330,7 @@ def _optiBET_brain_extraction_pipeline(self, in_file, **kwargs): outputs = [DatasetSpec('brain', nifti_gz_format), DatasetSpec('brain_mask', nifti_gz_format)] - if self.switch('optibet_gen_report'): + if self.parameter('optibet_gen_report'): outputs.append(DatasetSpec('optiBET_report', gif_format)) pipeline = self.create_pipeline( name='brain_extraction', diff --git a/nianalysis/study/mri/epi.py b/nianalysis/study/mri/epi.py index 618db1e6..b3b4c962 100644 --- a/nianalysis/study/mri/epi.py +++ b/nianalysis/study/mri/epi.py @@ -5,18 +5,19 @@ from arcana.dataset import DatasetSpec, FieldSpec from nianalysis.file_format import ( nifti_gz_format, text_matrix_format, directory_format, - par_format, motion_mats_format) + par_format, motion_mats_format, dicom_format) from nianalysis.citation import fsl_cite from nipype.interfaces import fsl from nianalysis.requirement import fsl509_req from arcana.study.base import StudyMetaClass from nianalysis.interfaces.custom.motion_correction import ( MergeListMotionMat, MotionMatCalculation) -from arcana.parameter import ParameterSpec +from arcana.parameter import ParameterSpec, SwitchSpec from nipype.interfaces.utility import Merge as merge_lists from nipype.interfaces.fsl.utils import Merge as fsl_merge from nipype.interfaces.fsl.epi import PrepareFieldmap from nipype.interfaces.fsl.preprocess import BET, FUGUE +from nianalysis.interfaces.custom.fmri import FieldMapTimeInfo class EPIStudy(MRIStudy, metaclass=StudyMetaClass): @@ -32,13 +33,19 @@ class EPIStudy(MRIStudy, metaclass=StudyMetaClass): DatasetSpec('align_mats', directory_format, 'intrascan_alignment_pipeline'), DatasetSpec('moco_par', par_format, - 'intrascan_alignment_pipeline')] + 'intrascan_alignment_pipeline'), + FieldSpec('field_map_delta_te', float, + 'field_map_time_info_pipeline')] add_parameter_specs = [ ParameterSpec('bet_robust', True), ParameterSpec('bet_f_threshold', 0.2), ParameterSpec('bet_reduce_bias', False), - ParameterSpec('linear_reg_method', 'epireg')] + ParameterSpec('fugue_echo_spacing', 0.000275)] + + add_switch_specs = [ + SwitchSpec('linear_reg_method', 'epireg', + choices=('flirt', 'spm', 'ants', 'epireg'))] def linear_coregistration_pipeline(self, **kwargs): if self.branch('linear_reg_method', 'epireg'): @@ -104,6 +111,25 @@ def intrascan_alignment_pipeline(self, **kwargs): return pipeline + def field_map_time_info_pipeline(self, **kwargs): + + pipeline = self.create_pipeline( + name='field_map_time_info_pipeline', + inputs=[DatasetSpec('field_map_mag', dicom_format)], + outputs=[FieldSpec('field_map_delta_te', float)], + desc=("Pipeline to extract delta TE from field map " + "images, if provided"), + version=1, + citations=[fsl_cite], + **kwargs) + + delta_te = pipeline.create_node(FieldMapTimeInfo(), + name='extract_delta_te') + pipeline.connect_input('field_map_mag', delta_te, 'fm_mag') + pipeline.connect_output('field_map_delta_te', delta_te, 'delta_te') + + return pipeline + def preproc_pipeline(self, **kwargs): if ('field_map_phase' in self.input_names and @@ -183,7 +209,8 @@ def _fugue_pipeline(self, **kwargs): name='preproc_pipeline', inputs=[DatasetSpec('primary', nifti_gz_format), DatasetSpec('field_map_mag', nifti_gz_format), - DatasetSpec('field_map_phase', nifti_gz_format)], + DatasetSpec('field_map_phase', nifti_gz_format), + FieldSpec('field_map_delta_te', float)], outputs=[DatasetSpec('preproc', nifti_gz_format)], desc=("Fugue distortion correction pipeline"), version=1, @@ -210,7 +237,9 @@ def _fugue_pipeline(self, **kwargs): create_fmap = pipeline.create_node( PrepareFieldmap(), name="prepfmap", wall_time=5, requirements=[fsl509_req]) - create_fmap.inputs.delta_TE = 2.46 +# create_fmap.inputs.delta_TE = 2.46 + pipeline.connect_input('field_map_delta_te', create_fmap, + 'delta_TE') pipeline.connect(bet, "out_file", create_fmap, "in_magnitude") pipeline.connect(fm_phase_reorient, 'out_file', create_fmap, 'in_phase') @@ -218,7 +247,7 @@ def _fugue_pipeline(self, **kwargs): fugue = pipeline.create_node(FUGUE(), name='fugue', wall_time=5, requirements=[fsl509_req]) fugue.inputs.unwarp_direction = 'x' - fugue.inputs.dwell_time = 0.000275 + fugue.inputs.dwell_time = self.parameter('fugue_echo_spacing') fugue.inputs.unwarped_file = 'example_func.nii.gz' pipeline.connect(create_fmap, 'out_fieldmap', fugue, 'fmap_in_file') diff --git a/nianalysis/study/mri/functional/fmri.py b/nianalysis/study/mri/functional/fmri.py index d755a2c2..b1ae8523 100644 --- a/nianalysis/study/mri/functional/fmri.py +++ b/nianalysis/study/mri/functional/fmri.py @@ -15,7 +15,7 @@ from nipype.interfaces.utility import Merge as NiPypeMerge import os from nipype.interfaces.utility.base import IdentityInterface -from arcana.parameter import ParameterSpec +from arcana.parameter import ParameterSpec, SwitchSpec from nianalysis.study.mri.epi import EPIStudy from nipype.interfaces.ants.resampling import ApplyTransforms from nianalysis.study.mri.structural.t1 import T1Study @@ -27,7 +27,6 @@ from nianalysis.interfaces.c3d import ANTs2FSLMatrixConversion import logging from arcana.exception import ArcanaNameError -from arcana.interfaces.utils import CopyToDir logger = logging.getLogger('nianalysis') @@ -42,16 +41,15 @@ class FunctionalMRIStudy(EPIStudy, metaclass=StudyMetaClass): - add_option_specs = [ + add_parameter_specs = [ ParameterSpec('component_threshold', 20), ParameterSpec('motion_reg', True), ParameterSpec('highpass', 0.01), ParameterSpec('brain_thresh_percent', 5), ParameterSpec('MNI_template', - os.path.join(atlas_path, 'MNI152_T1_2mm.nii.gz')), + os.path.join(atlas_path, 'MNI152_T1_2mm.nii.gz')), ParameterSpec('MNI_template_mask', os.path.join( atlas_path, 'MNI152_T1_2mm_brain_mask.nii.gz')), - ParameterSpec('linear_reg_method', 'ants'), ParameterSpec('group_ica_components', 15)] add_data_specs = [ @@ -74,6 +72,10 @@ class FunctionalMRIStudy(EPIStudy, metaclass=StudyMetaClass): DatasetSpec('smoothed_ts', nifti_gz_format, 'smoothing_pipeline')] + add_switch_specs = [ + SwitchSpec('linear_reg_method', 'ants', + choices=('flirt', 'spm', 'ants', 'epireg'))] + def rsfMRI_filtering_pipeline(self, **kwargs): pipeline = self.create_pipeline( @@ -226,77 +228,6 @@ def fix_preparation_pipeline(self, **kwargs): return pipeline -# def fix_training_pipeline(self, **kwargs): -# -# inputs = [] -# sub_study_names = [] -# for sub_study_spec in self.sub_study_specs(): -# try: -# spec = self.data_spec(sub_study_spec.inverse_map('fix_dir')) -# spec._format = directory_format -# inputs.append(spec) -# inputs.append( -# self.data_spec(sub_study_spec.inverse_map( -# 'hand_label_noise'))) -# sub_study_names.append(sub_study_spec.name) -# except ArcanaNameError: -# continue # Sub study doesn't have fix dir -# -# pipeline = self.create_pipeline( -# name='training_fix', -# inputs=inputs, -# outputs=[DatasetSpec('train_data', rfile_format)], -# desc=("Pipeline to create the training set for FIX given a group " -# "of subjects with the hand_label_noise.txt file within " -# "their fix_dir."), -# version=1, -# citations=[fsl_cite], -# **kwargs) -# -# num_fix_dirs = len(sub_study_names) -# merge_fix_dirs = pipeline.create_node(NiPypeMerge(num_fix_dirs), -# name='merge_fix_dirs') -# merge_label_files = pipeline.create_node(NiPypeMerge(num_fix_dirs), -# name='merge_label_files') -# for i, sub_study_name in enumerate(sub_study_names, start=1): -# spec = self.sub_study_spec(sub_study_name) -# pipeline.connect_input( -# spec.inverse_map('fix_dir'), merge_fix_dirs, 'in{}'.format(i)) -# pipeline.connect_input( -# spec.inverse_map('hand_label_noise'), merge_label_files, -# 'in{}'.format(i)) -# -# merge_visits = pipeline.create_join_visits_node( -# IdentityInterface(['list_dir', 'list_label_files']), -# joinfield=['list_dir', 'list_label_files'], name='merge_visits') -# merge_subjects = pipeline.create_join_subjects_node( -# NiPypeMerge(2), joinfield=['in1', 'in2'], name='merge_subjects') -# merge_subjects.inputs.ravel_inputs = True -# -# prepare_training = pipeline.create_node(PrepareFIXTraining(), -# name='prepare_training') -# prepare_training.inputs.epi_number = num_fix_dirs -# pipeline.connect(merge_fix_dirs, 'out', merge_visits, 'list_dir') -# pipeline.connect(merge_visits, 'list_dir', merge_subjects, 'in1') -# pipeline.connect(merge_label_files, 'out', merge_visits, -# 'list_label_files') -# pipeline.connect(merge_visits, 'list_label_files', merge_subjects, -# 'in2') -# pipeline.connect(merge_subjects, 'out', prepare_training, -# 'inputs_list') -# -# fix_training = pipeline.create_node( -# FSLFixTraining(), name='fix_training', -# wall_time=240, requirements=[fix_req]) -# fix_training.inputs.outname = 'FIX_training_set' -# fix_training.inputs.training = True -# pipeline.connect(prepare_training, 'prepared_dirs', fix_training, -# 'list_dir') -# -# pipeline.connect_output('train_data', fix_training, 'training_set') -# -# return pipeline - def fix_classification_pipeline(self, **kwargs): pipeline = self.create_pipeline( @@ -413,9 +344,8 @@ class FunctionalMRIMixin(MultiStudy, metaclass=MultiStudyMetaClass): add_data_specs = [ DatasetSpec('train_data', rfile_format, 'fix_training_pipeline', frequency='per_project'), -# DatasetSpec('fmri_pre-processeing_results', directory_format, -# 'gather_fmri_result_pipeline'), - DatasetSpec('group_melodic', directory_format, 'group_melodic_pipeline')] + DatasetSpec('group_melodic', directory_format, + 'group_melodic_pipeline')] def fix_training_pipeline(self, **kwargs): @@ -487,40 +417,6 @@ def fix_training_pipeline(self, **kwargs): pipeline.connect_output('train_data', fix_training, 'training_set') return pipeline - -# def gather_fmri_result_pipeline(self, **kwargs): -# -# inputs = [] -# sub_study_names = [] -# for sub_study_spec in self.sub_study_specs(): -# try: -# inputs.append( -# self.data_spec(sub_study_spec.inverse_map('smoothed_ts'))) -# sub_study_names.append(sub_study_spec.name) -# except ArcanaNameError: -# continue # Sub study doesn't have fix dir -# -# pipeline = self.create_pipeline( -# name='gather_fmri', -# inputs=inputs, -# outputs=[DatasetSpec('fmri_pre-processeing_results', directory_format)], -# desc=("Pipeline to gather together all the pre-processed fMRI images"), -# version=1, -# citations=[fsl_cite], -# **kwargs) -# -# merge_inputs = pipeline.create_node(NiPypeMerge(len(inputs)), -# name='merge_inputs') -# for i, sub_study_name in enumerate(sub_study_names, start=1): -# spec = self.sub_study_spec(sub_study_name) -# pipeline.connect_input( -# spec.inverse_map('smoothed_ts'), merge_inputs, 'in{}'.format(i)) -# -# copy2dir = pipeline.create_node(CopyToDir(), name='copy2dir') -# pipeline.connect(merge_inputs, 'out', copy2dir, 'in_files') -# -# pipeline.connect_output('fmri_pre-processeing_results', copy2dir, 'out_dir') -# return pipeline def group_melodic_pipeline(self, **kwargs): @@ -554,8 +450,8 @@ def group_melodic_pipeline(self, **kwargs): return pipeline -def create_fmri_study_class(name, t1, epis, epi_number, fm_mag=None, - fm_phase=None, run_regression=False): +def create_fmri_study_class(name, t1, epis, epi_number, echo_spacing, + fm_mag=None, fm_phase=None, run_regression=False): inputs = [] dct = {} @@ -586,18 +482,22 @@ def create_fmri_study_class(name, t1, epis, epi_number, fm_mag=None, epi_refspec.update({'t1_wm_seg': 'coreg_ref_wmseg', 't1_preproc': 'coreg_ref_preproc', 'train_data': 'train_data'}) - - study_specs.extend(SubStudySpec('epi_{}'.format(i), FunctionalMRIStudy, - epi_refspec) - for i in range(epi_number)) + study_specs.append(SubStudySpec('epi_0', FunctionalMRIStudy, epi_refspec)) + if epi_number > 1: + epi_refspec.update({'t1_wm_seg': 'coreg_ref_wmseg', + 't1_preproc': 'coreg_ref_preproc', + 'train_data': 'train_data', + 'epi_0_coreg_to_atlas_warp': 'coreg_to_atlas_warp', + 'epi_0_coreg_to_atlas_mat': 'coreg_to_atlas_mat'}) + study_specs.extend(SubStudySpec('epi_{}'.format(i), FunctionalMRIStudy, + epi_refspec) + for i in range(1, epi_number)) for i in range(epi_number): inputs.append(DatasetMatch('epi_{}_primary'.format(i), dicom_format, epis, order=i, is_regex=True)) -# inputs.extend(DatasetMatch( -# 'epi_{}_hand_label_noise'.format(i), text_format, -# 'hand_label_noise_{}'.format(i+1)) -# for i in range(epi_number)) + parameter_specs.append( + ParameterSpec('epi_{}_fugue_echo_spacing'.format(i), echo_spacing)) if distortion_correction: inputs.extend(DatasetMatch( @@ -621,4 +521,5 @@ def create_fmri_study_class(name, t1, epis, epi_number, fm_mag=None, dct['add_data_specs'] = data_specs dct['add_parameter_specs'] = parameter_specs dct['__metaclass__'] = MultiStudyMetaClass - return MultiStudyMetaClass(name, (FunctionalMRIMixin,), dct), inputs, output_files + return (MultiStudyMetaClass(name, (FunctionalMRIMixin,), dct), inputs, + output_files) diff --git a/nianalysis/study/mri/structural/diffusion.py b/nianalysis/study/mri/structural/diffusion.py index 40f11292..5dd15bd3 100644 --- a/nianalysis/study/mri/structural/diffusion.py +++ b/nianalysis/study/mri/structural/diffusion.py @@ -85,6 +85,7 @@ class DiffusionStudy(EPIStudy, metaclass=StudyMetaClass): SwitchSpec('preproc_denoise', False), SwitchSpec('response_algorithm', 'tax', ('tax', 'dhollander', 'msmt_5tt')), + SwitchSpec('fod_algorithm', 'csd', ('csd', 'msmt_csd')), SwitchSpec('brain_extract_method', 'mrtrix', ('mrtrix', 'fsl')), SwitchSpec('bias_correct_method', 'ants', @@ -142,11 +143,13 @@ def preproc_pipeline(self, **kwargs): # @UnusedVariable @IgnorePep8 # Run denoising denoise = pipeline.create_node(DWIDenoise(), name='denoise', requirements=[mrtrix3_req]) + denoise.inputs.out_file_ext = '.mif' # Calculate residual noise subtract_operands = pipeline.create_node(Merge(2), name='subtract_operands') subtract = pipeline.create_node(MRCalc(), name='subtract', requirements=[mrtrix3_req]) + subtract.inputs.out_ext = '.mif' subtract.inputs.operation = 'subtract' dwipreproc = pipeline.create_node( DWIPreproc(), name='dwipreproc', diff --git a/nianalysis/study/mri/structural/t1.py b/nianalysis/study/mri/structural/t1.py index 3d730946..fea1794b 100644 --- a/nianalysis/study/mri/structural/t1.py +++ b/nianalysis/study/mri/structural/t1.py @@ -9,7 +9,7 @@ from arcana.interfaces.utils import JoinPath from ..base import MRIStudy from arcana.study.base import StudyMetaClass -from arcana.parameter import ParameterSpec +from arcana.parameter import ParameterSpec, SwitchSpec class T1Study(MRIStudy, metaclass=StudyMetaClass): @@ -20,11 +20,13 @@ class T1Study(MRIStudy, metaclass=StudyMetaClass): DatasetSpec('brain', nifti_gz_format, 'brain_extraction_pipeline')] add_parameter_specs = [ - ParameterSpec('bet_method', 'optibet', - choices=MRIStudy.parameter_spec('bet_method').choices), ParameterSpec('bet_robust', True), ParameterSpec('bet_f_threshold', 0.57), ParameterSpec('bet_g_threshold', -0.1)] + + add_switch_specs = [ + SwitchSpec('bet_method', 'optibet', + choices=MRIStudy.switch_spec('bet_method').choices)] def freesurfer_pipeline(self, **kwargs): """ diff --git a/nianalysis/study/multimodal/mrpet.py b/nianalysis/study/multimodal/mrpet.py index a28d910d..84df11de 100644 --- a/nianalysis/study/multimodal/mrpet.py +++ b/nianalysis/study/multimodal/mrpet.py @@ -5,7 +5,7 @@ from nianalysis.interfaces.custom.motion_correction import ( MeanDisplacementCalculation, MotionFraming, PlotMeanDisplacementRC, AffineMatAveraging, PetCorrectionFactor, CreateMocoSeries, FixedBinning, - UmapAlign2Reference) + UmapAlign2Reference, ReorientUmap) from nianalysis.citation import fsl_cite from arcana.study.multi import ( MultiStudy, SubStudySpec, MultiStudyMetaClass) @@ -22,7 +22,7 @@ from nianalysis.interfaces.custom.pet import ( CheckPetMCInputs, PetImageMotionCorrection, StaticPETImageGeneration, PETFovCropping) -from arcana.option import OptionSpec +from arcana.parameter import ParameterSpec, SwitchSpec import os from nianalysis.interfaces.converters import Nii2Dicom from arcana.interfaces.utils import CopyToDir, ListDir, dicom_fname_sort_key @@ -43,12 +43,12 @@ logging.getLogger("urllib3").setLevel(logging.WARNING) reference_path = os.path.abspath( - os.path.join(os.path.dirname(__file__), '../', + os.path.join(os.path.dirname(__file__), '../../', 'reference_data')) template_path = os.path.abspath( - os.path.join(os.path.dirname(__file__), '../', - 'templates')) + os.path.join(os.path.dirname(__file__).split('nianalysis')[0], + 'nianalysis', 'nianalysis', 'templates')) class MotionDetectionMixin(MultiStudy, metaclass=MultiStudyMetaClass): @@ -99,6 +99,10 @@ class MotionDetectionMixin(MultiStudy, metaclass=MultiStudyMetaClass): 'motion_framing_pipeline'), DatasetSpec('mean_displacement_plot', png_format, 'plot_mean_displacement_pipeline'), + DatasetSpec('rotation_plot', png_format, + 'plot_mean_displacement_pipeline'), + DatasetSpec('translation_plot', png_format, + 'plot_mean_displacement_pipeline'), DatasetSpec('average_mats', directory_format, 'frame_mean_transformation_mats_pipeline'), DatasetSpec('correction_factors', text_format, @@ -120,27 +124,27 @@ class MotionDetectionMixin(MultiStudy, metaclass=MultiStudyMetaClass): FieldSpec('pet_start_time', dtype=str, pipeline_name='pet_header_info_extraction_pipeline')] - add_option_specs = [OptionSpec('framing_th', 2.0), - OptionSpec('framing_temporal_th', 30.0), - OptionSpec('framing_duration', 0), - OptionSpec('md_framing', True), - OptionSpec('align_pct', False), - OptionSpec('align_fixed_binning', False), - OptionSpec('moco_template', os.path.join( + add_parameter_specs = [ParameterSpec('framing_th', 2.0), + ParameterSpec('framing_temporal_th', 30.0), + ParameterSpec('framing_duration', 0), + ParameterSpec('md_framing', True), + ParameterSpec('align_pct', False), + ParameterSpec('align_fixed_binning', False), + ParameterSpec('moco_template', os.path.join( reference_path, 'moco_template.IMA')), - OptionSpec('PET_template_MNI', os.path.join( + ParameterSpec('PET_template_MNI', os.path.join( template_path, 'PET_template_MNI.nii.gz')), - OptionSpec('fixed_binning_n_frames', 0), - OptionSpec('pet_offset', 0), - OptionSpec('fixed_binning_bin_len', 60), - OptionSpec('crop_xmin', 100), - OptionSpec('crop_xsize', 130), - OptionSpec('crop_ymin', 100), - OptionSpec('crop_ysize', 130), - OptionSpec('crop_zmin', 20), - OptionSpec('crop_zsize', 100), - OptionSpec('PET2MNI_reg', False), - OptionSpec('dynamic_pet_mc', False)] + ParameterSpec('fixed_binning_n_frames', 0), + ParameterSpec('pet_offset', 0), + ParameterSpec('fixed_binning_bin_len', 60), + ParameterSpec('crop_xmin', 100), + ParameterSpec('crop_xsize', 130), + ParameterSpec('crop_ymin', 100), + ParameterSpec('crop_ysize', 130), + ParameterSpec('crop_zmin', 20), + ParameterSpec('crop_zsize', 100), + ParameterSpec('PET2MNI_reg', False), + ParameterSpec('dynamic_pet_mc', False)] def mean_displacement_pipeline(self, **kwargs): inputs = [DatasetSpec('ref_brain', nifti_gz_format)] @@ -276,9 +280,12 @@ def plot_mean_displacement_pipeline(self, **kwargs): pipeline = self.create_pipeline( name='plot_mean_displacement', inputs=[DatasetSpec('mean_displacement_rc', text_format), + DatasetSpec('motion_par_rc', text_format), DatasetSpec('offset_indexes', text_format), DatasetSpec('frame_start_times', text_format)], - outputs=[DatasetSpec('mean_displacement_plot', png_format)], + outputs=[DatasetSpec('mean_displacement_plot', png_format), + DatasetSpec('rotation_plot', png_format), + DatasetSpec('translation_plot', png_format)], desc=("Plot the mean displacement real clock"), version=1, citations=[fsl_cite], @@ -293,8 +300,14 @@ def plot_mean_displacement_pipeline(self, **kwargs): 'false_indexes') pipeline.connect_input('frame_start_times', plot_md, 'frame_start_times') + pipeline.connect_input('motion_par_rc', plot_md, + 'motion_par_rc') pipeline.connect_output('mean_displacement_plot', plot_md, 'mean_disp_plot') + pipeline.connect_output('rotation_plot', plot_md, + 'rot_plot') + pipeline.connect_output('translation_plot', plot_md, + 'trans_plot') return pipeline def frame_mean_transformation_mats_pipeline(self, **kwargs): @@ -386,6 +399,9 @@ def nifti2dcm_conversion_pipeline(self, **kwargs): **kwargs) list_niftis = pipeline.create_node(ListDir(), name='list_niftis') + reorient_niftis = pipeline.create_node( + ReorientUmap(), name='reorient_niftis', requirements=[mrtrix3_req]) + nii2dicom = pipeline.create_map_node( Nii2Dicom(), name='nii2dicom', iterfield=['in_file'], wall_time=20) @@ -395,12 +411,15 @@ def nifti2dcm_conversion_pipeline(self, **kwargs): copy2dir = pipeline.create_node(CopyToDir(), name='copy2dir') copy2dir.inputs.extension = 'Frame' # Connect nodes - pipeline.connect(list_niftis, 'files', nii2dicom, 'in_file') + pipeline.connect(list_niftis, 'files', reorient_niftis, 'niftis') + pipeline.connect(reorient_niftis, 'reoriented_umaps', nii2dicom, + 'in_file') pipeline.connect(list_dicoms, 'files', nii2dicom, 'reference_dicom') pipeline.connect(nii2dicom, 'out_file', copy2dir, 'in_files') # Connect inputs pipeline.connect_input('umaps_align2ref', list_niftis, 'directory') pipeline.connect_input('umap', list_dicoms, 'directory') + pipeline.connect_input('umap', reorient_niftis, 'umap') # Connect outputs pipeline.connect_output('umap_aligned_dicoms', copy2dir, 'out_dir') @@ -416,19 +435,17 @@ def umap_realignment_pipeline(self, **kwargs): inputs.append(DatasetSpec('umap', nifti_gz_format)) outputs.append(DatasetSpec('umaps_align2ref', directory_format)) pipeline = self.create_pipeline( - name='static_frame2ref_alignment', + name='umap_realignment', inputs=inputs, outputs=outputs, - desc=("Pipeline to create an affine mat to align each " - "detected frame to the reference. If umap is provided" - ", it will be also aligned to match the head position" - " in each frame and improve the static PET image " - "quality."), + desc=("Pipeline to align the original umap (if provided)" + "to match the head position in each frame and improve the " + "static PET image quality."), version=1, citations=[fsl_cite], **kwargs) frame_align = pipeline.create_node( - UmapAlign2Reference(), name='static_frame2ref_alignment', + UmapAlign2Reference(), name='umap2ref_alignment', requirements=[fsl509_req]) frame_align.inputs.pct = self.parameter('align_pct') pipeline.connect_input('umap_ref_coreg_matrix', frame_align, @@ -443,7 +460,10 @@ def umap_realignment_pipeline(self, **kwargs): return pipeline def create_moco_series_pipeline(self, **kwargs): - + """This pipeline is probably wrong as we still do not know how to + import back the new moco series into the scanner. This was just a first + attempt. + """ pipeline = self.create_pipeline( name='create_moco_series', inputs=[DatasetSpec('start_times', text_format), @@ -509,7 +529,7 @@ def motion_correction_pipeline(self, **kwargs): inputs = [DatasetSpec('pet_data_prepared', directory_format), DatasetSpec('ref_brain', nifti_gz_format), DatasetSpec('mean_displacement_plot', png_format)] - if self.option_spec('dynamic_pet_mc').value: + if self.parameter_spec('dynamic_pet_mc').value: inputs.append(DatasetSpec('fixed_binning_mats', directory_format)) outputs = [DatasetSpec('dynamic_motion_correction_results', directory_format)] @@ -713,7 +733,8 @@ def create_motion_correction_class(name, ref=None, ref_type=None, t1s=None, dct = {} data_specs = [] run_pipeline = False - option_specs = [OptionSpec('ref_preproc_resolution', [1])] + parameter_specs = [ParameterSpec('ref_preproc_resolution', [1])] + switch_specs = [] if struct2align is not None: struct_image = struct2align.split('/')[-1].split('.')[0] @@ -728,7 +749,7 @@ def create_motion_correction_class(name, ref=None, ref_type=None, t1s=None, DatasetMatch('struct2align', nifti_gz_format, struct_image)) if pet_data_dir is not None and pet_recon_dir is not None and dynamic: output_data = 'dynamic_motion_correction_results' - option_specs.append(OptionSpec('dynamic_pet_mc', True)) + parameter_specs.append(ParameterSpec('dynamic_pet_mc', True)) if struct2align is not None: inputs.append( DatasetMatch('struct2align', nifti_gz_format, struct_image)) @@ -826,9 +847,9 @@ def create_motion_correction_class(name, ref=None, ref_type=None, t1s=None, dwi_refspec.update({'ref_wm_seg': 'coreg_ref_wmseg', 'ref_preproc': 'coreg_ref_preproc'}) if dmris_main: - option_specs.extend( - OptionSpec('dwi_{}_brain_extract_method'.format(i), 'fsl') - for i in range(len(dmris_main))) + switch_specs.extend( + SwitchSpec('dwi_{}_brain_extract_method'.format(i), 'fsl', + ('mrtrix', 'fsl')) for i in range(len(dmris_main))) if dmris_main and not dmris_opposite: logger.warning( 'No opposite phase encoding direction b0 provided. DWI ' @@ -912,7 +933,8 @@ def create_motion_correction_class(name, ref=None, ref_type=None, t1s=None, dct['add_sub_study_specs'] = study_specs dct['add_data_specs'] = data_specs dct['__metaclass__'] = MultiStudyMetaClass - dct['add_option_specs'] = option_specs + dct['add_parameter_specs'] = parameter_specs + dct['add_switch_specs'] = switch_specs return (MultiStudyMetaClass(name, (MotionDetectionMixin,), dct), inputs, output_data) @@ -925,7 +947,7 @@ def create_motion_detection_class(name, ref=None, ref_type=None, t1s=None, dct = {} data_specs = [] run_pipeline = False - option_specs = [OptionSpec('ref_preproc_resolution', [1])] + parameter_specs = [ParameterSpec('ref_preproc_resolution', [1])] if pet_data_dir is not None: inputs.append(DatasetMatch('pet_data_dir', directory_format, @@ -1065,5 +1087,5 @@ def create_motion_detection_class(name, ref=None, ref_type=None, t1s=None, dct['add_sub_study_specs'] = study_specs dct['add_data_specs'] = data_specs dct['__metaclass__'] = MultiStudyMetaClass - dct['add_option_specs'] = option_specs + dct['add_parameter_specs'] = parameter_specs return MultiStudyMetaClass(name, (MotionDetectionMixin,), dct), inputs diff --git a/scripts/run_diffusion_preprocessing.py b/scripts/run_diffusion_preprocessing.py new file mode 100755 index 00000000..ce548f0a --- /dev/null +++ b/scripts/run_diffusion_preprocessing.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +import os.path +import errno +from nianalysis.study.mri.structural.diffusion import DiffusionStudy +from arcana.repository.xnat import XnatRepository +from nianalysis.file_format import dicom_format +import logging +import argparse +from arcana.dataset.match import DatasetMatch +from arcana.runner.linear import LinearRunner + + +if __name__ == "__main__": + + parser = argparse.ArgumentParser() + parser.add_argument('--subject', type=str, nargs='+', default=None, + help="Subject IDs to process") + parser.add_argument('--session', type=str, nargs='+', default=None, + help="Session IDs to process") + parser.add_argument('--study_name', type=str, default='diffusion', + help="Study name to be prepend to the output names " + "of all pre-processing results. Default is " + "'diffusion'.") + args = parser.parse_args() + + scratch_dir = os.path.expanduser('~/scratch') + + WORK_PATH = os.path.join(scratch_dir, 'mnd', 'diffusion') + try: + os.makedirs(WORK_PATH) + except OSError as e: + if e.errno != errno.EEXIST: + raise + + logger = logging.getLogger('NiAnalysis') + logger.setLevel(logging.DEBUG) + # Stream Handler + handler = logging.StreamHandler() + formatter = logging.Formatter("%(levelname)s - %(message)s") + handler.setFormatter(formatter) + logger.addHandler(handler) + # File Handler + handler = logging.FileHandler(os.path.join(WORK_PATH, 'out.log')) + formatter = logging.Formatter("%(levelname)s - %(message)s") + handler.setFormatter(formatter) + logger.addHandler(handler) + + inputs = [ + DatasetMatch('primary', dicom_format, + 'R-L_MRtrix_60_directions_interleaved_B0_ep2d_diff_p2'), + DatasetMatch('dwi_reference', dicom_format, + 'L-R_MRtrix_60_directions_interleaved_B0_ep2d_diff_p2')] + + study = DiffusionStudy( + name=args.study_name, + repository=XnatRepository( + project_id='MRH060', server='https://mbi-xnat.erc.monash.edu.au', + cache_dir=os.path.join(scratch_dir, 'xnat_cache-mnd')), + runner=LinearRunner(work_dir=os.path.join(scratch_dir, + 'xnat_working_dir-mnd')), + inputs=inputs, subject_ids=args.subject, visit_ids=args.session, + parameters={'preproc_pe_dir': 'RL'}, + switches={'preproc_denoise': True}) + + fods = study.data('fod') + # print(fods[0].path) + print('Done') diff --git a/scripts/run_fmri_preprocessing.py b/scripts/run_fmri_preprocessing.py old mode 100644 new mode 100755 index 5b47b45a..74119fc2 --- a/scripts/run_fmri_preprocessing.py +++ b/scripts/run_fmri_preprocessing.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python3 from nianalysis.study.mri.functional.fmri import create_fmri_study_class from arcana.repository.xnat import XnatRepository from arcana.runner.linear import LinearRunner @@ -18,15 +19,24 @@ 'and to carry out all the processing.') parser.add_argument('--hires_structural', '-struct', type=str, required=True, help='High resolution structural image ' - 'to used to improve the registation between fMRI and ' - 'MNI template (usually is a T1 weigthed image).') + 'to used to improve the registration between fMRI and ' + 'MNI template (usually is a T1 weighted image).') parser.add_argument('--fmri', type=str, required=True, help='Regular expression to match the name of the fMRI' ' images on XNAT to be pre-processed. If more than one' ' image match this expression, please provide the ' '--fmri_order as well. Also, in case of multiple fmri ' - ' images, the same field map images (if provided) will' + 'images, the same field map images (if provided) will' ' be used to perform B0 unwarping.') + parser.add_argument('--fmri_echo_spacing', type=float, required=True, + help='Echo spacing (in seconds) from the EPI ' + 'acquisition parameter. Please be aware that if you ' + 'acquired your EPI with acceleration factors ' + '(i.e. iPAT) you have to divide the echo spacing by ' + 'that factor. For example, if the iPAT was 2 than you ' + 'have to provide echo_spacing/2. ' + 'N.B. If you have multiband factor you DO NOT have to ' + 'divide the echo spacing by that factor!') parser.add_argument('--fmri_order', type=int, required=True, help='If more than one fmri image is going to match ' 'the --fmri regular expression provided, you can ' @@ -50,12 +60,12 @@ default='https://mbi-xnat.erc.monash.edu.au') parser.add_argument('--xnat_username', '-user', type=str, help='Username with which to connect to XNAT with. ' - 'This can be skiped if it has already been saved in ' + 'This can be skipped if it has already been saved in ' 'the .netrc in your home directory, otherwise it is ' 'mandatory.', default=None) parser.add_argument('--xnat_password', '-password', type=str, help='Password to connect to XNAt with. ' - 'This can be skiped if it has already been saved in ' + 'This can be skipped if it has already been saved in ' 'the .netrc in your home directory, otherwise it is ' 'mandatory.', default=None) parser.add_argument('--field_map_mag', '-mag', type=str, @@ -71,19 +81,20 @@ 'phase image is the output of a SIEMENS scanner. It ' 'does not support other vendors.', default=None) - parser.add_argument('--run_regression', '-regression', type=str, + parser.add_argument('--run_regression', '-regression', action='store_true', help='If ' 'provided, fix classification and regression of the ' 'noisy component will be performed and the final image' ' will be fully pre-processed. Otherwise, the ' 'pipeline will generate only the MELODIC L1 results ' 'with the right folder structure so that it can be ' - 'used with fsl FIX.', default=None) + 'used with fsl FIX.', default=False) args = parser.parse_args() fMRI, inputs, output_files = create_fmri_study_class( 'fMRI', args.hires_structural, args.fmri, args.fmri_order, - fm_mag=args.field_map_mag, fm_phase=args.field_map_phase, + args.fmri_echo_spacing, fm_mag=args.field_map_mag, + fm_phase=args.field_map_phase, run_regression=args.run_regression) CACHE_PATH = os.path.join(args.working_dir, 'xnat_cache') diff --git a/scripts/run_motion_correction.py b/scripts/run_motion_correction.py index 21c4c487..ae7ae9b4 100755 --- a/scripts/run_motion_correction.py +++ b/scripts/run_motion_correction.py @@ -55,10 +55,22 @@ def create_motion_correction_inputs(self): pr) = pkl.load(f) working_dir = ( input_dir+'/work_dir/work_sub_dir/work_session_dir/') - if pet_dir is not None and pd != pet_dir: + if pet_dir is not None and not pd and pd != pet_dir: shutil.copytree(pet_dir, working_dir+'/pet_data_dir') if pet_recon is not None and pr != pet_recon: - shutil.copytree(pet_dir, working_dir+'/pet_data_dir') + if pr: + print('Different PET recon dir, respect to that ' + 'provided in a previous run. The directory ' + 'pet_data_reconstructed in the working directory' + ' will be removed and substituted with the new ' + 'one.') + shutil.rmtree(working_dir+'/pet_data_reconstructed') + shutil.copytree(pet_recon, working_dir + + '/pet_data_reconstructed') + list_inputs = [ref, ref_type, t1s, epis, t2s, dmris, pd, + pet_recon] + with open(cache_input_path, 'wb') as f: + pkl.dump(list_inputs, f) cached_inputs = True except IOError as e: if e.errno == errno.ENOENT: