def direct_nifti_to_directory(base_directory, dicom_header, niftis, bvals=None, bvecs=None):
    """Move nifti file to the correct directory for the subject
    @param dicom_headers: dicom_header object created by dicom.read_file and stored in a pickle dump
    @param nifti_file: a list of nifti files to move
    @param subject_directory: string representing the main directory for the subject

    @return nifti_destination: string representing where the file moved to

    """
    import dicom
    import os
    
    import nipype.utils.filemanip
    
    ANATOMY = 0
    BOLD = 1
    DTI = 2
    FIELDMAP = 3
    LOCALIZER = 4
    REFERENCE = 5
    DERIVED = 6
    
    from nipype import logging
    iflogger = logging.getLogger('interface')
    wlogger = logging.getLogger('workflow')
    
    iflogger.debug('Enetering direct nifti with inputs {} {} {} '.format(dicom_header, niftis, base_directory))
    
    dicom_header_filename = dicom_header
    dicom_header = nipype.utils.filemanip.loadpkl(dicom_header)
    
    scan_keys = [['MPRAGE','FSE','T1w','T2w','PDT','PD-T2','tse2d','t2spc','t2_spc'],
        ['epfid'],
        ['ep_b'],
        ['fieldmap','field_mapping'],
        ['localizer','Scout'],
        ['SBRef']]
        
    file_type = None 

    destination = ""
    
    if type(niftis) is str:
        niftis = [niftis]
    
    try:
        for i in range(0,6):
            for key in scan_keys[i]:
                if (dicom_header.ProtocolName.lower().find(key.lower()) > -1) or \
                    (dicom_header.SeriesDescription.lower().find(key.lower()) > -1) or \
                    (dicom_header.SequenceName.lower().find(key.lower()) > -1):
                        file_type = i
                        
        if dicom_header.ImageType[0] != "ORIGINAL":
            file_type = DERIVED
                
    except AttributeError, e:
        iflogger.warning("Nifti File(s) {} not processed because Dicom header dataset throwing error {} \n ".format(niftis, e)
                         + "Check your dicom headers")
        return []

def calc_surface_potential_fn(dipole_file, dipole_row, leadfield, mesh_file, mesh_id):
    import os.path as op
    import numpy as np
    import h5py
    from forward.mesh import get_closest_element_to_point

    from nipype import logging
    iflogger = logging.getLogger('interface')

    dipole_data = np.loadtxt(dipole_file)
    if len(dipole_data.shape) == 1:
        dipole_data = [dipole_data]
    dip = dipole_data[int(dipole_row)]
    
    x, y, z = [float(i) for i in dip[2:5]]
    q_x, q_y, q_z = [float(j) for j in dip[6:9]]
    _, element_idx, centroid, element_data, lf_idx = get_closest_element_to_point(mesh_file, mesh_id, [[x, y, z]])

    lf_data_file = h5py.File(leadfield, "r")
    lf_data = lf_data_file.get("leadfield")
    leadfield_matrix = lf_data.value

    L = np.transpose(leadfield_matrix[lf_idx * 3:lf_idx * 3 + 3])
    J = np.array([q_x,q_y,q_z])
    potential = np.dot(L,J)
    out_potential = op.abspath("potential.npy")
    np.save(out_potential,potential)
    return out_potential

def World(in_file, some_parameter):
   from nipype import logging
   iflogger = logging.getLogger('interface')
   message = "World! " + 'some_parameter: ' + str(some_parameter)
   iflogger.info(message)
   with open(in_file, 'a') as fp:
       fp.write(message)

def return_subject_data(subject_id, data_file):
    import csv
    from nipype import logging
    iflogger = logging.getLogger('interface')

    f = open(data_file, 'r')
    csv_line_by_line = csv.reader(f)
    found = False
    # Must be stored in this order!:
    # 'subject_id', 'dose', 'weight', 'delay', 'glycemie', 'scan_time']
    for line in csv_line_by_line:
        if line[0] == subject_id:
            dose, weight, delay, glycemie, scan_time = [
                float(x) for x in line[1:]]
            iflogger.info('Subject %s found' % subject_id)
            iflogger.info('Dose: %s' % dose)
            iflogger.info('Weight: %s' % weight)
            iflogger.info('Delay: %s' % delay)
            iflogger.info('Glycemie: %s' % glycemie)
            iflogger.info('Scan Time: %s' % scan_time)
            found = True
            break
    if not found:
        raise Exception("Subject id %s was not in the data file!" % subject_id)
    return dose, weight, delay, glycemie, scan_time

