gin commit from rMBP-15-Malthe.local

Deleted files: 6

data/template_creation_pipline/README.md

@@ -1,47 +0,0 @@
-# DAPI template pipeline
-The original pipeline for pre-processing, reconstructing, and template
-creation of the DAPI template.
-
-## Description
-The reconstruction and template creation was run twice.
-First with the allen reference atlas as reference without non-linear
-registration steps resulting in an initial DAPI template. The initial DAPI
-template was used as reference (instead of the Allen reference atlas)
-in the second run of reconstruction and template creation, resulting in the
-final DAPI template. The two runs of the pipeline were manually modified.
-
-This pipeline is based on the original data structure of the slice images
-with a few sanity checks to make sure that all images are present. To run this
-pipeline with the data structure from the GIN-repository you need to remove
-these sanity checks, get the mouse ids from a different directory (like tif/),
-and get all the image paths in each individual tif folder.
-
-## Overview
-The main script calls the pre_processing, co_registration, and
-template_creation with appropriate parameters and paths. In pre_processing
-the original tif image is down-sampled to 8.6x8.6 um and masked according to
-the biggest (except background) connected component after smoothing and
-threshold to remove artefact outside of the brain slice.
-After pre-processing the slices for each brain is reconstructed
-into 3-dimensional brain volumes. This reconstruction is carried out by 
-iterative registrations (both between volumes and between slice-pairs) to
-a reference (Allen or DAPI from first iteration). Finally, each reconstructed
-brain volume is used in the creation of a population-based average in 
-template_creation.
-
-## Packages
-This pipeline uses the following non-standard packages (package name: version):
-
-```
-nibabel: 2.1.0
-nipype: 0.13.1
-numpy: 1.15.1
-PIL: 5.1.0
-scipy: 1.1.0
-skimage: 0.14.0
-sklearn: 0.19.2
-```
-
-## Registration software
-A precompiled binary of Advanced Normalization Tools (ANTs) V.2.1.0 for 
-linux/mac was used.

data/template_creation_pipline/co_registration.py

