unwarp_resting_state.py

from nipype.pipeline import engine as pe
from nipype.interfaces import fsl
from nipype.interfaces import io as nio
import os

from interfaces import MP2RageSkullStrip, MiconvResize

import pandas


DELTA_TE = 0.00102
SMOOTH_3D = 2.5
GRAPPA = 3


#ECHO_SPACING = 0.001 / GRAPPA #highres
ECHO_SPACING = 0.0013 / GRAPPA # lowres

workflow = pe.Workflow(name='unwarp_resting_state')

workflow.base_dir = os.path.abspath('/home/gdholla1/workflow_folders/')

dg = pe.Node(nio.DataGrabber(infields=['subject_id'], 
                             outfields=['fieldmap_scanner', 'magnitude_fieldmap', 'epi_resting_state', 'whole_brain_mp2rage',
                                        'whole_brain_t1', 'whole_brain_t1_weighted', 'whole_brain_inv2'], 
                             sort_filelist=True), name='datagrabber')

dg.inputs.base_directory = '/home/gdholla1/data/leipzig_2013/mri'
dg.inputs.template = '*'
dg.inputs.field_template = dict(fieldmap_scanner='%s/*_whole_brain/S*_fieldmap_DeltaTE1p02ms.nii',
                                magnitude_fieldmap='%s/*_whole_brain/S*_fieldmap_DeltaTE1p02ms_Te6.nii',
                                epi_resting_state='%s/*_whole_brain/S*_RS_*.nii',
								whole_brain_mp2rage='%s/*_zoomed_in/S2_mp2rage_whole_brain_INV1.nii',
								whole_brain_t1='%s/*_zoomed_in/S*_mp2rage_whole_brain_T1_Images.nii',
								whole_brain_t1_weighted='%s/*_zoomed_in/S*_mp2rage_whole_brain_UNI_Images.nii',
								whole_brain_inv2='%s/*_zoomed_in/S*_mp2rage_whole_brain_INV2.nii',)								
dg.inputs.template_args = dict(fieldmap_scanner=[['subject_id']],
                               magnitude_fieldmap=[['subject_id']],
                               epi_resting_state=[['subject_id']],
							   whole_brain_mp2rage=[['subject_id']],
							   	whole_brain_t1=[['subject_id']],
								whole_brain_t1_weighted=[['subject_id']],
								whole_brain_inv2=[['subject_id']],)

#dg.inputs.subject_id = 'SMVT'


subject_ids = list(pandas.read_csv('/home/gdholla1/data/leipzig_2013/subjects.tsv')['Leipzig-ID'])

subject_ids.pop(subject_ids.index('PF5T'))
subject_ids.pop(subject_ids.index('FMFT'))

dg.iterables = [('subject_id', subject_ids)]


ds = pe.Node(nio.DataSink(), name='datasink')
ds.inputs.base_directory = '/home/gdholla1/data/leipzig_2013/results'

ds.inputs.substitutions = [('_field_corrector_original0', 'original'),
                           ('_field_corrector_mc0', 'motion_corrected')]
                           
ds.inputs.regexp_substitutions = [('/_subject_id_(.*)/', '/\\1/')]


# returns the first three dimensions of nifti-file IN REVERSE ORDER
def get_dims(in_file):
    import nibabel as nb

    if type(in_file) == list:
        in_file = in_file[0]

    return nb.load(in_file).get_shape()[:3][::-1]


motion_correct = pe.MapNode(fsl.MCFLIRT(), name='motion_correct', iterfield=['in_file'])

# fslmaths $folder/$fieldmap_scanner -div 4096 -mul 3.14159 -div $DELTA_TE $folder/fieldmap.nii.gz -odt float
normalizer = pe.Node(fsl.maths.MathsCommand(), name='normalizer')
normalizer.inputs.args = '-div 4096 -mul 3.14159 -div %s' % DELTA_TE
normalizer.inputs.out_file = 'fieldmap.nii.gz'
normalizer.inputs.output_datatype = 'float'

# bet $folder/$Magnitude_fieldmap $folder/mag1_brain -m -f 0.5 -g 0
better = pe.Node(fsl.BET(), name='brain_extracter')
better.inputs.mask = True
better.inputs.frac = 0.5
better.inputs.vertical_gradient = 0.


resize_fieldmap = pe.Node(MiconvResize(), name='resize_fieldmap')
resize_brain_mask = pe.Node(MiconvResize(), name='resize_brain_mask')
resize_magnitude_fieldmap = pe.Node(MiconvResize(), name='resize_magnitude_fieldmap')

field_corrector_original = pe.MapNode(fsl.FUGUE(), name='field_corrector_original', iterfield=['in_file'])
field_corrector_mc = pe.MapNode(fsl.FUGUE(), name='field_corrector_mc', iterfield=['in_file'])

field_corrector_original.inputs.smooth3d = SMOOTH_3D
field_corrector_original.inputs.dwell_time = ECHO_SPACING
field_corrector_original.inputs.unwarp_direction = 'y-'

field_corrector_mc.inputs.smooth3d = SMOOTH_3D
field_corrector_mc.inputs.dwell_time = ECHO_SPACING
field_corrector_mc.inputs.unwarp_direction = 'y-'

workflow.connect(dg, 'epi_resting_state', motion_correct, 'in_file')
workflow.connect(dg, 'fieldmap_scanner', normalizer, 'in_file')
workflow.connect(dg, 'magnitude_fieldmap', better, 'in_file')

workflow.connect(normalizer, 'out_file', resize_fieldmap, 'in_file')
workflow.connect(better, 'mask_file', resize_brain_mask, 'in_file')
workflow.connect(dg, 'magnitude_fieldmap', resize_magnitude_fieldmap, 'in_file')

workflow.connect(dg, ('epi_resting_state', get_dims), resize_fieldmap, 'size')
workflow.connect(dg, ('epi_resting_state', get_dims), resize_brain_mask, 'size')
workflow.connect(dg, ('epi_resting_state', get_dims), resize_magnitude_fieldmap, 'size')

workflow.connect(dg, 'epi_resting_state', field_corrector_original, 'in_file')
workflow.connect(motion_correct, 'out_file', field_corrector_mc, 'in_file')

workflow.connect(resize_fieldmap, 'resized_file', field_corrector_original, 'fmap_in_file')
workflow.connect(resize_fieldmap, 'resized_file', field_corrector_mc, 'fmap_in_file')

workflow.connect(field_corrector_mc, 'unwarped_file', ds, 'whole_brain.resting_state.fieldmap_corrected.@mc')
workflow.connect(field_corrector_original, 'unwarped_file', ds, 'whole_brain.resting_state.fieldmap_corrected')

workflow.write_graph()
workflow.run(plugin='MultiProc', plugin_args={'n_procs' : 8})
#workflow.run()