def Hello():
   import os
   from nipype import logging
   iflogger = logging.getLogger('interface')
   message = "Hello "
   file_name =  'hello.txt'
   iflogger.info(message)
   with open(file_name, 'w') as fp:
       fp.write(message)
   return os.path.abspath(file_name)

def swap_element_ids(mesh_filename, elem_list, new_elem_id, new_phys_id):
    import time
    import os.path as op
    from nipype.utils.filemanip import split_filename
    from nipype import logging
    iflogger = logging.getLogger('interface')

    iflogger.info("Reading mesh file: %s" % mesh_filename)

    start_time = time.time()
    mesh_file = open(mesh_filename, 'r')
    _, name, _ = split_filename(mesh_filename)
    out_file = op.abspath(name + "_mask.msh")

    f = open(out_file, 'w')
    while True:
        line = mesh_file.readline()
        f.write(line)

        if line == '$Nodes\n':
            line = mesh_file.readline()
            f.write(line)
            number_of_nodes = int(line)
            iflogger.info("%d nodes in mesh" % number_of_nodes)

            for i in xrange(0, number_of_nodes):
                line = mesh_file.readline()
                f.write(line)

        elif line == '$Elements\n':
            line = mesh_file.readline()
            f.write(line)
            number_of_elements = int(line)
            iflogger.info("%d elements in mesh" % number_of_elements)
            elem_lines = []
            for i in xrange(0, number_of_elements):
                # -- If all elements were quads, each line has 10 numbers.
                #    The first one is a tag, and the last 4 are node numbers.
                line = mesh_file.readline()
                elem_lines.append(line)
                elem_data = line.split()
                if int(elem_data[0]) in elem_list:
                    elem_data[3] = str(new_phys_id)
                    elem_data[4] = str(new_elem_id)
                line = " ".join(elem_data) + "\n"
                f.write(line)

        elif line == '$EndElementData\n' or len(line) == 0:
            break

    mesh_file.close()
    f.close()
    elapsed_time = time.time() - start_time
    print(elapsed_time)
    return out_file

def parse_and_return_mats(one_d_file, mask_arr):
    '''
    '''

    # Import packages
    import numpy as np
    import scipy.sparse as sparse
    from nipype import logging

    # Init logger
    logger = logging.getLogger('workflow')

    # Parse out numbers
    logger.info('Parsing contents...')
    graph_arr = np.loadtxt(one_d_file, skiprows=6)

    # Cast as numpy arrays and extract i, j, w
    logger.info('Creating arrays...')
    one_d_rows = graph_arr.shape[0]

    # Extract 3d indices
    ijk1 = graph_arr[:, 2:5].astype('int32')
    ijk2 = graph_arr[:, 5:8].astype('int32')
    # Weighted array and binarized array
    w_arr = graph_arr[:,-1].astype('float32')
    del graph_arr
    b_arr = np.ones(w_arr.shape)

    # Non-zero elements from mask is size of similarity matrix
    mask_idx = np.argwhere(mask_arr)
    mask_voxs = mask_idx.shape[0]

    # Extract the ijw's from 1D file
    i_arr = [np.where((mask_idx == ijk1[ii]).all(axis=1))[0][0] \
             for ii in range(one_d_rows)]
    del ijk1
    j_arr = [np.where((mask_idx == ijk2[ii]).all(axis=1))[0][0] \
             for ii in range(one_d_rows)]
    del ijk2
    i_arr = np.array(i_arr, dtype='int32')
    j_arr = np.array(j_arr, dtype='int32')

    # Construct the sparse matrix
    logger.info('Constructing sparse matrix...')
    wmat_upper_tri = sparse.coo_matrix((w_arr, (i_arr, j_arr)),
                                       shape=(mask_voxs, mask_voxs))
    bmat_upper_tri = sparse.coo_matrix((b_arr, (i_arr, j_arr)),
                                       shape=(mask_voxs, mask_voxs))

    # Make symmetric
    w_similarity_matrix = wmat_upper_tri + wmat_upper_tri.T
    b_similarity_matrix = bmat_upper_tri + bmat_upper_tri.T

    # Return the symmetric matrices and affine
    return b_similarity_matrix, w_similarity_matrix