@@ -1,288 +0,0 @@
-from helper_functions import pretty, make_dirs, slice_volume, files_in_dir, files_ind_dir_wo_string, make_volume
-from nipype.interfaces import ants
-from os import path
-from shutil import rmtree
-import sys
-import getopt
-
-# Help message.
-tab = '     '
-usage_message = '\nCOMMAND:\n' + tab + sys.argv[0] + \
-                '\n\nOPTIONS:\n' + \
-                tab + '-h, help menu\n' + \
-                tab + '-e, --template, <input template>\n' + \
-                tab + '-s, --slices, <directory with input slices>\n' + \
-                tab + '-r, --rawvolume, <raw volume of slices>\n' + \
-                tab + '-t, --inittrans, <initial transform, roughly transforming raw slice volume to match input ' \
-                      'volume, manually created>\n' + \
-                tab + '-d, --id, <brain id>\n' + \
-                tab + '-o, --output, <output directory, will create new directory if not present>\n' + \
-                tab + '-l, --log, [bool]>\n'
-
-
-def registration(fixed, moving, result, reg_type, sigma, shrink, param, dim, initial_trans='', inverse_result=''):
-    """
-    Registration function. Registeres the moving image to the fixed with given parameters. If initial transform is not
-    given then the images are aligned by center of mass. Outputs inversely transformed result if given.
-    :param fixed: Path to fixed image.
-    :param moving: Path to moving image.
-    :param result: Path to result image.
-    :param reg_type: Registration type. List containing, Rigid, Affine, and/or SyN.
-    :param sigma: List of sigma parameters. Ordered according to the reg_type parameter.
-    :param shrink: List of shrink parameters. Ordered according to the reg_type parameter.
-    :param param: List of transform parameters. Ordered according to the reg_type parameter.
-    :param dim: Dimension of images. Typically 2 or 3.
-    :param initial_trans: Optional path to initial moving transform. If not given, the images will be matched initially
-    by aligning the center of mass.
-    :param inverse_result: Optional path to the inversely transformed fixed image.
-    :return: Returns nothing.
-    """
-    # Extract number of registrations.
-    steps = len(reg_type)
-
-    # Add 1000 iterations for each step.
-    if steps == 1:
-        iterations = [[1000] * len(sigma[0])]
-    else:
-        iterations = []
-        for i, reg in enumerate(reg_type):
-            iteration = [1000] * len(sigma[i])
-            iterations.append(iteration)
-
-    # Create registration instance.
-    reg_instance = ants.Registration(dimension=dim,
-                                     transforms=reg_type,
-                                     transform_parameters=param,
-                                     metric=['MeanSquares']*steps,
-                                     interpolation='BSpline',
-                                     fixed_image=[fixed],
-                                     moving_image=[moving],
-                                     metric_weight=[1]*steps,
-                                     radius_or_number_of_bins=[32]*steps,
-                                     number_of_iterations=iterations,
-                                     smoothing_sigmas=sigma,
-                                     shrink_factors=shrink,
-                                     convergence_threshold=[1.e-7],
-                                     sampling_strategy=['Regular']*steps,
-                                     use_histogram_matching=True,
-                                     winsorize_lower_quantile=0.05,
-                                     winsorize_upper_quantile=0.95,
-                                     write_composite_transform=True,
-                                     output_warped_image=result,
-                                     output_transform_prefix=result.split('.')[0] + '-Transformation',
-                                     num_threads=12)
-
-    # Add initial moving transform if given, else match by center of mass.
-    if not initial_trans == '':
-        reg_instance.inputs.initial_moving_transform = initial_trans
-    else:
-        reg_instance.inputs.initial_moving_transform_com = 1
-
-    # Output reverse results if path is given.
-    if not inverse_result == '':
-        reg_instance.inputs.output_inverse_warped_image = inverse_result
-
-    # Run registration.
-    reg_instance.run()
-
-
-def co_register(template, slices, raw_volume, initial_transform, out_dir, brain_id, r=3, a=2, nl=3, print_log=False):
-    """
-    Function for reconstructing brain volume by iteratively registering to given template. The slices are stacked into a
-    3-dimensional volume that is registered to the template. The template is inversely transformed and sliced, resulting
-    in slice pairs with a template slice and an input slice. For each pair, the input slice is registered to the
-    template slice, and all the resulting registered slices are stacked to a volume. This process is repeated r times
-    rigid, a times rigid and affine, and nl times rigid, affine, and non-linear. Returns the path to the final
-    reconstructed brain volume.
-    :param template: Path to template used as a reference for brain reconstruction.
-    :param slices: List of paths to slices that needs to be reconstructed.
-    :param raw_volume: Path to the initial stacked volume, needed for first rigid registration step.
-    :param initial_transform:Initial 3-dimensional transform of the raw volume to the template. If empty string given,
-    they will be matched by the center of mass.
-    :param out_dir: Output directory.
-    :param brain_id: Mouse id.
-    :param r: Int [Default: 3]. Number of rigid registrations.
-    :param a: Int [Default: 2]. Number of rigid and affine registrations.
-    :param nl: Int [Default: 3]. Number of rigid, affine, and non-linear registrations.
-    :param print_log: Boolean [Default: False]. Will print log if true.
-    :return: Return path to the final reconstructed brain volume.
-    """
-
-    # Determine whether the reference template is the Allen atlas (1. iteration) or the DAPI-template (2. iteration).
-    if 'allen' in template:
-        template_type = 'allen'
-    else:
-        template_type = 'dapi'
-    volume_postfix = '_hv_reg_' + template_type + '_' + brain_id
-
-    # List of slices.
-    raw_slices = files_in_dir(slices, '.nii')
-
-    if print_log:
-        print('Volume given interpreted as:', template_type)
-        print('Checking paths.')
-
-    if not path.exists(raw_volume):
-        print('The path:', raw_volume, 'does not exist.')
-    if not path.exists(template):
-        print('The path:', template, 'does not exist.')
-    if not path.exists(initial_transform):
-        print('The path:', initial_transform, 'does not exist.')
-
-    if print_log:
-        print('Defining iterations.')
-    reg_types = ['Rigid'] * r + ['Affine'] * a + ['SyN'] * nl
-
-    if print_log:
-        print(reg_types)
-
-    # Run through iterations.
-    for i, iteration in enumerate(reg_types):
-        cur_work_dir = out_dir + pretty(i) + '/'
-        make_dirs(cur_work_dir)
-
-        if print_log:
-            print('Running iteration', pretty(i))
-            print('Current work dir:', cur_work_dir)
-
-        if i == 0:
-            hv = raw_volume
-        else:
-            hv = out_dir + pretty(i) + volume_postfix + '.nii'
-            if print_log:
-                print('Deleting prev. iteration work dir:', out_dir + pretty(i - 1) + '/')
-            rmtree(out_dir + pretty(i - 1) + '/')
-
-        # Iteration specific registration parameters.
-        if iteration == 'Rigid':
-            reg = ['Rigid']
-            params = [(0.25,)]
-            vol_sigma = [[4, 4, 2, 2]]
-            vol_shrink = [[32, 16, 8, 4]]
-            slice_sigma = [[4, 4, 2, 2]]
-            slice_shrink = [[32, 16, 8, 4]]
-            if i == 0:
-                vol_sigma = [[4, 4, 4]]
-                vol_shrink = [[32, 16, 8]]
-                slice_sigma = [[4, 4, 2]]
-                slice_shrink = [[32, 16, 8]]
-        elif iteration == 'Affine':
-            reg = ['Rigid', 'Affine']
-            params = [(0.25,), (0.25,), ]
-            vol_sigma = [[4, 4, 2, 2], [4, 4, 2, 2]]
-            vol_shrink = [[32, 16, 8, 4], [32, 16, 8, 4]]
-            slice_sigma = [[4, 4, 2, 2], [4, 4, 2, 2]]
-            slice_shrink = [[32, 16, 8, 4], [32, 16, 8, 4]]
-        elif iteration == 'SyN':
-            reg = ['Rigid', 'Affine', 'SyN']
-            params = [(0.25,), (0.25,), (0.15, 3, 0.5)]
-            vol_sigma = [[4, 4, 2, 2], [4, 4, 2, 2], [4, 4, 4, 2]]
-            vol_shrink = [[32, 16, 8, 4], [32, 16, 8, 4], [64, 32, 16, 8]]
-            slice_sigma = [[4, 4, 2, 2], [4, 4, 2, 2], [4, 4, 4, 2, 2]]
-            slice_shrink = [[32, 16, 8, 4], [32, 16, 8, 4], [64, 32, 16, 8, 4]]
-            if i == len(reg_types)-1:
-                slice_sigma = [[4, 4, 2, 2], [4, 4, 2, 2], [4, 4, 4, 2, 2, 1]]
-                slice_shrink = [[32, 16, 8, 4], [32, 16, 8, 4], [64, 32, 16, 8, 4, 2]]
-        else:
-            sys.exit('The registration type "' + str(iteration) + '" does not exist. Reconstruction stopped.')
-
-        # Define name and path for registered template to hv.
-        cur_reg_name = pretty(i) + '-' + template_type + '_reg_to_hv.nii'
-        cur_reg_path = cur_work_dir + cur_reg_name
-
-        if print_log:
-            print('Registering hv to given volume - either MR or atlas')
-
-        if not path.exists(cur_reg_path):
-            registration(template, hv, cur_work_dir + 'hv_to_template.nii', reg, vol_sigma,
-                         vol_shrink, params, dim=3, initial_trans=initial_transform,
-                         inverse_result=cur_reg_path)
-
-        # Define and create directory for sliced template.
-        template_slice_dir = cur_work_dir + cur_reg_name.split('.')[0] + '-slices/'
-        make_dirs(template_slice_dir)
-
-        if print_log:
-            print('Slicing', cur_reg_path, 'into', template_slice_dir)
-
-        # Slice template.
-        slice_volume(cur_reg_path, template_slice_dir)
-        template_slices = files_in_dir(template_slice_dir, '.nii')
-
-        if len(raw_slices) == len(template_slices):
-            print('Number of images matches.')
-        else:
-            sys.exit('Mismatch between number of images.')
-
-        # Make directory for slice to slice registrations.
-        reg_slices_dir = cur_work_dir + pretty(i) + '-raw_reg_template/'
-        make_dirs(reg_slices_dir)
-
-        # Registered slice postfix.
-        reg_slice_postfix = '_reg'
-
-        # Perform slice to slice registrations.
-        for n, (raw_slice, template_slice) in enumerate(zip(raw_slices, template_slices)):
-            reg_slice = reg_slices_dir + pretty(i) + '_S' + pretty(n) + reg_slice_postfix + '.nii'
-            if not path.exists(reg_slice):
-                registration(template_slice, raw_slice, reg_slice, reg,
-                             slice_sigma, slice_shrink, params, 2)
-
-        # Stack registered slices into new volume.
-        reg_slices = files_ind_dir_wo_string(reg_slices_dir, exclusion_string='-Transformation')
-        final_volume = out_dir + pretty(i + 1) + volume_postfix + '.nii'
-        make_volume(reg_slices, final_volume)
-
-    # Return final reconstructed brain volume.
-    return final_volume
-
-
-if __name__ == '__main__':
-    arguments = sys.argv
-    log = False
-
-    if len(arguments) == 1:
-        print(usage_message)
-    else:
-        try:
-            opts, args = getopt.getopt(arguments[1:],
-                                       "he:s:r:t:d:o:l:",
-                                       ["template=", "slices=", "rawvolume", "inittrans", "id", "output", "log="])
-        except getopt.GetoptError:
-            print('\nSomething is not right in your syntax. Tak a look at your possibilities:\n', usage_message)
-            sys.exit(2)
-        for opt, arg in opts:
-            if opt == "-h":
-                print(usage_message)
-                sys.exit()
-            elif opt in ("-e", "--template"):
-                input_template = arg
-            elif opt in ("-s", "--slices"):
-                input_slices = arg
-            elif opt in ("-r", "--rawvolume"):
-                input_raw_volume = arg
-            elif opt in ("-t", "--inittrans"):
-                input_init_trans = arg
-            elif opt in ("-d", "--id"):
-                input_brain_id = arg
-            elif opt in ("-o", "--output"):
-                output_dir = arg
-            elif opt in ("-l", "--log"):
-                if arg == 'True' or arg == 'true':
-                    input_print_log = True
-                elif arg == 'False' or arg == 'false':
-                    input_print_log = False
-                else:
-                    print('\nSomething is not right in your syntax. Tak a look at your possibilities:\n', usage_message)
-                    sys.exit(2)
-        try:
-            co_register(input_template,
-                        input_slices,
-                        input_raw_volume,
-                        input_init_trans,
-                        output_dir,
-                        input_brain_id,
-                        print_log=input_print_log)
-        except NameError:
-            print('\nYou are not (properly) defining all the input variables. Tak a look at what is needed:\n',
-                    usage_message)

