def anat_preproc(file_to_register, register_to, warp_back, pipeline_dir):
    # DATA CONFIGURATION. FOLLOWING OPENFMRI STANDARD.
    save_to = os.path.join(pipeline_dir,
                           file_to_register.split('/')[-1].split('.')[0])
    # Run pipeline imperatively with caching (without workflow object)
    mem = Memory(pipeline_dir)
    antsreg = mem.cache(Registration)
    transform = mem.cache(ApplyTransforms)
    save_list = []
    # nodes manual parameter configuration and run
    reg = antsreg(args='--float',
                  collapse_output_transforms=True,
                  moving_image=file_to_register,
                  fixed_image=register_to,
                  initial_moving_transform_com=True,
                  num_threads=n_proc,
                  output_inverse_warped_image=True,
                  output_warped_image=True,
                  sigma_units=['vox']*3,
                  transforms=['Rigid', 'Affine', 'SyN'],
                  terminal_output='file',
                  winsorize_lower_quantile=0.005,
                  winsorize_upper_quantile=0.995,
                  convergence_threshold=[1e-06],
                  convergence_window_size=[10],
                  metric=['MI', 'MI', 'CC'],
                  metric_weight=[1.0]*3,
                  number_of_iterations=[[1000, 500, 250, 100],
                                        [1000, 500, 250, 100],
                                        [100, 70, 50, 20]],
                  radius_or_number_of_bins=[32, 32, 4],
                  sampling_percentage=[0.25, 0.25, 1],
                  sampling_strategy=['Regular',
                                     'Regular',
                                     'None'],
                  shrink_factors=[[8, 4, 2, 1]]*3,
                  smoothing_sigmas=[[3, 2, 1, 0]]*3,
                  transform_parameters=[(0.1,),
                                        (0.1,),
                                        (0.1, 3.0, 0.0)],
                  use_histogram_matching=True,
                  write_composite_transform=True)
    save_list.append([reg.outputs.composite_transform, save_to])
    save_list.append([reg.outputs.warped_image, save_to])
    save_list.append([reg.outputs.inverse_composite_transform, save_to])
    save_list.append([reg.outputs.inverse_warped_image, save_to])
    transformed = transform(args='--float',
                            input_image_type=3,
                            interpolation='NearestNeighbor',
                            invert_transform_flags=[False],
                            num_threads=n_proc,
                            reference_image=file_to_register,
                            terminal_output='file',
                            transforms=reg.outputs.inverse_composite_transform,
                            input_image=warp_back)
    save_list.append([transformed.outputs.output_image, save_to])
    return save_list

def convert_rawdata(base_directory, input_dir, out_prefix):
    os.environ['UNPACK_MGH_DTI'] = '0'
    file_list = os.listdir(input_dir)

    # If RAWDATA folder contains one (and only one) gunzipped nifti file -> copy it
    first_file = os.path.join(input_dir, file_list[0])
    if len(file_list) == 1 and first_file.endswith('nii.gz'):
        copyfile(first_file, os.path.join(base_directory, 'NIFTI', out_prefix+'.nii.gz'), False, False, 'content') # intelligent copy looking at input's content
    else:
        mem = Memory(base_dir=os.path.join(base_directory,'NIPYPE'))
        mri_convert = mem.cache(fs.MRIConvert)
        res = mri_convert(in_file=first_file, out_file=os.path.join(base_directory, 'NIFTI', out_prefix + '.nii.gz'))
        if len(res.outputs.get()) == 0:
            return False

    return True

def test_caching():
    temp_dir = mkdtemp(prefix='test_memory_')
    old_rerun = config.get('execution', 'stop_on_first_rerun')
    try:
        # Prevent rerun to check that evaluation is computed only once
        config.set('execution', 'stop_on_first_rerun', 'true')
        mem = Memory(temp_dir)
        first_nb_run = nb_runs
        results = mem.cache(SideEffectInterface)(input1=2, input2=1)
        assert_equal(nb_runs, first_nb_run + 1)
        assert_equal(results.outputs.output1, [1, 2])
        results = mem.cache(SideEffectInterface)(input1=2, input2=1)
        # Check that the node hasn't been rerun
        assert_equal(nb_runs, first_nb_run + 1)
        assert_equal(results.outputs.output1, [1, 2])
        results = mem.cache(SideEffectInterface)(input1=1, input2=1)
        # Check that the node hasn been rerun
        assert_equal(nb_runs, first_nb_run + 2)
        assert_equal(results.outputs.output1, [1, 1])
    finally:
        rmtree(temp_dir)
        config.set('execution', 'stop_on_first_rerun', old_rerun)