def calc_tpm_fn(tracks):
   import os
   from nipype import logging
   from nipype.utils.filemanip import split_filename
   path, name, ext = split_filename(tracks)
   file_name = os.path.abspath(name + 'TPM.nii')
   iflogger = logging.getLogger('interface')
   iflogger.info(tracks)
   import subprocess
   iflogger.info(" ".join(["tracks2prob","-vox", "1.0", "-totallength", tracks, file_name]))
   subprocess.call(["tracks2prob", "-vox", "1.0", "-totallength", tracks, file_name])
   return file_name

def split_warp_volumes_fn(in_file):
    from nipype import logging
    from nipype.utils.filemanip import split_filename
    import nibabel as nb
    import os.path as op
    iflogger = logging.getLogger('interface')
    iflogger.info(in_file)
    path, name, ext = split_filename(in_file)
    image = nb.load(in_file)
    x_img, y_img, z_img = nb.four_to_three(image)
    x = op.abspath(name + '_x' + ".nii.gz")
    y = op.abspath(name + '_y' + ".nii.gz")
    z = op.abspath(name + '_z' + ".nii.gz")
    nb.save(x_img, x)
    nb.save(y_img, y)
    nb.save(z_img, z)
    return x, y, z

def setup_logger(logger_name, file_path, level, to_screen=False):
    '''
    Function to initialize and configure a logger that can write to file
    and (optionally) the screen.

    Parameters
    ----------
    logger_name : string
        name of the logger
    file_path : string
        file path to the log file on disk
    level : integer
        indicates the level at which the logger should log; this is
        controlled by integers that come with the python logging
        package. (e.g. logging.INFO=20, logging.DEBUG=10)
    to_screen : boolean (optional)
        flag to indicate whether to enable logging to the screen

    Returns
    -------
    logger : logging.Logger object
        Python logging.Logger object which is capable of logging run-
        time information about the program to file and/or screen
    '''

    # Import packages
    import logging

    # Init logger, formatter, filehandler, streamhandler
    logger = logging.getLogger(logger_name)
    logger.setLevel(level)
    formatter = logging.Formatter('%(asctime)s : %(message)s')

    # Write logs to file
    fileHandler = logging.FileHandler(file_path)
    fileHandler.setFormatter(formatter)
    logger.addHandler(fileHandler)

    # Write to screen, if desired
    if to_screen:
        streamHandler = logging.StreamHandler()
        streamHandler.setFormatter(formatter)
        logger.addHandler(streamHandler)

    # Return the logger
    return logger

def get_dicom_headers(dicom_file):
    import dicom
    import nipype.utils.filemanip
    from nipype import logging
    
    if type(dicom_file) is list:
        dicom_file = dicom_file.pop()
    
    iflogger = logging.getLogger('interface')
    iflogger.debug('Getting headers from {}'.format(dicom_file))
    headers = dicom.read_file(dicom_file)
    
    
    
    iflogger.debug('Returning headers from {}'.format(dicom_file))
    nipype.utils.filemanip.savepkl(dicom_file + '.pklz', headers)
    
    return dicom_file + '.pklz'

def extract_frommask_component(realigned_file, mask_file):
    import os
    import nibabel as nb
    import numpy as np
    from utils import mean_roi_signal
    from nipype import logging
    iflogger = logging.getLogger('interface')

    data = nb.load(realigned_file).get_data().astype('float64')
    mask = nb.load(mask_file).get_data().astype('float64')
    iflogger.info('Data and mask loaded.')
    mask_comp = mean_roi_signal(data, mask.astype('bool'))

    components_file = os.path.join(os.getcwd(), 'mask_mean_component.txt')
    iflogger.info('Saving components file:' + components_file)
    np.savetxt(components_file, mask_comp)

    return components_file