data/template_creation_pipline/helper_functions.py

@@ -1,324 +0,0 @@
-import os
-import nibabel as nib
-import numpy as np
-from nipype.interfaces.ants.segmentation import N4BiasFieldCorrection
-import math
-from sklearn.metrics import mutual_info_score
-from copy import copy
-
-
-# Makes directory if not already present.
-def make_dirs(paths):
-    if isinstance(paths, str):
-        if not os.path.exists(paths):
-            os.makedirs(paths)
-            print(paths + ' directory made.')
-    elif isinstance(paths, list):
-        for path in paths:
-            if not os.path.exists(path):
-                os.makedirs(path)
-                print(path + ' directory made.')
-    else:
-        print('The given path was neither a single path or a list of paths')
-
-
-# Get list of paths of files in dir with specific file type.
-def files_in_dir(path, file_type=''):
-    path_list = []
-    for file in os.listdir(path):
-        if file.startswith('.'):
-            continue
-        if file.endswith(file_type):
-            path_list.append(path + file)
-    return sorted(path_list)
-
-
-def files_ind_dir_wo_string(path, exclusion_string):
-    path_list = []
-    for file in os.listdir(path):
-        if exclusion_string not in file:
-            path_list.append(path + file)
-    return sorted(path_list)
-
-
-# Get list of images with specific file type in folders. Also returns list of folders.
-def images_in_folders(path, file_type=''):
-    path_list = []
-    folder_list = []
-    for folder in os.listdir(path):
-        folder_list.append(path + folder)
-        for file in os.listdir(path + folder + '/'):
-            if file.startswith('.'):
-                continue
-            if file.endswith(file_tyoe):
-                path_list.append(path + folder + '/' + file)
-    return sorted(path_list), sorted(folder_list)  # sorted added by FFS 26/6/2019.
-
-
-# Make .nii image from raw data. Will give specific pixel dimensions, if given.
-def make_nii(data, pixdim=[]):
-    nii = nib.Nifti1Image(data, np.eye(4))
-    if len(pixdim) != 0:
-        hdr = nii.header
-        hdr['pixdim'] = pixdim
-        nii = nib.Nifti1Pair(data, None, hdr)
-    return nii
-
-
-# Make volume from list of consecutive images.
-def make_volume(images, output):
-    init_image = nib.load(images[0])
-    width, height = init_image.shape
-    depth = len(images)
-    init_pixdim = init_image.header['pixdim']
-    volume_array = np.zeros((width, depth, height))
-
-    for i, image in enumerate(images):
-        image_data = np.rot90(nib.load(image).get_data(), 2)
-        volume_array[:, depth - i - 1, :] = image_data
-    nib.save(make_nii(volume_array, [1., init_pixdim[1], 0.100, init_pixdim[2], 1., 1., 1., 1.]), output)
-
-
-# Get parent directory.
-def get_parent_dir(path):
-    folders = path.split('/')
-    if folders[-1] == '':
-        parent_dir = '/'.join(folders[0:-2]) + '/'
-    else:
-        parent_dir = '/'.join(folders[0:-1]) + '/'
-    return parent_dir
-
-
-# Reorient MR
-def reorient_mr(MR, result, inverted=False):
-    original = nib.load(MR)
-    original_data = original.get_data()
-    pixdim = original.header['pixdim']
-    if inverted:
-        reoriented_data = np.fliplr(np.rot90(original_data, 1, axes=(1, 2)))
-    else:
-        reoriented_data = np.rot90(original_data, 1, axes=(1, 2))
-    nib.save(make_nii(reoriented_data, [1., pixdim[1], pixdim[3], pixdim[2], 1., 1., 1., 1.]), result)
-
-
-# Reorient atlas
-def reorient_atlas(atlas, result):
-    original = nib.load(atlas)
-    original_data = original.get_data()
-    pixdim = original.header['pixdim']
-    print(original_data.shape)
-    reoriented_data = np.reshape(original_data, (456, 528, 320), order='A')
-    print(reoriented_data.shape)
-    nib.save(make_nii(reoriented_data, [1., pixdim[1], pixdim[3], pixdim[2], 1., 1., 1., 1.]), result)
-
-
-def slice_volume(volume, slices_path):
-    current_mr = nib.load(volume)
-    current_data = current_mr.get_data()
-    number_of_slices = current_data.shape[1]
-    pixdim = current_mr.header['pixdim']
-
-    for i in range(number_of_slices):
-        image_nr = i + 1
-        image_nr_beauty = '0' + str(image_nr) if image_nr < 10 else str(image_nr)
-        nib.save(make_nii(np.rot90(current_data[:, number_of_slices - image_nr, :], 2),
-                          [1., pixdim[1], pixdim[3], pixdim[2], 1., 1., 1., 1.]),
-                 slices_path + 'Template_S' + image_nr_beauty + '.nii')
-
-
-def bias_field_correction(image, output, dim=2):
-    N4BiasFieldCorrection(
-        input_image=image,
-        output_image=output,
-        dimension=dim
-    ).run()
-
-
-def flip_volume_lr(image, output):
-    nii = nib.load(image)
-    hdr = nii.header
-    data = nii.get_data()
-    flipped = np.flip(data, axis=0)
-    nib.save(nib.Nifti1Pair(flipped, None, hdr), output)
-
-
-def add_slash_to_path(path):
-    return path + '/' if not path[-1] == '/' else path
-
-
-def identity_transform(work_dir, dimensions):
-    path = work_dir + 'identity_transform.txt'
-    identity_transform = open(path, 'w+')
-    if dimensions == 3:
-        identity_transform.write('Transform: MatrixOffsetTransformBase_double_3_3'
-                                 '\nParameters: 1 0 0 0 1 0 0 0 1 0 0 0'
-                                 '\nFixedParameters: 0 0 0\n')
-    elif dimensions == 2:
-        identity_transform.write('Transform: MatrixOffsetTransformBase_double_2_2'
-                                 '\nParameters: 1 0 0 1 0 0'
-                                 '\nFixedParameters: 0 0\n')
-    else:
-        print('Helper function identity_transform does only support 2d and 3d.')
-    return path
-
-
-# Returns True i given path exists.
-def path_exists(path):
-    return True if os.path.exists(path) else False
-
-
-def robust_avg_volumes(volumes, output):
-    # Determine images matrix dimensions.
-    depth_list = []
-    init_nii = nib.load(volumes[0])
-    pix_dim = init_nii.header['pixdim']
-    a, b, c = init_nii.shape
-    d = len(volumes)
-
-    for volume_shape in volumes:
-        depth_list.append(nib.load(volume_shape).shape[1])
-
-    b = max(depth_list)
-
-    matrix_shape = (a, b, c, d)
-    images_matrix = np.zeros(matrix_shape)
-    print('Average:')
-    for i, volume in enumerate(volumes):
-        print('Adding', volume)
-        nii = nib.load(volume)
-        images_matrix[:, 0:nii.shape[1], :, i] = nii.get_data()
-
-    # Calculate mean and std.
-    images_mean = np.mean(images_matrix, axis=3)
-    images_std = images_matrix.std(3)
-    images_std_scaled = images_std * 2
-
-    # Duplicate mean and std to full dimensions for easy matching.
-    mean_matrix = np.repeat(images_mean[:, :, :, np.newaxis], d, axis=3)
-    scaled_std_matrix = np.repeat(images_std_scaled[:, :, :, np.newaxis], d, axis=3)
-
-    robust_matrix = np.where(
-        (images_matrix < mean_matrix - scaled_std_matrix) | (images_matrix > mean_matrix + scaled_std_matrix),
-        np.nan,
-        images_matrix)
-
-    robust_mean = np.nanmedian(robust_matrix, axis=3)
-    nib.save(make_nii(robust_mean, pixdim=pix_dim), output)
-
-
-def split_images(images, output_folder, max_depth=50):
-    temp_folder = output_folder
-    make_dirs(temp_folder)
-
-    for image in images:
-        # Extract name from path.
-        name = image.split('/')[-1].split('.')[0]
-
-        # Load nii, header and data.
-        nii = nib.load(image)
-        header = nii.header
-        data = nii.get_data()
-
-        # Get depth from header and calculate number of parts.
-        depth = nii.shape[1]
-        nr_of_parts = math.ceil(depth / max_depth)
-
-        # If first iteration create empty list of list of parts.
-        if image == images[0]:
-            parts = [[] for _ in range(nr_of_parts)]
-            output_shape = nii.shape
-            output_header = copy(header)
-
-        # Make parts.
-        for i in range(nr_of_parts):
-            # Define form and to index.
-            from_i = i * max_depth
-            to_i = min((i + 1) * max_depth - 1, depth)
-
-            # Calculate depth of current part. Set in nii header.
-            cur_depth = to_i - from_i + 1
-            header['dim'][2] = cur_depth
-
-            # Slice data.
-            cur_data = data[:, from_i:to_i+1, :]
-            cur_path = temp_folder + name + '_' + str(i) + '.nii'
-
-            # Append to parts list.
-            parts[i] += [cur_path]
-
-            # Save part.
-            nib.save(nib.Nifti1Image(cur_data, None, header), cur_path)
-
-    return parts, output_shape, output_header
-
-
-def robust_large_images(images, output, max_depth=50):
-    temp_folder = get_parent_dir(images[0]) + 'avg_temp_volumes/'
-    parts, output_shape, output_header = split_images(images, temp_folder, max_depth=max_depth)
-
-    avgs = []
-    for n, part in enumerate(parts):
-        print('Averaging: ' + str(part))
-        cur_avg = temp_folder + str(n) + '_avg.nii'
-        robust_avg_volumes(part, cur_avg)
-        avgs += [cur_avg]
-
-    final_avg = np.zeros(output_shape)
-
-    for m, avg in enumerate(avgs):
-        print('Appending avg: ' + str(avg))
-        final_avg[:, m*max_depth:min((m+1)*max_depth, output_shape[1]), :] = nib.load(avg).get_data()
-
-    nib.save(nib.Nifti1Image(final_avg, None, output_header), output)
-
-
-def avg_template(image, template, template_weight, output):
-    print(image, template)
-    a = nib.load(image)
-    b = nib.load(template)
-    pix_dim = a.header['pixdim']
-
-    depth = max([a.shape[1], b.shape[1]])
-
-    a_scaled = np.zeros((a.shape[0], depth, a.shape[2]))
-    b_scaled = np.zeros((a.shape[0], depth, a.shape[2]))
-
-    a_scaled[:a.shape[0], :a.shape[1], :a.shape[2]] = a.get_data() * (1 - template_weight)
-    b_scaled[:b.shape[0], :b.shape[1], :b.shape[2]] = b.get_data() * template_weight
-
-    weighted = (a_scaled + b_scaled) / 2
-    nib.save(make_nii(weighted, pixdim=pix_dim), output)
-
-
-def normal_list_1d(sigma, n, mu=0.5, max_is_one=False):
-    """
-    Normal distribution.
-    :param sigma: sigma.
-    :param n: number of values.
-    :param mu: [0.0:1.0] a float determining the relative position of the peak. 0.0 left, 1.0 right. Default: 0.5.
-    :return: [x] a list of the natural numbers from 0 to n. [y]: normal value
-    """
-    mu = mu*n - 0.5
-    x = [x for x in range(n)]
-    y = [1 / (math.sqrt(2 * math.pi * sigma ** 2)) * math.e ** ((-(xi - mu) ** 2) / (2 * sigma ** 2)) for xi in x]
-
-    if max_is_one:
-        max_y = max(y)
-        y_adjusted = [yi * 1/max_y for yi in y]
-        return x, y_adjusted
-
-    return x, y
-
-
-def mutual_information(x, y, bins):
-    c_xy = np.histogram2d(x.ravel(), y.ravel(), bins)[0]
-    mi = mutual_info_score(None, None, contingency=c_xy)
-    return mi
-
-
-def immediate_sub_dirs(directory, full_path=True):
-    sub_dirs = sorted(next(os.walk(directory))[1])
-    if full_path:
-        return [directory + x + '/' for x in sub_dirs]
-    else:
-        return sub_dirs

