!nib-ls /data/ds102/sub-01/*/*.nii.gz

"""
Explanation: Example 1: Preprocessing Workflow
This is meant as a very simple example for a preprocessing workflow. In this workflow we will conduct the following steps:

Motion Correction of functional images with FSL's MCFLIRT
Image Resampling of anatomical image to a resolution of 3mmx3mmx3mm voxel size with AFNI's 3dresample
Image Correction of motion corrected functional images to the resampled anatomical image with FSL's FLIRT
Smoothing of coregistrated functional images with FWHM set to 4mm and 8mm

For this example, we have 3 subjects with 2 different runs each. As a short recap, the image properties are:
End of explanation
"""


%pylab inline
from os.path import join as opj
from nipype.interfaces.fsl import MCFLIRT, FLIRT
from nipype.interfaces.afni import Resample
from nipype.interfaces.spm import Smooth
from nipype.interfaces.utility import IdentityInterface
from nipype.interfaces.io import SelectFiles, DataSink
from nipype.pipeline.engine import Workflow, Node

"""
Explanation: So, let's start!
Imports
First, let's import all modules we later will be needing.
End of explanation
"""


experiment_dir = '/output'
output_dir = 'datasink'
working_dir = 'workingdir'

# list of subject identifiers
subject_list = ['sub-01', 'sub-02', 'sub-03', 'sub-04', 'sub-05']

# list of session identifiers
session_list = ['run-1', 'run-2']

# Smoothing widths to apply
fwhm = [4, 8]

"""
Explanation: Experiment parameters
It's always a good idea to specify all parameters that might change between experiments at the beginning of your script.
End of explanation
"""


# MCFLIRT - motion correction
mcflirt = Node(MCFLIRT(mean_vol=True,
                       save_plots=True,
                       output_type='NIFTI'),
               name="mcflirt")

# Resample - resample anatomy to 3x3x3 voxel resolution
resample = Node(Resample(voxel_size=(3, 3, 3),
                         outputtype='NIFTI'),
                name="resample")

# FLIRT - coregister functional images to anatomical images
coreg_step1 = Node(FLIRT(output_type='NIFTI'), name="coreg_step1")
coreg_step2 = Node(FLIRT(output_type='NIFTI',
                         apply_xfm=True), name="coreg_step2")

# Smooth - image smoothing
smooth = Node(Smooth(), name="smooth")
smooth.iterables = ("fwhm", fwhm)

"""
Explanation: Specify Nodes
Initiate all the different interfaces (represented as nodes) that you want to use in your workflow.
End of explanation
"""


# Infosource - a function free node to iterate over the list of subject names
infosource = Node(IdentityInterface(fields=['subject_id', 'session_id']),
                  name="infosource")
infosource.iterables = [('subject_id', subject_list),
                        ('session_id', session_list)]

# SelectFiles - to grab the data (alternativ to DataGrabber)
anat_file = opj('{subject_id}', 'anat', '{subject_id}_T1w.nii.gz')
func_file = opj('{subject_id}', 'func',
                '{subject_id}_task-flanker_{session_id}_bold.nii.gz')

templates = {'anat': anat_file,
             'func': func_file}
selectfiles = Node(SelectFiles(templates,
                               base_directory='/data/ds102'),
                   name="selectfiles")

# Datasink - creates output folder for important outputs
datasink = Node(DataSink(base_directory=experiment_dir,
                         container=output_dir),
                name="datasink")

# Use the following DataSink output substitutions
substitutions = [('_subject_id', ''),
                 ('_session_id_', ''),
                 ('_task-flanker', ''),
                 ('_mcf.nii_mean_reg', '_mean'),
                 ('.nii.par', '.par'),
                 ]
subjFolders = [('%s_%s/' % (sess, sub), '%s/%s' % (sub, sess))
               for sess in session_list
               for sub in subject_list]
subjFolders += [('%s_%s' % (sub, sess), '')
                for sess in session_list
                for sub in subject_list]
subjFolders += [('%s%s_' % (sess, sub), '')
                for sess in session_list
                for sub in subject_list]
substitutions.extend(subjFolders)
datasink.inputs.substitutions = substitutions

"""
Explanation: Specify input & output stream
Specify where the input data can be found & where and how to save the output data.
End of explanation
"""


# Create a preprocessing workflow
preproc = Workflow(name='preproc')
preproc.base_dir = opj(experiment_dir, working_dir)

# Connect all components of the preprocessing workflow
preproc.connect([(infosource, selectfiles, [('subject_id', 'subject_id'),
                                            ('session_id', 'session_id')]),
                 (selectfiles, mcflirt, [('func', 'in_file')]),
                 (selectfiles, resample, [('anat', 'in_file')]),

                 (mcflirt, coreg_step1, [('mean_img', 'in_file')]),
                 (resample, coreg_step1, [('out_file', 'reference')]),

                 (mcflirt, coreg_step2, [('out_file', 'in_file')]),
                 (resample, coreg_step2, [('out_file', 'reference')]),
                 (coreg_step1, coreg_step2, [('out_matrix_file',
                                              'in_matrix_file')]),

                 (coreg_step2, smooth, [('out_file', 'in_files')]),

                 (mcflirt, datasink, [('par_file', 'preproc.@par')]),
                 (resample, datasink, [('out_file', 'preproc.@resample')]),
                 (coreg_step1, datasink, [('out_file', 'preproc.@coregmean')]),
                 (smooth, datasink, [('smoothed_files', 'preproc.@smooth')]),
                 ])

"""
Explanation: Specify Workflow
Create a workflow and connect the interface nodes and the I/O stream to each other.
End of explanation
"""


# Create preproc output graph
preproc.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=opj(preproc.base_dir, 'preproc', 'graph.dot.png'))

# Visualize the detailed graph
preproc.write_graph(graph2use='flat', format='png', simple_form=True)
Image(filename=opj(preproc.base_dir, 'preproc', 'graph_detailed.dot.png'))

"""
Explanation: Visualize the workflow
It always helps to visualize your workflow.
End of explanation
"""


preproc.run('MultiProc', plugin_args={'n_procs': 4})

"""
Explanation: Run the Workflow
Now that everything is ready, we can run the preprocessing workflow. Change n_procs to the number of jobs/cores you want to use.
End of explanation
"""


!tree /output/datasink/preproc/sub-01/

"""
Explanation: Inspect output
let's check the structure of the output folder, to see if we have everything we wanted to save.
End of explanation
"""


preproc.run('MultiProc', plugin_args={'n_procs': 4})
from nilearn import image, plotting
plotting.plot_epi(
    '/output/datasink/preproc/sub-01/run-1_bold_mean_flirt.nii', title="fwhm = 0mm",
    display_mode='ortho', annotate=False, draw_cross=False, cmap='gray')

mean_img = image.mean_img('/output/datasink/preproc/sub-01/run-1_fwhm_4/s_bold_mcf_flirt.nii')
plotting.plot_epi(mean_img, title="fwhm = 4mm", display_mode='ortho',
                  annotate=False, draw_cross=False, cmap='gray')

mean_img = image.mean_img('/output/datasink/preproc/sub-01/run-1_fwhm_8/s_bold_mcf_flirt.nii')
plotting.plot_epi(mean_img, title="fwhm = 8mm", display_mode='ortho',
                  annotate=False, draw_cross=False, cmap='gray')

"""
Explanation: Visualize results
let's check the effect of the different smoothing kernels.
End of explanation
"""


par = np.loadtxt('/output/datasink/preproc/sub-01/run-1_bold_mcf.par')
fig, axes = plt.subplots(2, 1, figsize=(15, 5))
axes[0].set_ylabel('rotation (radians)')
axes[0].plot(par[0:, :3])
axes[1].plot(par[0:, 3:])
axes[1].set_xlabel('time (TR)')
axes[1].set_ylabel('translation (mm)')

"""
Explanation: How do the motion parameters look like?
End of explanation
"""