def degree_centrality(corr_matrix, r_value, method, out=None):
    """
    Calculate centrality for the rows in the corr_matrix using
    a specified correlation threshold. The centrality output can 
    be binarized or weighted.
    
    Paramaters
    ---------
    corr_matrix : numpy.ndarray
    r_value : float
    method : str
        Can be 'binarize' or 'weighted'
    out : numpy.ndarray (optional)
        If specified then should have shape of `corr_matrix.shape[0]`
    
    Returns
    -------
    out : numpy.ndarray
    """

    # Import packages
    from nipype import logging

    # Init logger
    logger = logging.getLogger('workflow')

    if method not in ["binarize", "weighted"]:
        raise Exception("Method must be one of binarize or weighted and not %s" % method)
    
    if corr_matrix.dtype.itemsize == 8:
        dtype   = "double"
        r_value = np.float64(r_value)
    else:
        dtype   = "float"
        r_value = np.float32(r_value)
    
    if out is None:
        out = np.zeros(corr_matrix.shape[0], dtype=corr_matrix.dtype)
    logger.info('about to call thresh_and_sum')
    func_name   = "centrality_%s_%s" % (method, dtype)
    func        = globals()[func_name]
    func(corr_matrix, out, r_value)
    
    return out

def clean_warp_field_fn(combined_warp_x, combined_warp_y, combined_warp_z, default_value):
    import os.path as op
    from nipype import logging
    from nipype.utils.filemanip import split_filename
    import nibabel as nb
    import numpy as np
    path, name, ext = split_filename(combined_warp_x)
    out_file = op.abspath(name + 'CleanedWarp.nii')
    iflogger = logging.getLogger('interface')
    iflogger.info(default_value)
    imgs = []
    filenames = [combined_warp_x, combined_warp_y, combined_warp_z]
    for fname in filenames:
        img = nb.load(fname)
        data = img.get_data()
        data[data==default_value] = np.NaN
        new_img = nb.Nifti1Image(data=data, header=img.get_header(), affine=img.get_affine())
        imgs.append(new_img)
    image4d = nb.concat_images(imgs, check_affines=True)
    nb.save(image4d, out_file)
    return out_file

def read_mesh_elem_data(mesh_filename, view_name="Cost function"):
    import numpy as np

    from nipype import logging
    iflogger = logging.getLogger('interface')

    iflogger.info("Reading mesh file: %s" % mesh_filename)

    mesh_file = open(mesh_filename, 'r')
    while True:
        line = mesh_file.readline()
        if '$ElementData' in line:
            line = mesh_file.readline()
            name = mesh_file.readline()
            if name == '"' + view_name + '"\n':
                line = mesh_file.readline()
                line = mesh_file.readline()
                line = mesh_file.readline()
                line = mesh_file.readline()
                line = mesh_file.readline()
                number_of_elements = int(mesh_file.readline().replace("\n",""))

                elem_data = []
                iflogger.info("%d elements in element data" % number_of_elements)
                for i in xrange(0, number_of_elements):
                    line = mesh_file.readline()
                    line_data = line.split()
                    polygon = {}
                    polygon["element_id"] = int(line_data[0])
                    polygon["data"] = float(line_data[1])
                    elem_data.append(polygon)

                iflogger.info("Done reading element data")
                break
    
    # Loop through and assign points to each polygon, save as a dictionary
    num_polygons = len(elem_data)
    iflogger.info("%d polygons found" % num_polygons)

    return elem_data

def register_template(hemi, sphere_file, transform, templates_path, template):
    """
    Register surface to template with FreeSurfer's mris_register.

    Transform the labels from multiple atlases via a template
    (using FreeSurfer's mris_register).

    """
    from os import path
    from nipype.interfaces.base import CommandLine
    from nipype import logging

    logger = logging.getLogger("interface")

    template_file = path.join(templates_path, hemi + "." + template)
    output_file = hemi + "." + transform
    cli = CommandLine(command="mris_register")
    cli.inputs.args = " ".join(["-curv", sphere_file, template_file, output_file])
    logger.info(cli.cmdline)
    cli.run()

    return transform