data/template_creation_pipline/main.py

@@ -1,105 +0,0 @@
-from pre_processing import pre_process
-from co_registration import co_register
-from template_creation import create_template_from_images
-import helper_functions as hf
-import sys
-import os
-from shutil import copyfile
-
-
-def pre_process_brain_tiffs(tif_path, result_path, mouse_id):
-    """
-    Utility function pre-processing all tiffs in folder, outputting pre-processed nifties in result_path.
-    :param tif_path: The tiffs directory.
-    :param result_path: The output directory.
-    :param mouse_id: The mouse id.
-    :return: Returns nothing.
-    """
-    # Get all images and the enclosing tiff folders.
-    images, folders = hf.images_in_folders(tif_path)
-
-    # Ensure that the number of images matches the naming of the folders
-    if len(images) == int(folders[-1][-2:]):
-        print('Number of images matches the naming on the folders')
-    else:
-        sys.exit('Mismatch between number of images and naming of last TIF folder.')
-
-    # Make output directory.
-    hf.make_dirs([result_path])
-
-    # Loop through each image and pre-process.
-    for image_nr, image_path in enumerate(images, 1):
-        image_nr_beauty = str(image_nr).zfill(2)
-        output = result_path + mouse_id + '_S' + image_nr_beauty + '.nii'
-        if not os.path.exists(output):
-            pre_process(image_path, output)
-
-
-def main(preprocess=False, coreg=False, make_template=False):
-    """
-    Main function performing pre-processing of tiffs, 3-dimensional reconstruction to reference and population-based
-    template creation.
-    :param preprocess: Boolean [Default: False]. Will perform pre-processing of tiffs if set to true.
-    :param coreg: Boolean [Default: False]. Will perform brain reconstruction based on the sequential slices if set to
-    true.
-    :param make_template: Boolean [Default: False]. Will perform population-based template creation based on the
-    reconstructed brain volumes if set to true.
-    :return: Returns nothing.
-    """
-    # Working directory
-    sandbox = '/run/media/frederik/S_T1/frederik_filip_mqr375/auto-seg/' # Should be removed for final submission.
-
-    # Get list of mice ids.
-    sandbox_mice = sandbox + 'mice/'
-    mice = os.listdir(sandbox_mice)
-
-    if coreg or preprocess:
-        # Loop through mice ids.
-        for mouse in sorted(mice):
-            # Current mouse directory and subdirectories.
-            work_dir = sandbox_mice + mouse + '/'
-            raw_dir = work_dir + 'TIF/'
-            pre_processed_dir = work_dir + 'preprocessed/'
-            out_dir = work_dir + 'volumes/'
-
-            # Make volumes directory.
-            hf.make_dirs(out_dir)
-
-            # Define raw volume and initial transform path (if present).
-            initial_vol_trans = out_dir + 'initial_transform_' + mouse + '.txt'
-            raw_vol_path = out_dir + '00_hv_raw_' + mouse + '.nii'
-
-            # Pre process tiff images.
-            if preprocess:
-                pre_process_brain_tiffs(raw_dir, pre_processed_dir, mouse)
-
-            # Create raw volume if not already present.
-            if not os.path.exists(raw_vol_path):
-                hf.make_volume(hf.files_in_dir(pre_processed_dir, '.nii'), raw_vol_path)
-
-            # Check that an initial transform is present.
-            if not os.path.exists(initial_vol_trans):
-                sys.exit('An initial moving transform was not found for at least brain with id : ' + mouse +
-                         '\nMake sure it is named and placed exactly: ' + initial_vol_trans)
-
-            # -- This call should be manually modified depending on first or second iteration -- #
-            # Perform co-registration to reference (Allen in first iteration. DAPI template in second iteration).
-            # Manually created initial transform is only used in the first iteration.
-            # First iteration had 0 non-linear steps (parameter nl=0).
-            if coreg:
-                mouse_coreg = co_register(sandbox + 'second-iteration_average_allen_0_15.nii.gz',
-                                          pre_processed_dir, raw_vol_path, '', out_dir, mouse,
-                                          print_log=True, nl=2)
-                copyfile(mouse_coreg, sandbox + 'volumes/' + mouse + '.nii')
-
-    # Population based template creation using symmetric modelling.
-    if make_template:
-        create_template_from_images(
-            sandbox + 'volumes/',
-            sandbox + 'template/',
-            symmetric=True
-        )
-
-
-if __name__ == '__main__':
-    main(preprocess=True, coreg=True, make_template=True)