from procasl import preprocessing, _utils
current_directory = os.getcwd()
for (func_file, anat_file) in zip(
        heroes['func ASL'], heroes['anat']):
    # Create a memory context
    subject_directory = os.path.relpath(anat_file, subjects_parent_directory)
    subject_directory = subject_directory.split(os.sep)[0]
    cache_directory = os.path.join(os.path.expanduser('~/CODE/process-asl'),
                                   'procasl_cache', 'heroes',
                                   subject_directory)
    if not os.path.exists(cache_directory):
        os.mkdir(cache_directory)

    # nipype saves .m scripts into cwd
    os.chdir(cache_directory)
    mem = Memory(cache_directory)

    # Get Tag/Control sequence
    get_tag_ctl = mem.cache(preprocessing.RemoveFirstScanControl)
    out_get_tag_ctl = get_tag_ctl(in_file=func_file)

    # Rescale
    rescale = mem.cache(preprocessing.Rescale)
    out_rescale = rescale(in_file=out_get_tag_ctl.outputs.tag_ctl_file,
                          ss_tr=35.4, t_i_1=800., t_i_2=1800.)

    # Realign to first scan
    realign = mem.cache(preprocessing.ControlTagRealign)
    out_realign = realign(
        in_file=out_rescale.outputs.rescaled_file,
        register_to_mean=False,

current_directory = os.getcwd()

# Loop over subjects
for (func_file, anat_file) in zip(
        heroes['BOLD EPI'], heroes['anat']):
    # Create a memory context
    subject_directory = os.path.relpath(anat_file, subjects_parent_directory)
    subject_directory = subject_directory.split(os.sep)[0]
    cache_directory = os.path.join(os.path.expanduser('~/CODE/process-asl'),
                                   'procasl_cache', 'heroes',
                                   subject_directory)
    if not os.path.exists(cache_directory):
        os.mkdir(cache_directory)

    os.chdir(cache_directory)  # nipype saves .m scripts in current directory
    mem = Memory(cache_directory)

    # Realign EPIs
    realign = mem.cache(spm.Realign)
    out_realign = realign(
        in_files=func_file,
        register_to_mean=True)

    # Coregister anat to mean EPIs
    coregister = mem.cache(spm.Coregister)
    out_coregister = coregister(
        target=out_realign.outputs.mean_image,
        source=anat_file,
        write_interp=3,
        jobtype='estimate')

            print nifti_file, anat_image
            shutil.move(nifti_file, anat_image)
        else:
            print nifti_file, fmri_sessions[session_id]
            shutil.move(nifti_file, fmri_sessions[session_id])

        # remove the dicom dirs
        for x in glob.glob(os.path.join(dicom_dir, '*')):
            os.remove(x)
        os.removedirs(dicom_dir)

    ##############################################################
    # Preprocessing
    ##############################################################

    mem = Memory(base_dir=subject_dir)

    ##############################################################
    # Anatomical segmentation (White/Grey matter)

    seg = mem.cache(spm.Segment)
    
    out_seg = seg(data=anat_image,
                  gm_output_type=[True, True, True],
                  wm_output_type=[True, True, True],
                  csf_output_type=[True, True, True])
    sn_file = out_seg.outputs.transformation_mat
    inv_sn_file = out_seg.outputs.inverse_transformation_mat
    gm_image = out_seg.outputs.normalized_gm_image
    native_gm_image = out_seg.outputs.native_gm_image

import numpy as np
from cfutils import get_subjects, get_subject_data

X = get_subjects()
_, pdata = get_subject_data(X)
X = pdata.subject
y = pdata.lsas_pre - pdata.lsas_post

lgroup,_ = get_subject_data(X[y<=np.median(y)])
hgroup,_ = get_subject_data(X[y>np.median(y)])

import nipype.interfaces.spm as spm

from nipype.caching import Memory
os.makedirs('/mindhive/scratch/satra/sadfigures/nipype_mem')
mem = Memory('/mindhive/scratch/satra/sadfigures')

designer = mem.cache(spm.OneSampleTTestDesign)
estimator = mem.cache(spm.EstimateModel)
cestimator = mem.cache(spm.EstimateContrast)

ldesres =  designer(in_files = lgroup)
lestres = estimator(spm_mat_file=ldesres.outputs.spm_mat_file,
                    estimation_method={'Classical':None})
lcestres = cestimator(spm_mat_file=lestres.outputs.spm_mat_file,
                      beta_images=lestres.outputs.beta_images,
                      residual_image=lestres.outputs.residual_image,
                      group_contrast=True,
                      contrasts=[('LGroup', 'T', ['mean'], [1])])

hdesres =  designer(in_files = hgroup)

# Loop over subjects
for (func_file, anat_file) in zip(
        heroes['basal ASL'], heroes['anat']):
    # Create a memory context
    subject_directory = os.path.relpath(anat_file, subjects_parent_directory)
    subject_directory = subject_directory.split(os.sep)[0]
    cache_directory = os.path.join(os.path.expanduser('~/CODE/process-asl'),
                                   'procasl_cache', 'heroes',
                                   subject_directory)
    if not os.path.exists(cache_directory):
        os.mkdir(cache_directory)

    # nipype saves .m scripts into cwd
    os.chdir(cache_directory)
    mem = Memory(cache_directory)

    # Get Tag/Control sequence
    get_tag_ctl = mem.cache(preprocessing.RemoveFirstScanControl)
    out_get_tag_ctl = get_tag_ctl(in_file=func_file)

    # Rescale
    rescale = mem.cache(preprocessing.Rescale)
    out_rescale = rescale(in_file=out_get_tag_ctl.outputs.tag_ctl_file,
                          ss_tr=35.4, t_i_1=800., t_i_2=1800.)

    # Realign to first scan
    realign = mem.cache(preprocessing.Realign)
    out_realign = realign(
        in_file=out_rescale.outputs.rescaled_file,
        register_to_mean=False,

    out.runtime.cwd
"""

from nipype.interfaces import fsl
fsl.FSLCommand.set_default_output_type('NIFTI')

from nipype.caching import Memory

import glob

# First retrieve the list of files that we want to work upon
in_files = glob.glob('data/*/f3.nii')

# Create a memory context
mem = Memory('.')

# Apply an arbitrary (and pointless, here) threshold to the files)
threshold = [mem.cache(fsl.Threshold)(in_file=f, thresh=i)
                        for i, f in enumerate(in_files)]

# Merge all these files along the time dimension
out_merge = mem.cache(fsl.Merge)(dimension="t",
                            in_files=[t.outputs.out_file for t in threshold],
                        )
# And finally compute the mean
out_mean = mem.cache(fsl.MeanImage)(in_file=out_merge.outputs.merged_file)

# To avoid having increasing disk size we can keep only what was touched
# in this run
#mem.clear_previous_runs()

# Load functional ASL image of HEROES dataset first subject
import os
from procasl import datasets

heroes = datasets.load_heroes_dataset(
    subjects=(0,),
    subjects_parent_directory=os.path.join(os.path.expanduser("~/procasl_data"), "heroes"),
    paths_patterns={"raw ASL": "fMRI/acquisition1/vismot1_rawASL*.nii"},
)
raw_asl_file = heroes["raw ASL"][0]

# Create a memory context
from nipype.caching import Memory

cache_directory = "/tmp"
mem = Memory("/tmp")
os.chdir(cache_directory)
# Rescale
from procasl import preprocessing

rescale = mem.cache(preprocessing.Rescale)
out_rescale = rescale(in_file=raw_asl_file, ss_tr=35.4, t_i_1=800.0, t_i_2=1800.0)

# Plot the first volume before and after rescaling
from nilearn import plotting
import matplotlib.pylab as plt

for filename, title in zip([raw_asl_file, out_rescale.outputs.rescaled_file], ["raw", "rescaled"]):
    figure = plt.figure(figsize=(5, 4))
    first_scan_file = preprocessing.save_first_scan(filename)
    plotting.plot_img(first_scan_file, figure=figure, display_mode="z", cut_coords=(65,), title=title, colorbar=True)

def from_native_to_mni(img, sub_id, include_trans=[True, True, True],
                       interpolation='Linear'):
    '''Maps image from native space to mni.

    WARNING THERE IS A CLEAR PROBLEM IN THE UNDERSTANDING OF TRANSFORM ORDER
    WHEN ONLY USING THE LAST TWO TRANSFORMS THE ORDER SHOULD BE INVERTED

    We assume that the transformation files already exist for the mappings
    between:
    1) mean bold and anatomy
    2) anatomy and oasis template
    3) oasis template and mni template

    The transforms to include are:
    1) From bold to anat
    2) From anat to oasis
    3) From oasis to mni

    The include transforms should be sequential to have meaninful output,
    which means that transformations sequence [True, False, True] is invalid.
    '''
    check = (include_trans == [True, False, True])
    if check:
        raise Exception('Invalid transformation sequence')
    pipeline_dir = 'pipelines/transformations'
    if not os.path.exists(pipeline_dir):
        os.makedirs(pipeline_dir)
    mem = Memory(pipeline_dir)
    transform = mem.cache(ApplyTransforms)

    anat = os.path.join('pipelines',
                        'preprocessing',
                        'sub{0}'.format(sub_id),
                        'highres001.nii')
    oasis_template = os.path.join('pipelines',
                                  'OASIS-30_Atropos_template',
                                  'T_template0.nii.gz')
    mni_template = os.path.join('pipelines',
                                'mni_icbm152_nlin_asym_09a_nifti',
                                'mni_icbm152_nlin_asym_09a',
                                'mni_icbm152_t1_tal_nlin_asym_09a.nii')
    bold_to_anat = os.path.join('pipelines', 'preprocessing',
                                'sub{0}'.format(sub_id),
                                'bold_to_anat.txt')
    anat_to_oasis = os.path.join('pipelines', 'preprocessing',
                                 'sub{0}'.format(sub_id),
                                 'anat_to_oasis.h5')
    oasis_to_mni = os.path.join('pipelines', 'preprocessing',
                                'registered_templates', 'oasis_to_mni.h5')
    all_references = [anat, oasis_template, mni_template]
    all_trans = [bold_to_anat, anat_to_oasis, oasis_to_mni]
    all_inv_trans = [False, False, False]
    transforms = []
    inv_trans_flags = []
    reference = None
    for idx, flag in enumerate(include_trans):
        if flag:
            transforms.append(all_trans[idx])
            inv_trans_flags.append(all_inv_trans[idx])
            # Use latest transformation as reference
            reference = all_references[idx]

    trans = transform(args='--float',
                      input_image_type=3,
                      interpolation=interpolation,
                      invert_transform_flags=inv_trans_flags[::-1],
                      num_threads=n_proc,
                      reference_image=reference,
                      terminal_output='file',
                      transforms=transforms[::-1],
                      input_image=img)

    return trans.outputs.output_image

Realignment demo
================

This example compares standard realignement to realignement with tagging
correction.
"""
# Load 4D ASL image of KIRBY dataset first subject
import os
from procasl import datasets
kirby = datasets.fetch_kirby(subjects=[4])
raw_asl_file = kirby.asl[0]

# Create a memory context
from nipype.caching import Memory
cache_directory = '/tmp'
mem = Memory('/tmp')
os.chdir(cache_directory)

# Realign with and without tagging correction
from procasl import preprocessing
import numpy as np
realign = mem.cache(preprocessing.ControlTagRealign)
x_translation = {}
for correct_tagging in [True, False]:
    out_realign = realign(in_file=raw_asl_file,
                          correct_tagging=correct_tagging)
    x_translation[correct_tagging] = np.loadtxt(
        out_realign.outputs.realignment_parameters)[:, 2]

# Plot x-translation parameters with and without tagging correction
import matplotlib.pylab as plt

    def check_input(self, gui=True):
        print "**** Check Inputs ****"
        diffusion_available = False
        t1_available = False
        t2_available = False
        valid_inputs = False

        mem = Memory(base_dir=os.path.join(self.base_directory, "NIPYPE"))
        swap_and_reorient = mem.cache(SwapAndReorient)

        # Check for (and if existing, convert) diffusion data
        diffusion_model = []
        for model in ["DSI", "DTI", "HARDI"]:
            input_dir = os.path.join(self.base_directory, "RAWDATA", model)
            if len(os.listdir(input_dir)) > 0:
                if convert_rawdata(self.base_directory, input_dir, model):
                    diffusion_available = True
                    diffusion_model.append(model)

        # Check for (and if existing, convert)  T1
        input_dir = os.path.join(self.base_directory, "RAWDATA", "T1")
        if len(os.listdir(input_dir)) > 0:
            if convert_rawdata(self.base_directory, input_dir, "T1_orig"):
                t1_available = True

        # Check for (and if existing, convert)  T2
        input_dir = os.path.join(self.base_directory, "RAWDATA", "T2")
        if len(os.listdir(input_dir)) > 0:
            if convert_rawdata(self.base_directory, input_dir, "T2_orig"):
                t2_available = True

        if diffusion_available:
            # project.stages['Diffusion'].config.imaging_model_choices = diffusion_model
            if t2_available:
                swap_and_reorient(
                    src_file=os.path.join(self.base_directory, "NIFTI", "T2_orig.nii.gz"),
                    ref_file=os.path.join(self.base_directory, "NIFTI", diffusion_model[0] + ".nii.gz"),
                    out_file=os.path.join(self.base_directory, "NIFTI", "T2.nii.gz"),
                )
            if t1_available:
                swap_and_reorient(
                    src_file=os.path.join(self.base_directory, "NIFTI", "T1_orig.nii.gz"),
                    ref_file=os.path.join(self.base_directory, "NIFTI", diffusion_model[0] + ".nii.gz"),
                    out_file=os.path.join(self.base_directory, "NIFTI", "T1.nii.gz"),
                )
                valid_inputs = True
                input_message = "Inputs check finished successfully.\nDiffusion and morphological data available."
            else:
                input_message = "Error during inputs check.\nMorphological data (T1) not available."
        elif t1_available:
            input_message = "Error during inputs check. \nDiffusion data not available (DSI/DTI/HARDI)."
        else:
            input_message = (
                "Error during inputs check. No diffusion or morphological data available in folder "
                + os.path.join(self.base_directory, "RAWDATA")
                + "!"
            )

        imaging_model = diffusion_model[0]

        if gui:
            input_notification = Check_Input_Notification(
                message=input_message, imaging_model_options=diffusion_model, imaging_model=imaging_model
            )
            input_notification.configure_traits()
            self.global_conf.imaging_model = input_notification.imaging_model
            diffusion_file = os.path.join(self.base_directory, "NIFTI", input_notification.imaging_model + ".nii.gz")
            n_vol = nib.load(diffusion_file).shape[3]
            if (
                self.stages["Preprocessing"].config.end_vol == 0
                or self.stages["Preprocessing"].config.end_vol == self.stages["Preprocessing"].config.max_vol
                or self.stages["Preprocessing"].config.end_vol >= n_vol - 1
            ):
                self.stages["Preprocessing"].config.end_vol = n_vol - 1
            self.stages["Preprocessing"].config.max_vol = n_vol - 1
            self.stages["Registration"].config.imaging_model = input_notification.imaging_model
            self.stages["Diffusion"].config.imaging_model = input_notification.imaging_model
        else:
            print input_message
            self.global_conf.imaging_model = imaging_model
            diffusion_file = os.path.join(self.base_directory, "NIFTI", imaging_model + ".nii.gz")
            n_vol = nib.load(diffusion_file).shape[3]
            if (
                self.stages["Preprocessing"].config.end_vol == 0
                or self.stages["Preprocessing"].config.end_vol == self.stages["Preprocessing"].config.max_vol
                or self.stages["Preprocessing"].config.end_vol >= n_vol - 1
            ):
                self.stages["Preprocessing"].config.end_vol = n_vol - 1
            self.stages["Preprocessing"].config.max_vol = n_vol - 1
            self.stages["Registration"].config.imaging_model = imaging_model
            self.stages["Diffusion"].config.imaging_model = imaging_model

        if t2_available:
            self.stages["Registration"].config.registration_mode_trait = [
                "Linear (FSL)",
                "BBregister (FS)",
                "Nonlinear (FSL)",
            ]

        self.fill_stages_outputs()
#.........这里部分代码省略.........

        else:
            print nifti_file, fmri_sessions[session_id]
            shutil.move(nifti_file, fmri_sessions[session_id])
        
        # remove the dicom dirs
        for x in glob.glob(os.path.join(dicom_dir, '*')):
            os.remove(x)
        os.removedirs(dicom_dir)
    
    ##############################################################
    # Preprocessing
    ##############################################################

    ##############################################################
    # Anatomical segmentation (White/Grey matter)
    mem = Memory(base_dir=subject_dir)
    seg = mem.cache(spm.Segment)
    out_seg = seg(data=anat_image,
                  gm_output_type=[True, True, True],
                  wm_output_type=[True, True, True],
                  csf_output_type=[True, True, True])
    sn_file = out_seg.outputs.transformation_mat
    inv_sn_file = out_seg.outputs.inverse_transformation_mat
    gm_image = out_seg.outputs.normalized_gm_image
    native_gm_image = out_seg.outputs.native_gm_image

    shutil.copyfile(native_gm_image, os.path.join(t1_dir,
        '%s_gm_image.nii' % subject))

    ##############################################################
    #  Slice timing correction

"""
================
Realignment demo
================

This example compares standard realignement to realignement with tagging
correction.
"""
# Create a memory context
from nipype.caching import Memory

mem = Memory("/tmp")

# Give the path to the 4D ASL image
raw_asl_file = "/tmp/func.nii"

# Realign with and without tagging correction
from procasl import preprocessing
import numpy as np

realign = mem.cache(preprocessing.Realign)
x_translation = {}
for correct_tagging in [True, False]:
    out_realign = realign(in_file=raw_asl_file, correct_tagging=correct_tagging)
    x_translation[correct_tagging] = np.loadtxt(out_realign.outputs.realignment_parameters)[:, 2]

# Plot x-translation parameters with and without tagging correction
import matplotlib.pylab as plt

plt.figure(figsize=(10, 5))
for correct_tagging, label, color in zip([True, False], ["corrected", "uncorrected"], "rb"):

import nipype.interfaces.spm as spm
from nipype.caching import Memory

from procasl import preprocessing, quantification

# Create a memory context
mem = Memory('/tmp/no_workflow')

# Give data location
func_file = '/tmp/func.nii'
anat_file = '/tmp/anat.nii'

# Set spm paths
matlab_cmd = '/i2bm/local/spm8-standalone/run_spm8.sh ' +\
    '/i2bm/local/spm8-standalone/mcr/v713 script'
spm.SPMCommand.set_mlab_paths(matlab_cmd=matlab_cmd, use_mcr=True)
paths = ['/i2bm/local/spm8-standalone/spm8_mcr/spm8/']  # TODO: check needed

# Get Tag/Control sequence
get_tag_ctl = mem.cache(preprocessing.GetTagControl)
out_get_tag_ctl = get_tag_ctl(in_file=func_file)

# Rescale
rescale = mem.cache(preprocessing.Rescale)
out_rescale = rescale(in_file=out_get_tag_ctl.outputs.tag_ctl_file,
                      ss_tr=35.4, t_i_1=800., t_i_2=1800.)

# Realign to first scan
realign = mem.cache(preprocessing.Realign)
out_realign = realign(
    in_file=out_rescale.outputs.rescaled_file,

def from_mni_to_native(img, sub_id, include_trans=[True, True, True],
                       interpolation='Linear'):
    '''Maps image from native space to mni.

    We assume that the transformation files already exist for the mappings
    between:
    1) mean bold and anatomy
    2) anatomy and oasis template
    3) oasis template and mni template

    The transforms to include are:
    1) From mni to oasis
    2) From oasis to anat
    3) From anat to bold

    The include transforms should be sequential to have meaninful output,
    which means that transformations sequence [True, False, True] is invalid.
    '''
    check = (include_trans == [True, False, True])
    if check:
        raise Exception('Invalid transformation sequence')
    pipeline_dir = 'pipelines/transformations'
    if not os.path.exists(pipeline_dir):
        os.makedirs(pipeline_dir)
    mem = Memory(pipeline_dir)
    transform = mem.cache(ApplyTransforms)

    oasis_template = os.path.join('pipelines',
                                  'OASIS-30_Atropos_template',
                                  'T_template0.nii.gz')
    anat = os.path.join('pipelines',
                        'preprocessing',
                        'sub{0}'.format(sub_id),
                        'highres001.nii')
    mean_bold = os.path.join('pipelines', 'preprocessing',
                             'sub{0}'.format(sub_id),
                             'mean_bold.nii')
    mni_to_oasis = os.path.join('pipelines', 'preprocessing',
                                'registered_templates', 'mni_to_oasis.h5')
    oasis_to_anat = os.path.join('pipelines', 'preprocessing',
                                 'sub{0}'.format(sub_id),
                                 'oasis_to_anat.h5')
    bold_to_anat = os.path.join('pipelines', 'preprocessing',
                                'sub{0}'.format(sub_id),
                                'bold_to_anat.txt')

    all_references = [oasis_template, anat, mean_bold]
    all_trans = [mni_to_oasis, oasis_to_anat, bold_to_anat]
    all_inv_trans = [False, False, True]
    transforms = []
    inv_trans_flags = []
    reference = None
    for idx, flag in enumerate(include_trans):
        if flag:
            transforms.append(all_trans[idx])
            inv_trans_flags.append(all_inv_trans[idx])
            # Use latest transformation as reference
            reference = all_references[idx]

    trans = transform(args='--float',
                      input_image_type=3,
                      interpolation=interpolation,
                      invert_transform_flags=inv_trans_flags[::-1],
                      num_threads=n_proc,
                      reference_image=reference,
                      terminal_output='file',
                      transforms=transforms[::-1],
                      input_image=img)

    return trans.outputs.output_image

    img.to_filename('/home/ys218403/Data/dartel_cache/sub_%s' % name)
    func_file = '/home/ys218403/Data/dartel_cache/sub_%s' % name
    func_niimg = gaelmem.cache(resample_img)(
        nb.load(func_file),
        target_affine=anat_niimg.get_affine(),
        target_shape=anat_niimg.shape)
    func_niimg.to_filename(
        '/home/ys218403/Data/dartel_cache/oversampled_%s' % name)
    resampled_func.append(
        '/home/ys218403/Data/dartel_cache/oversampled_%s' % name)

print resampled_func

cache_dir = "/home/ys218403/Data/dartel_cache"
if not os.path.exists(cache_dir): os.makedirs(cache_dir)
mem = Memory(cache_dir)

tricky_kwargs = {}

dartelnorm2mni_result = mem.cache(spm.DARTELNorm2MNI)(
    apply_to_files=resampled_func[:1],
    flowfield_files=[flow_fields],
    template_file=template_file,
    ignore_exception=False,
    modulate=False,  # don't modulate
    fwhm=0.,  # don't smooth
    **tricky_kwargs)

normalized_func = dartelnorm2mni_result.outputs.normalized_files

# createwarped_result = mem.cache(spm.CreateWarped)(

def do_subject_preproc(subject_id,
                       output_dir,
                       func,
                       anat,
                       do_bet=True,
                       do_mc=True,
                       do_coreg=True,
                       do_normalize=True,
                       cmd_prefix="fsl5.0-",
                       **kwargs
                       ):
    """
    Preprocesses subject data using FSL.

    Parameters
    ----------

    """

    output = {'func': func,
              'anat': anat
              }

    # output dir
    subject_output_dir = os.path.join(output_dir, subject_id)
    if not os.path.exists(subject_output_dir):
        os.makedirs(subject_output_dir)

    # prepare for smart-caching
    cache_dir = os.path.join(output_dir, "cache_dir")
    if not os.path.exists(cache_dir):
        os.makedirs(cache_dir)

    nipype_mem = NipypeMemory(base_dir=cache_dir)
    joblib_mem = JoblibMemory(cache_dir, verbose=100)

    # sanitize input files
    if not isinstance(output['func'], basestring):
        output['func'] = joblib_mem.cache(do_fsl_merge)(
            func, subject_output_dir, output_prefix='Merged',
            cmd_prefix=cmd_prefix)

    ######################
    #  Brain Extraction
    ######################
    if do_bet:
        if not fsl.BET._cmd.startswith("fsl"):
            fsl.BET._cmd = cmd_prefix + fsl.BET._cmd

        bet = nipype_mem.cache(fsl.BET)
        bet_results = bet(in_file=output['anat'],
                          )

        output['anat'] = bet_results.outputs.out_file

    #######################
    #  Motion correction
    #######################
    if do_mc:
        if not fsl.MCFLIRT._cmd.startswith("fsl"):
            fsl.MCFLIRT._cmd = cmd_prefix + fsl.MCFLIRT._cmd

        mcflirt = nipype_mem.cache(fsl.MCFLIRT)
        mcflirt_results = mcflirt(in_file=output['func'],
                                  cost='mutualinfo',
                                  save_mats=True,  # save mc matrices
                                  save_plots=True  # save mc params
                                  )

        output['motion_parameters'] = mcflirt_results.outputs.par_file
        output['motion_matrices'] = mcflirt_results.outputs.mat_file
        output['func'] = mcflirt_results.outputs.out_file

    ###################
    # Coregistration
    ###################
    if do_coreg:
        if not fsl.FLIRT._cmd.startswith("fsl"):
            fsl.FLIRT._cmd = cmd_prefix + fsl.FLIRT._cmd

        flirt1 = nipype_mem.cache(fsl.FLIRT)
        flirt1_results = flirt1(in_file=output['func'],
                                reference=output['anat']
                                )

        if not do_normalize:
            output['func'] = flirt1_results.outputs.out_file

    ##########################
    # Spatial normalization
    ##########################
    if do_normalize:
        if not fsl.FLIRT._cmd.startswith("fsl"):
            fsl.FLIRT._cmd = cmd_prefix + fsl.FLIRT._cmd

        # T1 normalization
        flirt2 = nipype_mem.cache(fsl.FLIRT)
        flirt2_results = flirt2(in_file=output['anat'],
                                reference=FSL_T1_TEMPLATE)

#.........这里部分代码省略.........

    def check_input(self, gui=True):
        print '**** Check Inputs ****'
        fMRI_available = False
        t1_available = False
        t2_available = False
        valid_inputs = False

        mem = Memory(base_dir=os.path.join(self.base_directory,'NIPYPE'))
        swap_and_reorient = mem.cache(SwapAndReorient)

        # Check for (and if existing, convert) functional data
        input_dir = os.path.join(self.base_directory,'RAWDATA','fMRI')
        if len(os.listdir(input_dir)) > 0:
            if convert_rawdata(self.base_directory, input_dir, 'fMRI'):
                fMRI_available = True

        # Check for (and if existing, convert)  T1
        input_dir = os.path.join(self.base_directory,'RAWDATA','T1')
        if len(os.listdir(input_dir)) > 0:
            if convert_rawdata(self.base_directory, input_dir, 'T1_orig'):
                t1_available = True

        # Check for (and if existing, convert)  T2
        input_dir = os.path.join(self.base_directory,'RAWDATA','T2')
        if len(os.listdir(input_dir)) > 0:
            if convert_rawdata(self.base_directory, input_dir, 'T2_orig'):
                t2_available = True   

        if fMRI_available:
            if t2_available:
                swap_and_reorient(src_file=os.path.join(self.base_directory,'NIFTI','T2_orig.nii.gz'),
                                  ref_file=os.path.join(self.base_directory,'NIFTI','fMRI.nii.gz'),
                                  out_file=os.path.join(self.base_directory,'NIFTI','T2.nii.gz'))
            if t1_available:
                swap_and_reorient(src_file=os.path.join(self.base_directory,'NIFTI','T1_orig.nii.gz'),
                                  ref_file=os.path.join(self.base_directory,'NIFTI','fMRI.nii.gz'),
                                  out_file=os.path.join(self.base_directory,'NIFTI','T1.nii.gz'))
                valid_inputs = True
                input_message = 'Inputs check finished successfully.\nfMRI and morphological data available.'
            else:
                input_message = 'Error during inputs check.\nMorphological data (T1) not available.'
        elif t1_available:
            input_message = 'Error during inputs check. \nfMRI data not available (fMRI).'
        else:
            input_message = 'Error during inputs check. No fMRI or morphological data available in folder '+os.path.join(self.base_directory,'RAWDATA')+'!'

        if gui: 
            input_notification = Check_Input_Notification(message=input_message, imaging_model='fMRI')
            input_notification.configure_traits()
            self.global_conf.imaging_model = input_notification.imaging_model
            self.stages['Registration'].config.imaging_model = input_notification.imaging_model
        else:
            print input_message
            self.global_conf.imaging_model = 'fMRI'
            self.stages['Registration'].config.imaging_model = 'fMRI'
       
        if t2_available:
            self.stages['Registration'].config.registration_mode_trait = ['Linear (FSL)','BBregister (FS)','Nonlinear (FSL)']
       
        self.fill_stages_outputs()
       
        return valid_inputs