def rewrite_mesh_from_binary_mask(mask_file, mesh_file, mesh_id, new_phys_id=1015):
    # Takes about 8 minutes
    import numpy as np
    import nibabel as nb
    import os.path as op
    from forward.mesh import read_mesh
    from nipype import logging

    iflogger = logging.getLogger('interface')

    mask_img = nb.load(mask_file)
    mask_data = mask_img.get_data()
    mask_affine = mask_img.get_affine()
    mask_header = mask_img.get_header()
    #tensor_data = np.flipud(tensor_image.get_data())

    # Define various constants
    elements_to_consider = [1002] #Use only white matter
    vx, vy, vz = mask_header.get_zooms()[0:3]
    max_x, max_y, max_z = np.shape(mask_data)[0:3]
    halfx, halfy, halfz = np.array((vx*max_x, vy*max_y, vz*max_z))/2.0

    mesh_data, _, _, _ = read_mesh(mesh_file, elements_to_consider)

    elem_list = []
    for idx, poly in enumerate(mesh_data):
        i = np.round((poly['centroid'][0]+halfx)/vx).astype(int)
        j = np.round((poly['centroid'][1]+halfy)/vy).astype(int)
        k = np.round((poly['centroid'][2]+halfz)/vz).astype(int)
        if mask_data[i,j,k]:
            elem_list.append(poly['element_id'])
        #iflogger.info("%3.3f%%" % (float(idx)/num_polygons*100.0))
    new_elem_id = new_phys_id+2000
    out_file = swap_element_ids(mesh_file, elem_list, new_elem_id, new_phys_id)
    print('Writing binary masked mesh file to %s' % out_file)
    return out_file