data/template_creation_pipline/slice_to_template.py

@@ -1,227 +0,0 @@
-import argparse
-from nipype.interfaces import ants
-import nibabel as nib
-import numpy as np
-import os
-from skimage.filters import threshold_otsu
-import matplotlib.pyplot as plt
-from helper_functions import make_nii, make_dirs, normal_list_1d, mutual_information
-from pre_processing import pre_process
-from skimage import morphology
-from scipy.ndimage import gaussian_filter
-import sys
-
-
-def slice_template(template, template_seg, out_dir, pos, padding=5):
-    # Load template data and header.
-    template = nib.load(template)
-    template_data = template.get_data()
-    template_header = template.header
-
-    # Load template segmentation data.
-    template_seg = nib.load(template_seg)
-    template_seg_data = template_seg.get_data()
-
-    # Set pixel dimensions.
-    slice_pixdim = [1., template_header['pixdim'][1], template_header['pixdim'][3], 1., 1., 1., 1., 1.]
-
-    images = []
-
-    for i in range(pos-padding, pos+padding+1):
-        output = out_dir + 'S' + str(i).zfill(2) + '.nii'
-        seg_output = out_dir + 'S' + str(i).zfill(2) + '_seg.nii'
-        if not os.path.exists(output):
-            nib.save(make_nii(np.rot90(template_data[:, i, :], 2), pixdim=slice_pixdim), output)
-        if not os.path.exists(seg_output):
-            nib.save(make_nii(np.rot90(template_seg_data[:, i, :], 2), pixdim=slice_pixdim), seg_output)
-        images.append(output)
-
-    # Returns a list of image paths.
-    return images
-
-
-def reg_slice(fixed, moving, output):
-    registration = ants.Registration(
-        dimension=2,
-        transforms=['Affine'],
-        transform_parameters=[(0.15, ), ],
-        initial_moving_transform_com=1,
-        metric=['MI'],
-        interpolation='BSpline',
-        fixed_image=[fixed],
-        moving_image=[moving],
-        metric_weight=[1],
-        radius_or_number_of_bins=[32],
-        number_of_iterations=[[1000, 1000, 1000, 1000, 1000]],
-        smoothing_sigmas=[[4, 4, 4, 2, 2]],
-        shrink_factors=[[64, 32, 16, 8, 4]],
-        convergence_threshold=[1.e-7],
-        sampling_strategy=['Regular'],
-        use_histogram_matching=True,
-        winsorize_lower_quantile=0.05,
-        winsorize_upper_quantile=0.95,
-        output_warped_image=output,
-        output_transform_prefix=output.split('.')[0] + '-Transformation',
-        num_threads=12
-    )
-    registration.run()
-
-
-def full_reg_slice(fixed, moving, output, seg):
-    output_pre = output.split('.')[0]
-    registration = ants.Registration(
-        dimension=2,
-        transforms=['Rigid', 'Affine', 'SyN'],
-        transform_parameters=[(0.15,), (0.15, ), (0.3, 0.3, 0.5)],
-        metric=['MeanSquares']*3,
-        initial_moving_transform_com=1,
-        interpolation='BSpline',
-        fixed_image=[fixed],
-        moving_image=[moving],
-        metric_weight=[1]*3,
-        radius_or_number_of_bins=[32]*3,
-        number_of_iterations=[[1000, 1000, 1000, 1000], [1000, 1000, 1000, 1000, 1000], [1000, 1000, 1000, 1000, 1000, 1000]],
-        smoothing_sigmas=[[4, 4, 2, 2], [4, 4, 4, 2, 2], [4, 4, 4, 2, 2, 2]],
-        shrink_factors=[[32, 16, 8, 4], [64, 32, 16, 8, 4], [64, 32, 16, 8, 4, 2]],
-        convergence_threshold=[1.e-7],
-        sampling_strategy=['Regular']*3,
-        use_histogram_matching=True,
-        winsorize_lower_quantile=0.05,
-        winsorize_upper_quantile=0.95,
-        output_warped_image=output,
-        output_transform_prefix=output_pre + '-Transformation-',
-        num_threads=12
-    )
-    registration.run()
-
-    # Define transformation to be applied on segmentation. Dependent on the system running ANTs.
-    if os.path.exists(output_pre + '-Transformation-0GenericAffine.mat'):
-        transform_paths = [output_pre + '-Transformation-1Warp.nii.gz',
-                           output_pre + '-Transformation-0GenericAffine.mat']
-    elif os.path.exists(output_pre + '-Transformation-0DerivedInitialMovingTranslation.mat'):
-        transform_paths = [output_pre + '-Transformation-3Warp.nii.gz',
-                           output_pre + '-Transformation-2Affine.mat',
-                           output_pre + '-Transformation-1Rigid.mat',
-                           output_pre + '-Transformation-0DerivedInitialMovingTranslation.mat']
-    else:
-        sys.exit('The ANTs-output transformation names/postfix was not recognised.')
-
-    # Apply recognised transformations to slice segmentation.
-    ants.ApplyTransforms(
-        dimension=2,
-        input_image=seg,
-        reference_image=output,
-        transforms=transform_paths,
-        interpolation='NearestNeighbor',
-        output_image=output.split('.')[0] + '_seg.nii',
-        num_threads=12
-    ).run()
-
-
-def main(input_slice_path, template_path, template_seg_path, output_path,
-         bregma, preprocess_input, inverse_output, padding=6, return_result_dict=False):
-    # Directory variables.
-    work_dir = os.path.dirname(output_path) + '/'
-    temp_dir = work_dir + 'temp/'
-    template_slice_dir = temp_dir + 'template_slices/'
-    make_dirs([temp_dir, template_slice_dir])
-
-    # Inverse bool
-    output_inverse = inverse_output != None
-    print(output_inverse)
-
-    # Slice template
-    coor = int(9.7431 * bregma + 38.7044)
-    template_slices = slice_template(template_path, template_seg_path, template_slice_dir, coor, padding=padding)
-
-    # Pre-process input image.
-    if preprocess_input in ['true', 'True']:
-        pre_processed_output = input_slice_path.split('.')[0] + '.nii'
-        pre_process(input_slice_path, pre_processed_output, print_log=True, shrink=1)
-        slice_path = pre_processed_output
-    else:
-        slice_path = input_slice_path
-
-    # Rigid registration
-    registration_dir = temp_dir + 'rigid_registrations/'
-    registered_slices = []
-    make_dirs(registration_dir)
-    for template_slice in template_slices:
-        name = template_slice.split('/')[-1].split('.')[0]
-        registered_slice = registration_dir + name + '.nii'
-        reg_slice(slice_path, template_slice, registered_slice)
-        registered_slices.append(registered_slice)
-
-    # Similarity extrema
-    mutual_informations = []
-    input_nii = nib.load(slice_path)
-    input_data = input_nii.get_data()
-    input_data_smooth = gaussian_filter(input_data, 2)
-
-    for image in registered_slices:
-        nii = nib.load(image)
-        data = nii.get_data()
-        mutual_informations.append(mutual_information(input_data_smooth, data, 32)-0.5)
-
-    # Weighted similarity.
-    _, weights = normal_list_1d(10, len(mutual_informations), max_is_one=True)
-    weighted_mi = [mi*w for mi, w in zip(mutual_informations, weights)]
-    max_y = max(weighted_mi)
-    max_index = weighted_mi.index(max_y)
-    maximum_slice = template_slices[max_index]
-    maximum_slice_seg = maximum_slice.split('.')[0] + '_seg.nii'
-
-    # Save pdf of similarities, with minimum marked.
-    f = plt.figure()
-    range_min = coor - padding
-    range_max = coor + padding + 1
-    max_x = range_min + max_index
-    plt.plot(range(range_min, range_max), weighted_mi)
-    plt.plot(range(range_min, range_max), mutual_informations)
-    plt.plot(max_x, max_y, 'rx')
-    plt.plot(range_min + mutual_informations.index(max(mutual_informations)), max(mutual_informations), 'rx')
-    plt.title('Mutual information of template slice registered to input slice.')
-    f.savefig(temp_dir + 'mutual_information.pdf')
-
-    # Run full (rigid, affine and nonlinear) registration on selected template slice.
-    # Transform segmentation accordingly.
-    full_reg_slice(slice_path, maximum_slice, output_path, maximum_slice_seg)
-
-    # Load segmentation, threshold, and close potential holes.
-    slice_data = nib.load(slice_path).get_data()
-    slice_seg_path = output_path.split('.nii')[0] + '_seg.nii'
-    slice_seg_nii = nib.load(slice_seg_path)
-    slice_seg_header = slice_seg_nii.header
-    slice_seg_data = slice_seg_nii.get_data()
-
-    # Calculate threshold.
-    threshold = threshold_otsu(slice_data)/6
-    foreground = np.where(slice_data > threshold, 1, 0)
-    foreground = morphology.binary_closing(foreground, selem=morphology.square(4))
-    slice_seg_data[foreground == 0] = 0
-
-    # Save threshold segmentation.
-    nib.save(nib.Nifti1Image(slice_seg_data, None, slice_seg_header), output_path.split('.nii')[0] + '_seg_thres.nii')
-    if return_result_dict:
-        return {'image': slice_path,
-                'segmentation': slice_seg_path,
-                'threshold_segmentation': output_path.split('.nii')[0] + '_seg_thres.nii',
-                'index': max_x,
-                'bregma': round((max_x - 38.7044)/9.7431, 2)}
-
-if __name__ == '__main__':
-    parser = argparse.ArgumentParser(description='Register a 2 dimensional coronal brain slice to a template.')
-
-    parser.add_argument("-s", "--slice", required=True, metavar='\b', help="Existing path to the slice.")
-    parser.add_argument("-t", "--template", required=True, metavar='\b', help="Existing path to the template.")
-    parser.add_argument("-e", "--segmentation", required=True, metavar='\b', help="Existing path to the template segmentation.")
-    parser.add_argument("-o", "--output", required=True, metavar='\b', help="Output.")
-    parser.add_argument("-b", "--bregma", required=True, metavar='\b', help="Approximate bregma coordinate of slice.")
-    parser.add_argument("-p", "--preprocess", required=True, metavar='\b', help="""Pre-process input image. Includes
-    finding biggest connected component, normalizing and biasfieldcorrection.""")
-    parser.add_argument("-i", "--inverse", required=True, metavar='\b', help="""Path to input slice inversely
-    transformed to template space. Will not transform inversely if not set.""")
-
-    args = parser.parse_args()
-
-    main(args.slice, args.template, args.segmentation, args.output, args.bregma, args.preprocess, args.inverse)

