gin commit from rMBP-15-Malthe.local

New files: 6

template_creation_pipline/README.md

@@ -0,0 +1,47 @@
+# 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.

template_creation_pipline/co_registration.py

@@ -0,0 +1,288 @@
+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)

template_creation_pipline/helper_functions.py

@@ -0,0 +1,324 @@
+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

template_creation_pipline/main.py

@@ -0,0 +1,105 @@
+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)

template_creation_pipline/slice_to_template.py

@@ -0,0 +1,227 @@
+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)

template_creation_pipline/template_creation.py

@@ -0,0 +1,266 @@
+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)