def include_gmsh_tensor_elements(mesh_file, tensor_file, mask_file, mask_threshold=0.5, lower_triangular=True):
    import numpy as np
    import nibabel as nb
    from nipype.utils.filemanip import split_filename
    import os.path as op
    from forward.mesh import read_mesh
    from nipype import logging
    import shutil

    iflogger = logging.getLogger("interface")

    # Load 4D (6 volume upper or lower triangular) conductivity tensor image
    tensor_image = nb.load(tensor_file)
    tensor_data = np.flipud(tensor_image.get_data())

    # Correct the tensors after flipping in the X direction
    tensor_data[:, :, :, 1] *= -1
    if lower_triangular:
        tensor_data[:, :, :, 3] *= -1
    else:
        tensor_data[:, :, :, 2] *= -1

    # Load mask (usually fractional anisotropy) image
    mask_image = nb.load(mask_file)
    mask_data = np.flipud(mask_image.get_data())
    header = tensor_image.get_header()

    # Make sure the files share the same (3D) dimensions before continuing
    assert np.shape(tensor_data)[0:3] == np.shape(mask_data)[0:3]

    # Define various constants
    elements_to_consider = [1001]  # Use only white matter
    vx, vy, vz = header.get_zooms()[0:3]
    max_x, max_y, max_z = np.shape(tensor_data)[0:3]
    halfx, halfy, halfz = np.array((vx * max_x, vy * max_y, vz * max_z)) / 2.0

    mesh_data, _, _, _ = read_mesh(mesh_file, elements_to_consider)

    # Create the output mesh file
    path, name, ext = split_filename(mesh_file)
    out_file = op.abspath(name + "_cond.msh")
    iflogger.info("Copying current mesh file to %s" % out_file)
    shutil.copyfile(mesh_file, out_file)

    f = open(out_file, "a")  # Append to the end of the file
    iflogger.info("Appending Conductivity tensors to %s" % out_file)

    # Write the tag information to the file:
    num_polygons = len(mesh_data)
    f.write("$ElementData\n")
    str_tag = '"Conductivity"'
    timestep = 0.0001

    f.write("1\n")  # Num String tags
    f.write(str_tag + "\n")
    f.write("1\n")  # Num Real tags
    f.write("%f\n" % timestep)

    # Three integer tags: timestep, num field components, num elements
    f.write("3\n")  # Three int tags
    f.write("0\n")  # Time step index
    f.write("9\n")  # Num field components

    # Get the centroid of all white matter elements
    # Find out which voxel they lie inside
    iflogger.info("Getting tensor for each element")
    # ipdb.set_trace()
    nonzero = 0
    elem_list = []

    if lower_triangular:
        for idx, poly in enumerate(mesh_data):
            i = np.round((poly["centroid"][0] + halfx) / vx).astype(int)
            j = np.round((poly["centroid"][1] + halfy) / vy).astype(int)
            k = np.round((poly["centroid"][2] + halfz) / vz).astype(int)
            T = tensor_data[i, j, k]
            if not (all(T == 0) and mask_data[i, j, k] >= mask_threshold):
                elementdata_str = "%d %e %e %e %e %e %e %e %e %e\n" % (
                    poly["element_id"],
                    T[0],
                    T[1],
                    T[3],
                    T[1],
                    T[2],
                    T[4],
                    T[3],
                    T[4],
                    T[5],
                )
                elem_list.append(elementdata_str)
                nonzero += 1
            # iflogger.info("%3.3f%%" % (float(idx)/num_polygons*100.0))
    else:
        for idx, poly in enumerate(mesh_data):
            i = np.round((poly["centroid"][0] + halfx) / vx).astype(int)
            j = np.round((poly["centroid"][1] + halfy) / vy).astype(int)
            k = np.round((poly["centroid"][2] + halfz) / vz).astype(int)
            T = tensor_data[i, j, k]
            if not (all(T == 0) and mask_data[i, j, k] >= mask_threshold):
                elementdata_str = "%d %e %e %e %e %e %e %e %e %e\n" % (
#.........这里部分代码省略.........

import nipype.interfaces.utility as util  # utility
import nipype.pipeline.engine as pe  # pypeline engine
import nipype.interfaces.fsl as fsl
import nipype.interfaces.freesurfer as fs
import nipype.interfaces.dipy as dipy
import os, os.path as op

from nipype import logging

iflogger = logging.getLogger("interface")

fsl.FSLCommand.set_default_output_type("NIFTI_GZ")


def nifti_tensors_to_gmsh(in_file, fa_file, threshold=0.8, lower_triangular=True):
    import numpy as np
    import nibabel as nb
    from nipype.utils.filemanip import split_filename
    import os.path as op

    tensor_image = nb.load(in_file)
    fa_image = nb.load(fa_file)
    path, name, ext = split_filename(in_file)
    out_file = op.abspath(name + ".pos")
    f = open(out_file, "w")
    print "Writing tensors to {f}".format(f=out_file)
    tensor_data = tensor_image.get_data()
    orig_t = tensor_data.copy()
    tensor_data = np.flipud(tensor_data)

    tensor_data[:, :, :, 1] *= -1

import nipype.pipeline.engine as pe
import nipype.interfaces.utility as util
import nipype.interfaces.diffusion_toolkit as dtk
import nipype.interfaces.fsl as fsl
import nipype.interfaces.freesurfer as fs
import nipype.interfaces.mrtrix as mrtrix
import nipype.interfaces.camino as camino
from nipype.utils.filemanip import split_filename
import nibabel as nib

from nipype.interfaces.base import CommandLine, CommandLineInputSpec,\
    traits, TraitedSpec, BaseInterface, BaseInterfaceInputSpec
import nipype.interfaces.base as nibase

from nipype import logging
iflogger = logging.getLogger('interface')

# Reconstruction configuration
    
class DTK_recon_config(HasTraits):
    imaging_model = Str
    maximum_b_value = Int(1000)
    gradient_table_file = Enum('siemens_06',['mgh_dti_006','mgh_dti_018','mgh_dti_030','mgh_dti_042','mgh_dti_060','mgh_dti_072','mgh_dti_090','mgh_dti_120','mgh_dti_144',
                          'siemens_06','siemens_12','siemens_20','siemens_30','siemens_64','siemens_256','Custom...'])
    gradient_table = Str
    custom_gradient_table = File
    flip_table_axis = List(editor=CheckListEditor(values=['x','y','z'],cols=3))
    dsi_number_of_directions = Enum([514,257,124])
    number_of_directions = Int(514)
    number_of_output_directions = Int(181)
    recon_matrix_file = Str('DSI_matrix_515x181.dat')