data/template_creation_pipline/template_creation.py

@@ -1,266 +0,0 @@
-from helper_functions import files_in_dir, make_dirs, add_slash_to_path, flip_volume_lr, robust_avg_volumes,\
-    files_ind_dir_wo_string, avg_template
-import sys
-import getopt
-import os
-import shutil
-from nipype.interfaces import ants
-
-
-# Help message.
-tab = '     '
-usage_message = '\nCOMMAND:\n' + tab + sys.argv[0] + \
-                '\n\nOPTIONS:\n' + \
-                tab + '-h, help menu\n' + \
-                tab + '-i, --images, <list of input images> or path to directory with only .nii images.\n' + \
-                tab + '-o, --output, <path to output directory>\n' + \
-                tab + '-s, --symmetric, [bool] use symmetric modelling. Default: False. Currently only works for 3d images.\n' + \
-                tab + '-t, --template, <path to template> can include template. Will not be included in final average. \n' + \
-                tab + '-w, --weight, [Double] specifies the weight a given template will have in averages. From 0 to 1. Default: 1/(number of images + 1)\n' + \
-                tab + '-l, --log, [bool]>\n'
-
-
-# Registration
-def registration(fixed, moving, result, reg_type, sigma, shrink, param, dim, initial_trans='', inverse_result=''):
-    """
-    Registration function. Registeres the moving image to the fixed with given parameters. If initial transform is not
-    given then the images are aligned by center of mass. Outputs inversely transformed result if given.
-    :param fixed: Path to fixed image.
-    :param moving: Path to moving image.
-    :param result: Path to result image.
-    :param reg_type: Registration type. List containing, Rigid, Affine, and/or SyN.
-    :param sigma: List of sigma parameters. Ordered according to the reg_type parameter.
-    :param shrink: List of shrink parameters. Ordered according to the reg_type parameter.
-    :param param: List of transform parameters. Ordered according to the reg_type parameter.
-    :param dim: Dimension of images. Typically 2 or 3.
-    :param initial_trans: Optional path to initial moving transform. If not given, the images will be matched initially
-    by aligning the center of mass.
-    :param inverse_result: Optional path to the inversely transformed fixed image.
-    :return: Returns nothing.
-    """
-    # Extract number of registrations.
-    steps = len(reg_type)
-
-    # Add 1000 iterations for each step.
-    if steps == 1:
-        iterations = [[1000] * len(sigma[0])]
-    else:
-        iterations = []
-        for i, reg in enumerate(reg_type):
-            iteration = [1000] * len(sigma[i])
-            iterations.append(iteration)
-
-    # Create registration instance.
-    reg_instance = ants.Registration(dimension=dim,
-                                     transforms=reg_type,
-                                     transform_parameters=param,
-                                     metric=['MI']*steps,
-                                     interpolation='BSpline',
-                                     fixed_image=[fixed],
-                                     moving_image=[moving],
-                                     metric_weight=[1]*steps,
-                                     radius_or_number_of_bins=[32]*steps,
-                                     number_of_iterations=iterations,
-                                     smoothing_sigmas=sigma,
-                                     shrink_factors=shrink,
-                                     convergence_threshold=[1.e-7],
-                                     sampling_strategy=['Regular']*steps,
-                                     use_histogram_matching=True,
-                                     winsorize_lower_quantile=0.05,
-                                     winsorize_upper_quantile=0.95,
-                                     output_warped_image=result,
-                                     write_composite_transform=True,
-                                     output_transform_prefix=result.split('.')[0] + '-Transformation',
-                                     num_threads=12)
-    # Add initial moving transform if given, else match by center of mass.
-    if not initial_trans == '':
-        reg_instance.inputs.initial_moving_transform = initial_trans
-    else:
-        reg_instance.inputs.initial_moving_transform_com = 1
-
-    # Output reverse results if path is given.
-    if not inverse_result == '':
-        reg_instance.inputs.output_inverse_warped_image = inverse_result
-
-    # Run registration.
-    reg_instance.run()
-
-
-def create_template_from_images(images_input, results_dir, symmetric=False, template='', template_weight=1,
-                                r=3, a=2, nl=3, print_log=False):
-    """
-    Population based template creation from input images. If input template is given this will be part of the averaging
-    with weight template_weight. Will perform r rigid, a rigid and affine, and nl rigid, affine, and non-linear
-    iterations. Will flip all input images in the left-right plane if symmetric is set to true.
-    :param images_input: List of paths to images.
-    :param results_dir: Working directory.
-    :param symmetric: Boolean [Default: False]. Will flip all input images in the left-right plane if symmetric is set
-    to true.
-    :param template: str [Default: '']. Will include template in averaging if path is specified.
-    :param template_weight: float [Default: 1]. Will weight the template with given weight if set to different from 1.
-    Can be in the range of 0 < template_weight < 1.
-    :param r: Int [Default: 3]. Number of rigid registrations.
-    :param a: Int [Default: 2]. Number of rigid and affine registrations.
-    :param nl: Int [Default: 3]. Number of rigid, affine, and non-linear registrations.
-    :param print_log: Boolean [Default: False]. Will print log if true.
-    :return: Returns nothing.
-    """
-    out_dir = add_slash_to_path(results_dir)
-    make_dirs(out_dir)
-    template_present = not template == ''
-
-    if print_log:
-        print('Determining images-variable format and defining input images.')
-    if isinstance(images_input, list):
-        images = images_input
-    elif isinstance(images_input, str):
-        images_input = add_slash_to_path(images_input)
-        images = files_in_dir(images_input, '.nii')
-
-    all_images = images
-
-    if symmetric:
-        flipped_path = out_dir + 'raw_flipped/'
-        flipped_images = []
-        make_dirs(flipped_path)
-        if print_log:
-            print('Creating flipped images.')
-        for image in images:
-            flipped_image = flipped_path + image.split('/')[-1].split('.')[0] + '_flipped.nii'
-            if not os.path.exists(flipped_image):
-                flip_volume_lr(image, flipped_image)
-            else:
-                if print_log:
-                    print('Flipped image', flipped_image, 'already present.')
-            flipped_images.append(flipped_image)
-
-        all_images.extend(flipped_images)
-
-    if print_log:
-        print('Defining iterations.')
-    reg_types = ['Rigid'] * r + ['Affine'] * a + ['SyN'] * nl
-
-    if print_log:
-        print(reg_types)
-
-    for i, iteration in enumerate(reg_types):
-        iteration_nr = str(i).zfill(2)
-        cur_it_dir = out_dir + iteration_nr + '/'
-        make_dirs(cur_it_dir)
-
-        if i == 0:
-            robust_avg_volumes(all_images, out_dir + '00-average.nii')
-            if template_present:
-                avg_template(out_dir + '00-average.nii',
-                             template,
-                             template_weight,
-                             out_dir + '00-average-template.nii')
-        else:
-            shutil.rmtree(out_dir + str(i-1).zfill(2))
-
-        average = out_dir + iteration_nr + '-average-template.nii' if template_present else out_dir + iteration_nr + '-average.nii'
-
-        # Iteration specific registration parameters.
-        if iteration == 'Rigid':
-            reg = ['Rigid']
-            params = [(0.25,)]
-            sigma = [[4, 4, 2, 2]]
-            shrink = [[32, 16, 8, 4]]
-            if i == 0:
-                sigma = [[4, 4, 4]]
-                shrink = [[32, 16, 8]]
-        elif iteration == 'Affine':
-            reg = ['Rigid', 'Affine']
-            params = [(0.25,), (0.25,), ]
-            sigma = [[4, 4, 2, 2], [4, 4, 2, 2]]
-            shrink = [[32, 16, 8, 4], [32, 16, 8, 4]]
-        elif iteration == 'SyN':
-            reg = ['Rigid', 'Affine', 'SyN']
-            params = [(0.25,), (0.25,), (0.15, 3, 0.5)]
-            sigma = [[4, 4, 2, 2], [4, 4, 2, 2], [4, 4, 4, 2]]
-            shrink = [[32, 16, 8, 4], [32, 16, 8, 4], [64, 32, 16, 8]]
-            if i == len(reg_types) - 1:
-                sigma = [[4, 4, 2, 2], [4, 4, 2, 2], [4, 4, 4, 2, 1]]
-                shrink = [[32, 16, 8, 4], [32, 16, 8, 4], [64, 32, 16, 8, 4]]
-        else:
-            sys.exit('The registration type "' + str(iteration) + '" does not exist. Reconstruction stopped.')
-
-        if template_present:
-            all_images.extend([template])
-
-        for image in all_images:
-            reg_output = cur_it_dir + image.split('/')[-1]
-            registration(average, image, reg_output, reg, sigma, shrink, params, dim=3)
-
-        registered_images = files_ind_dir_wo_string(cur_it_dir, '-Transformation')
-        registered_images_wo_allen = []
-
-        for reg_image in registered_images:
-            if 'allen' in reg_image:
-                registered_allen = reg_image
-            else:
-                registered_images_wo_allen.append(reg_image)
-
-        robust_avg_volumes(registered_images_wo_allen, out_dir + str(i+1).zfill(2) + '-average.nii')
-        if template_present:
-            avg_template(out_dir + str(i+1).zfill(2) + '-average.nii',
-                         registered_allen,
-                         template_weight,
-                         out_dir + str(i + 1).zfill(2) + '-average-template.nii')
-
-
-if __name__ == '__main__':
-    arguments = sys.argv
-
-    input_symmetric = False
-    input_robust = True
-    input_template = ''
-    input_weight = 1
-    show_log = False
-
-    if len(arguments) == 1:
-        print(usage_message)
-    else:
-        try:
-            opts, args = getopt.getopt(arguments[1:], "hi:o:s:t:w:l:", ["images=", "output=", "symmetric=",
-                                                                        "template=", "weight=", "log="])
-        except getopt.GetoptError:
-            print('\nSomething is not right in your syntax. Tak a look at your possibilities:\n', usage_message)
-            sys.exit(2)
-        for opt, arg in opts:
-            if opt == "-h":
-                print(usage_message)
-                sys.exit()
-            elif opt in ("-i", "--images"):
-                input_images = arg
-            elif opt in ("-o", "--output"):
-                output_directory = arg
-            elif opt in ("-s", "--symmetric"):
-                if arg == 'True' or arg == 'true':
-                    input_symmetric = True
-                elif arg == 'False' or arg == 'false':
-                    input_symmetric = False
-                else:
-                    print('\nSomething is not right in your syntax. Tak a look at your possibilities:\n', usage_message)
-            elif opt in ("-t", "--template"):
-                input_template = arg
-            elif opt in ("-w", "--weight"):
-                input_weight = arg
-            elif opt in ("-l", "--log"):
-                if arg == 'True' or arg == 'true':
-                    show_log = True
-                elif arg == 'False' or arg == 'false':
-                    show_log = False
-                else:
-                    print('\nSomething is not right in your syntax. Tak a look at your possibilities:\n', usage_message)
-                    sys.exit(2)
-        try:
-            create_template_from_images(input_images,
-                                        output_directory,
-                                        symmetric=input_symmetric,
-                                        template=input_template,
-                                        template_weight=input_weight,
-                                        print_log=show_log)
-        except NameError:
-            print('\nYou are not (properly) defining all the input variables. Tak a look at what is needed:\n',
-                  usage_message)