import os
import glob
import nibabel as nib
import numpy as np
import shutil
from nipype.interfaces.ants import N4BiasFieldCorrection
from sklearn.feature_extraction.image import extract_patches as sk_extract_patches
from sklearn.utils import shuffle
import scipy.misc

num_mod = 2

def get_filename(set_name, case_idx, input_name, loc):
    pattern = '{0}/{1}/{3}/subject-{2}-{3}.nii'
    return pattern.format(loc, set_name, case_idx, input_name)

def get_set_name(case_idx):
    return 'Training' if case_idx < 11 else 'Testing'

def read_data(case_idx, input_name, loc):
    set_name = get_set_name(case_idx)
    image_path = get_filename(set_name, case_idx, input_name, loc)
    print(image_path)
    return nib.load(image_path)

def read_vol(case_idx, input_name, dir):
    image_data = read_data(case_idx, input_name, dir)
    return image_data.get_data()

def correct_bias(in_file, out_file):
    correct = N4BiasFieldCorrection()
    correct.inputs.input_image = in_file
    correct.inputs.output_image = out_file
    done = correct.run()
    return done.outputs.output_image

def normalise(case_idx, input_name, in_dir, out_dir,copy=False):
	set_name = get_set_name(case_idx)
	image_in_path = get_filename(set_name, case_idx, input_name, in_dir)
	image_out_path = get_filename(set_name, case_idx, input_name, out_dir)
	if copy:
		shutil.copy(image_in_path, image_out_path)
	else:
		correct_bias(image_in_path, image_out_path)
	print(image_in_path + " done.")

"""
To extract patches from a 3D image
"""
def extract_patches(volume, patch_shape, extraction_step,datype='float32'):
    patch_h, patch_w, patch_d = patch_shape[0], patch_shape[1], patch_shape[2]
    stride_h, stride_w, stride_d = extraction_step[0], extraction_step[1], extraction_step[2]
    img_h, img_w, img_d = volume.shape[0],volume.shape[1],volume.shape[2]
    N_patches_h = (img_h-patch_h)//stride_h+1
    N_patches_w = (img_w-patch_w)//stride_w+1
    N_patches_d = (img_d-patch_d)//stride_d+1
    N_patches_img = N_patches_h * N_patches_w * N_patches_d
    raw_patch_martrix = np.zeros((N_patches_img,patch_h,patch_w,patch_d),dtype=datype)
    k=0

    #iterator over all the patches
    for h in range((img_h-patch_h)//stride_h+1):
        for w in range((img_w-patch_w)//stride_w+1):
            for d in range((img_d-patch_d)//stride_d+1):
                raw_patch_martrix[k]=volume[h*stride_h:(h*stride_h)+patch_h,\
                                			w*stride_w:(w*stride_w)+patch_w,\
                                  				d*stride_d:(d*stride_d)+patch_d]
                k+=1
    assert(k==N_patches_img)
    return raw_patch_martrix


"""
To extract labeled patches from array of 3D labeled images
"""
def get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
                    patch_shape,validating,testing,num_images_training):
    patch_shape_1d=patch_shape[0]
    # Extract patches from input volumes and ground truth
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod),dtype="float32")
    y = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d),dtype="uint8")
    for idx in range(len(T1_vols)) :
        y_length = len(y)
        if testing:
            print(("Extracting Patches from Image %2d ....")%(num_images_training+idx+2))
        elif validating:
            print(("Extracting Patches from Image %2d ....")%(num_images_training+idx+1))
        else:
            print(("Extracting Patches from Image %2d ....")%(1+idx))
        label_patches = extract_patches(label_vols[idx], patch_shape, extraction_step,
        														datype="uint8")

        # Select only those who are important for processing
        if testing or validating:
            valid_idxs = np.where(np.sum(label_patches, axis=(1, 2, 3)) != -1)
        else:
            valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3)) > 6000)

        # Filtering extracted patches
        label_patches = label_patches[valid_idxs]

        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
                                                patch_shape_1d, patch_shape_1d, num_mod),dtype="float32")))
        y = np.vstack((y, np.zeros((len(label_patches), patch_shape_1d,
                                                patch_shape_1d, patch_shape_1d),dtype="uint8")))

        y[y_length:, :, :, :] = label_patches

        # Sampling strategy: reject samples which labels are mostly 0 and have less than 6000 nonzero elements
        T1_train = extract_patches(T1_vols[idx], patch_shape, extraction_step,datype="float32")
        x[y_length:, :, :, :, 0] = T1_train[valid_idxs]

        # Sampling strategy: reject samples which labels are mostly 0 and have less than 6000 nonzero elements
        T2_train = extract_patches(T2_vols[idx], patch_shape, extraction_step,datype="float32")
        x[y_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x, y

"""
To preprocess the labeled training data
"""
def preprocess_dynamic_lab(dir, seed, num_classes, extraction_step,patch_shape,num_images_training=2,
                                validating=False,testing=False,num_images_testing=7):
    x = list(range(1,11))
    if testing:
        print("Testing")
        index_start = num_images_training + 2
        index_end = index_start + num_images_testing
        T1_vols = np.empty((num_images_testing, 144, 192, 256),dtype="float32")
        T2_vols = np.empty((num_images_testing, 144, 192, 256),dtype="float32")
        label_vols = np.empty((num_images_testing, 144, 192, 256),dtype="uint8")
    elif validating:
        print("Validating")
        index_start = num_images_training + 1
        index_end = index_start + 1
        T1_vols = np.empty((1, 144, 192, 256),dtype="float32")
        T2_vols = np.empty((1, 144, 192, 256),dtype="float32")
        label_vols = np.empty((1, 144, 192, 256),dtype="uint8")
    else:
        print("Training")
        index_start = 1
        index_end = index_start + num_images_training
        T1_vols = np.empty((num_images_training, 144, 192, 256),dtype="float32")
        T2_vols = np.empty((num_images_training, 144, 192, 256),dtype="float32")
        label_vols = np.empty((num_images_training, 144, 192, 256),dtype="uint8")
    i = 0
    for index in range(index_start, index_end) :
        print(x[index-1])
        T1_vols[i, :, :, :] = read_vol(x[index-1], 'T1', dir)
        T2_vols[i, :, :, :] = read_vol(x[index-1], 'T2', dir)
        label_vols[i, :, :, :] = read_vol(x[index-1], 'label', dir)
        i = i + 1
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std

    for i in range(T1_vols.shape[0]):
        T1_vols[i] = ((T1_vols[i] - np.min(T1_vols[i])) /
                                    (np.max(T1_vols[i])-np.min(T1_vols[i])))*255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = ((T2_vols[i] - np.min(T2_vols[i])) /
                                    (np.max(T2_vols[i])-np.min(T2_vols[i])))*255
    T1_vols = T1_vols/127.5 -1.
    T2_vols = T2_vols/127.5 -1.
    x,y=get_patches_lab(T1_vols,T2_vols,label_vols,extraction_step,patch_shape,validating=validating,
                                testing=testing,num_images_training=num_images_training)
    print("Total Extracted Labelled Patches Shape:",x.shape,y.shape)
    if testing:
        return np.rollaxis(x, 4, 1), label_vols
    elif validating:
        return np.rollaxis(x, 4, 1), y, label_vols
    else:
        return np.rollaxis(x, 4, 1), y

"""
To extract labeled patches from array of 3D ulabeled images
"""
def get_patches_unlab(T1_vols, T2_vols, extraction_step,patch_shape,dir):
    patch_shape_1d=patch_shape[0]
    # Extract patches from input volumes and ground truth
    label_ref= np.empty((1, 144, 192, 256),dtype="uint8")
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod))
    label_ref = read_vol(1, 'label', dir)
    for idx in range(len(T1_vols)) :

        x_length = len(x)
        print(("Processing the Image %2d ....")%(idx+11))
        label_patches = extract_patches(label_ref, patch_shape, extraction_step)

        # Select only those who are important for processing
        # Sampling strategy: reject samples which labels are mostly 0 and have less than 6000 nonzero elements
        valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3)) > 6000)

        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
                                            patch_shape_1d, patch_shape_1d, num_mod))))

        T1_train = extract_patches(T1_vols[idx], patch_shape, extraction_step,datype="float32")
        x[x_length:, :, :, :, 0] = T1_train[valid_idxs]

        T2_train = extract_patches(T2_vols[idx], patch_shape, extraction_step,datype="float32")
        x[x_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x

"""
To preprocess the unlabeled training data
"""
def preprocess_dynamic_unlab(dir,extraction_step,patch_shape,num_images_training_unlab):
    T1_vols = np.empty((num_images_training_unlab, 144, 192, 256),dtype="float32")
    T2_vols = np.empty((num_images_training_unlab, 144, 192, 256),dtype="float32")
    for case_idx in range(11, 11+num_images_training_unlab) :
        T1_vols[(case_idx - 11), :, :, :] = read_vol(case_idx, 'T1', dir)
        T2_vols[(case_idx - 11), :, :, :] = read_vol(case_idx, 'T2', dir)
        #print(read_vol(case_idx, 'T2', dir).shape)
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = ((T1_vols[i] - np.min(T1_vols[i])) /
                                        (np.max(T1_vols[i])-np.min(T1_vols[i])))*255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = ((T2_vols[i] - np.min(T2_vols[i])) /
                                        (np.max(T2_vols[i])-np.min(T2_vols[i])))*255
    T1_vols = T1_vols/127.5 -1.
    T2_vols = T2_vols/127.5 -1.
    x=get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape,dir)
    print("Total Extracted Unlabeled Patches Shape:",x.shape)
    return np.rollaxis(x, 4, 1)

b44668422a047f675a229a8eb585daf33e970592
<mask token>


def get_filename(set_name, case_idx, input_name, loc):
    pattern = '{0}/{1}/{3}/subject-{2}-{3}.nii'
    return pattern.format(loc, set_name, case_idx, input_name)


def get_set_name(case_idx):
    return 'Training' if case_idx < 11 else 'Testing'


def read_data(case_idx, input_name, loc):
    set_name = get_set_name(case_idx)
    image_path = get_filename(set_name, case_idx, input_name, loc)
    print(image_path)
    return nib.load(image_path)


def read_vol(case_idx, input_name, dir):
    image_data = read_data(case_idx, input_name, dir)
    return image_data.get_data()


def correct_bias(in_file, out_file):
    correct = N4BiasFieldCorrection()
    correct.inputs.input_image = in_file
    correct.inputs.output_image = out_file
    done = correct.run()
    return done.outputs.output_image


<mask token>


def get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
    patch_shape, validating, testing, num_images_training):
    patch_shape_1d = patch_shape[0]
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d,
        num_mod), dtype='float32')
    y = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d), dtype
        ='uint8')
    for idx in range(len(T1_vols)):
        y_length = len(y)
        if testing:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 2))
        elif validating:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 1))
        else:
            print('Extracting Patches from Image %2d ....' % (1 + idx))
        label_patches = extract_patches(label_vols[idx], patch_shape,
            extraction_step, datype='uint8')
        if testing or validating:
            valid_idxs = np.where(np.sum(label_patches, axis=(1, 2, 3)) != -1)
        else:
            valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 
                2, 3)) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod), dtype='float32')))
        y = np.vstack((y, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d), dtype='uint8')))
        y[y_length:, :, :, :] = label_patches
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x, y


<mask token>


def preprocess_dynamic_lab(dir, seed, num_classes, extraction_step,
    patch_shape, num_images_training=2, validating=False, testing=False,
    num_images_testing=7):
    x = list(range(1, 11))
    if testing:
        print('Testing')
        index_start = num_images_training + 2
        index_end = index_start + num_images_testing
        T1_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        T2_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        label_vols = np.empty((num_images_testing, 144, 192, 256), dtype=
            'uint8')
    elif validating:
        print('Validating')
        index_start = num_images_training + 1
        index_end = index_start + 1
        T1_vols = np.empty((1, 144, 192, 256), dtype='float32')
        T2_vols = np.empty((1, 144, 192, 256), dtype='float32')
        label_vols = np.empty((1, 144, 192, 256), dtype='uint8')
    else:
        print('Training')
        index_start = 1
        index_end = index_start + num_images_training
        T1_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        T2_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        label_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'uint8')
    i = 0
    for index in range(index_start, index_end):
        print(x[index - 1])
        T1_vols[i, :, :, :] = read_vol(x[index - 1], 'T1', dir)
        T2_vols[i, :, :, :] = read_vol(x[index - 1], 'T2', dir)
        label_vols[i, :, :, :] = read_vol(x[index - 1], 'label', dir)
        i = i + 1
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = (T1_vols[i] - np.min(T1_vols[i])) / (np.max(T1_vols[i]
            ) - np.min(T1_vols[i])) * 255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = (T2_vols[i] - np.min(T2_vols[i])) / (np.max(T2_vols[i]
            ) - np.min(T2_vols[i])) * 255
    T1_vols = T1_vols / 127.5 - 1.0
    T2_vols = T2_vols / 127.5 - 1.0
    x, y = get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
        patch_shape, validating=validating, testing=testing,
        num_images_training=num_images_training)
    print('Total Extracted Labelled Patches Shape:', x.shape, y.shape)
    if testing:
        return np.rollaxis(x, 4, 1), label_vols
    elif validating:
        return np.rollaxis(x, 4, 1), y, label_vols
    else:
        return np.rollaxis(x, 4, 1), y


<mask token>


def get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape, dir):
    patch_shape_1d = patch_shape[0]
    label_ref = np.empty((1, 144, 192, 256), dtype='uint8')
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod))
    label_ref = read_vol(1, 'label', dir)
    for idx in range(len(T1_vols)):
        x_length = len(x)
        print('Processing the Image %2d ....' % (idx + 11))
        label_patches = extract_patches(label_ref, patch_shape, extraction_step
            )
        valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3
            )) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod))))
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x


<mask token>

<mask token>


def get_filename(set_name, case_idx, input_name, loc):
    pattern = '{0}/{1}/{3}/subject-{2}-{3}.nii'
    return pattern.format(loc, set_name, case_idx, input_name)


def get_set_name(case_idx):
    return 'Training' if case_idx < 11 else 'Testing'


def read_data(case_idx, input_name, loc):
    set_name = get_set_name(case_idx)
    image_path = get_filename(set_name, case_idx, input_name, loc)
    print(image_path)
    return nib.load(image_path)


def read_vol(case_idx, input_name, dir):
    image_data = read_data(case_idx, input_name, dir)
    return image_data.get_data()


def correct_bias(in_file, out_file):
    correct = N4BiasFieldCorrection()
    correct.inputs.input_image = in_file
    correct.inputs.output_image = out_file
    done = correct.run()
    return done.outputs.output_image


def normalise(case_idx, input_name, in_dir, out_dir, copy=False):
    set_name = get_set_name(case_idx)
    image_in_path = get_filename(set_name, case_idx, input_name, in_dir)
    image_out_path = get_filename(set_name, case_idx, input_name, out_dir)
    if copy:
        shutil.copy(image_in_path, image_out_path)
    else:
        correct_bias(image_in_path, image_out_path)
    print(image_in_path + ' done.')


<mask token>


def extract_patches(volume, patch_shape, extraction_step, datype='float32'):
    patch_h, patch_w, patch_d = patch_shape[0], patch_shape[1], patch_shape[2]
    stride_h, stride_w, stride_d = extraction_step[0], extraction_step[1
        ], extraction_step[2]
    img_h, img_w, img_d = volume.shape[0], volume.shape[1], volume.shape[2]
    N_patches_h = (img_h - patch_h) // stride_h + 1
    N_patches_w = (img_w - patch_w) // stride_w + 1
    N_patches_d = (img_d - patch_d) // stride_d + 1
    N_patches_img = N_patches_h * N_patches_w * N_patches_d
    raw_patch_martrix = np.zeros((N_patches_img, patch_h, patch_w, patch_d),
        dtype=datype)
    k = 0
    for h in range((img_h - patch_h) // stride_h + 1):
        for w in range((img_w - patch_w) // stride_w + 1):
            for d in range((img_d - patch_d) // stride_d + 1):
                raw_patch_martrix[k] = volume[h * stride_h:h * stride_h +
                    patch_h, w * stride_w:w * stride_w + patch_w, d *
                    stride_d:d * stride_d + patch_d]
                k += 1
    assert k == N_patches_img
    return raw_patch_martrix


<mask token>


def get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
    patch_shape, validating, testing, num_images_training):
    patch_shape_1d = patch_shape[0]
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d,
        num_mod), dtype='float32')
    y = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d), dtype
        ='uint8')
    for idx in range(len(T1_vols)):
        y_length = len(y)
        if testing:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 2))
        elif validating:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 1))
        else:
            print('Extracting Patches from Image %2d ....' % (1 + idx))
        label_patches = extract_patches(label_vols[idx], patch_shape,
            extraction_step, datype='uint8')
        if testing or validating:
            valid_idxs = np.where(np.sum(label_patches, axis=(1, 2, 3)) != -1)
        else:
            valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 
                2, 3)) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod), dtype='float32')))
        y = np.vstack((y, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d), dtype='uint8')))
        y[y_length:, :, :, :] = label_patches
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x, y


<mask token>


def preprocess_dynamic_lab(dir, seed, num_classes, extraction_step,
    patch_shape, num_images_training=2, validating=False, testing=False,
    num_images_testing=7):
    x = list(range(1, 11))
    if testing:
        print('Testing')
        index_start = num_images_training + 2
        index_end = index_start + num_images_testing
        T1_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        T2_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        label_vols = np.empty((num_images_testing, 144, 192, 256), dtype=
            'uint8')
    elif validating:
        print('Validating')
        index_start = num_images_training + 1
        index_end = index_start + 1
        T1_vols = np.empty((1, 144, 192, 256), dtype='float32')
        T2_vols = np.empty((1, 144, 192, 256), dtype='float32')
        label_vols = np.empty((1, 144, 192, 256), dtype='uint8')
    else:
        print('Training')
        index_start = 1
        index_end = index_start + num_images_training
        T1_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        T2_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        label_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'uint8')
    i = 0
    for index in range(index_start, index_end):
        print(x[index - 1])
        T1_vols[i, :, :, :] = read_vol(x[index - 1], 'T1', dir)
        T2_vols[i, :, :, :] = read_vol(x[index - 1], 'T2', dir)
        label_vols[i, :, :, :] = read_vol(x[index - 1], 'label', dir)
        i = i + 1
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = (T1_vols[i] - np.min(T1_vols[i])) / (np.max(T1_vols[i]
            ) - np.min(T1_vols[i])) * 255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = (T2_vols[i] - np.min(T2_vols[i])) / (np.max(T2_vols[i]
            ) - np.min(T2_vols[i])) * 255
    T1_vols = T1_vols / 127.5 - 1.0
    T2_vols = T2_vols / 127.5 - 1.0
    x, y = get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
        patch_shape, validating=validating, testing=testing,
        num_images_training=num_images_training)
    print('Total Extracted Labelled Patches Shape:', x.shape, y.shape)
    if testing:
        return np.rollaxis(x, 4, 1), label_vols
    elif validating:
        return np.rollaxis(x, 4, 1), y, label_vols
    else:
        return np.rollaxis(x, 4, 1), y


<mask token>


def get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape, dir):
    patch_shape_1d = patch_shape[0]
    label_ref = np.empty((1, 144, 192, 256), dtype='uint8')
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod))
    label_ref = read_vol(1, 'label', dir)
    for idx in range(len(T1_vols)):
        x_length = len(x)
        print('Processing the Image %2d ....' % (idx + 11))
        label_patches = extract_patches(label_ref, patch_shape, extraction_step
            )
        valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3
            )) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod))))
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x


<mask token>


def preprocess_dynamic_unlab(dir, extraction_step, patch_shape,
    num_images_training_unlab):
    T1_vols = np.empty((num_images_training_unlab, 144, 192, 256), dtype=
        'float32')
    T2_vols = np.empty((num_images_training_unlab, 144, 192, 256), dtype=
        'float32')
    for case_idx in range(11, 11 + num_images_training_unlab):
        T1_vols[case_idx - 11, :, :, :] = read_vol(case_idx, 'T1', dir)
        T2_vols[case_idx - 11, :, :, :] = read_vol(case_idx, 'T2', dir)
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = (T1_vols[i] - np.min(T1_vols[i])) / (np.max(T1_vols[i]
            ) - np.min(T1_vols[i])) * 255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = (T2_vols[i] - np.min(T2_vols[i])) / (np.max(T2_vols[i]
            ) - np.min(T2_vols[i])) * 255
    T1_vols = T1_vols / 127.5 - 1.0
    T2_vols = T2_vols / 127.5 - 1.0
    x = get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape, dir)
    print('Total Extracted Unlabeled Patches Shape:', x.shape)
    return np.rollaxis(x, 4, 1)

<mask token>
num_mod = 2


def get_filename(set_name, case_idx, input_name, loc):
    pattern = '{0}/{1}/{3}/subject-{2}-{3}.nii'
    return pattern.format(loc, set_name, case_idx, input_name)


def get_set_name(case_idx):
    return 'Training' if case_idx < 11 else 'Testing'


def read_data(case_idx, input_name, loc):
    set_name = get_set_name(case_idx)
    image_path = get_filename(set_name, case_idx, input_name, loc)
    print(image_path)
    return nib.load(image_path)


def read_vol(case_idx, input_name, dir):
    image_data = read_data(case_idx, input_name, dir)
    return image_data.get_data()


def correct_bias(in_file, out_file):
    correct = N4BiasFieldCorrection()
    correct.inputs.input_image = in_file
    correct.inputs.output_image = out_file
    done = correct.run()
    return done.outputs.output_image


def normalise(case_idx, input_name, in_dir, out_dir, copy=False):
    set_name = get_set_name(case_idx)
    image_in_path = get_filename(set_name, case_idx, input_name, in_dir)
    image_out_path = get_filename(set_name, case_idx, input_name, out_dir)
    if copy:
        shutil.copy(image_in_path, image_out_path)
    else:
        correct_bias(image_in_path, image_out_path)
    print(image_in_path + ' done.')


<mask token>


def extract_patches(volume, patch_shape, extraction_step, datype='float32'):
    patch_h, patch_w, patch_d = patch_shape[0], patch_shape[1], patch_shape[2]
    stride_h, stride_w, stride_d = extraction_step[0], extraction_step[1
        ], extraction_step[2]
    img_h, img_w, img_d = volume.shape[0], volume.shape[1], volume.shape[2]
    N_patches_h = (img_h - patch_h) // stride_h + 1
    N_patches_w = (img_w - patch_w) // stride_w + 1
    N_patches_d = (img_d - patch_d) // stride_d + 1
    N_patches_img = N_patches_h * N_patches_w * N_patches_d
    raw_patch_martrix = np.zeros((N_patches_img, patch_h, patch_w, patch_d),
        dtype=datype)
    k = 0
    for h in range((img_h - patch_h) // stride_h + 1):
        for w in range((img_w - patch_w) // stride_w + 1):
            for d in range((img_d - patch_d) // stride_d + 1):
                raw_patch_martrix[k] = volume[h * stride_h:h * stride_h +
                    patch_h, w * stride_w:w * stride_w + patch_w, d *
                    stride_d:d * stride_d + patch_d]
                k += 1
    assert k == N_patches_img
    return raw_patch_martrix


<mask token>


def get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
    patch_shape, validating, testing, num_images_training):
    patch_shape_1d = patch_shape[0]
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d,
        num_mod), dtype='float32')
    y = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d), dtype
        ='uint8')
    for idx in range(len(T1_vols)):
        y_length = len(y)
        if testing:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 2))
        elif validating:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 1))
        else:
            print('Extracting Patches from Image %2d ....' % (1 + idx))
        label_patches = extract_patches(label_vols[idx], patch_shape,
            extraction_step, datype='uint8')
        if testing or validating:
            valid_idxs = np.where(np.sum(label_patches, axis=(1, 2, 3)) != -1)
        else:
            valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 
                2, 3)) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod), dtype='float32')))
        y = np.vstack((y, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d), dtype='uint8')))
        y[y_length:, :, :, :] = label_patches
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x, y


<mask token>


def preprocess_dynamic_lab(dir, seed, num_classes, extraction_step,
    patch_shape, num_images_training=2, validating=False, testing=False,
    num_images_testing=7):
    x = list(range(1, 11))
    if testing:
        print('Testing')
        index_start = num_images_training + 2
        index_end = index_start + num_images_testing
        T1_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        T2_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        label_vols = np.empty((num_images_testing, 144, 192, 256), dtype=
            'uint8')
    elif validating:
        print('Validating')
        index_start = num_images_training + 1
        index_end = index_start + 1
        T1_vols = np.empty((1, 144, 192, 256), dtype='float32')
        T2_vols = np.empty((1, 144, 192, 256), dtype='float32')
        label_vols = np.empty((1, 144, 192, 256), dtype='uint8')
    else:
        print('Training')
        index_start = 1
        index_end = index_start + num_images_training
        T1_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        T2_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        label_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'uint8')
    i = 0
    for index in range(index_start, index_end):
        print(x[index - 1])
        T1_vols[i, :, :, :] = read_vol(x[index - 1], 'T1', dir)
        T2_vols[i, :, :, :] = read_vol(x[index - 1], 'T2', dir)
        label_vols[i, :, :, :] = read_vol(x[index - 1], 'label', dir)
        i = i + 1
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = (T1_vols[i] - np.min(T1_vols[i])) / (np.max(T1_vols[i]
            ) - np.min(T1_vols[i])) * 255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = (T2_vols[i] - np.min(T2_vols[i])) / (np.max(T2_vols[i]
            ) - np.min(T2_vols[i])) * 255
    T1_vols = T1_vols / 127.5 - 1.0
    T2_vols = T2_vols / 127.5 - 1.0
    x, y = get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
        patch_shape, validating=validating, testing=testing,
        num_images_training=num_images_training)
    print('Total Extracted Labelled Patches Shape:', x.shape, y.shape)
    if testing:
        return np.rollaxis(x, 4, 1), label_vols
    elif validating:
        return np.rollaxis(x, 4, 1), y, label_vols
    else:
        return np.rollaxis(x, 4, 1), y


<mask token>


def get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape, dir):
    patch_shape_1d = patch_shape[0]
    label_ref = np.empty((1, 144, 192, 256), dtype='uint8')
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod))
    label_ref = read_vol(1, 'label', dir)
    for idx in range(len(T1_vols)):
        x_length = len(x)
        print('Processing the Image %2d ....' % (idx + 11))
        label_patches = extract_patches(label_ref, patch_shape, extraction_step
            )
        valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3
            )) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod))))
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x


<mask token>


def preprocess_dynamic_unlab(dir, extraction_step, patch_shape,
    num_images_training_unlab):
    T1_vols = np.empty((num_images_training_unlab, 144, 192, 256), dtype=
        'float32')
    T2_vols = np.empty((num_images_training_unlab, 144, 192, 256), dtype=
        'float32')
    for case_idx in range(11, 11 + num_images_training_unlab):
        T1_vols[case_idx - 11, :, :, :] = read_vol(case_idx, 'T1', dir)
        T2_vols[case_idx - 11, :, :, :] = read_vol(case_idx, 'T2', dir)
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = (T1_vols[i] - np.min(T1_vols[i])) / (np.max(T1_vols[i]
            ) - np.min(T1_vols[i])) * 255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = (T2_vols[i] - np.min(T2_vols[i])) / (np.max(T2_vols[i]
            ) - np.min(T2_vols[i])) * 255
    T1_vols = T1_vols / 127.5 - 1.0
    T2_vols = T2_vols / 127.5 - 1.0
    x = get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape, dir)
    print('Total Extracted Unlabeled Patches Shape:', x.shape)
    return np.rollaxis(x, 4, 1)

import os
import glob
import nibabel as nib
import numpy as np
import shutil
from nipype.interfaces.ants import N4BiasFieldCorrection
from sklearn.feature_extraction.image import extract_patches as sk_extract_patches
from sklearn.utils import shuffle
import scipy.misc
num_mod = 2


def get_filename(set_name, case_idx, input_name, loc):
    pattern = '{0}/{1}/{3}/subject-{2}-{3}.nii'
    return pattern.format(loc, set_name, case_idx, input_name)


def get_set_name(case_idx):
    return 'Training' if case_idx < 11 else 'Testing'


def read_data(case_idx, input_name, loc):
    set_name = get_set_name(case_idx)
    image_path = get_filename(set_name, case_idx, input_name, loc)
    print(image_path)
    return nib.load(image_path)


def read_vol(case_idx, input_name, dir):
    image_data = read_data(case_idx, input_name, dir)
    return image_data.get_data()


def correct_bias(in_file, out_file):
    correct = N4BiasFieldCorrection()
    correct.inputs.input_image = in_file
    correct.inputs.output_image = out_file
    done = correct.run()
    return done.outputs.output_image


def normalise(case_idx, input_name, in_dir, out_dir, copy=False):
    set_name = get_set_name(case_idx)
    image_in_path = get_filename(set_name, case_idx, input_name, in_dir)
    image_out_path = get_filename(set_name, case_idx, input_name, out_dir)
    if copy:
        shutil.copy(image_in_path, image_out_path)
    else:
        correct_bias(image_in_path, image_out_path)
    print(image_in_path + ' done.')


<mask token>


def extract_patches(volume, patch_shape, extraction_step, datype='float32'):
    patch_h, patch_w, patch_d = patch_shape[0], patch_shape[1], patch_shape[2]
    stride_h, stride_w, stride_d = extraction_step[0], extraction_step[1
        ], extraction_step[2]
    img_h, img_w, img_d = volume.shape[0], volume.shape[1], volume.shape[2]
    N_patches_h = (img_h - patch_h) // stride_h + 1
    N_patches_w = (img_w - patch_w) // stride_w + 1
    N_patches_d = (img_d - patch_d) // stride_d + 1
    N_patches_img = N_patches_h * N_patches_w * N_patches_d
    raw_patch_martrix = np.zeros((N_patches_img, patch_h, patch_w, patch_d),
        dtype=datype)
    k = 0
    for h in range((img_h - patch_h) // stride_h + 1):
        for w in range((img_w - patch_w) // stride_w + 1):
            for d in range((img_d - patch_d) // stride_d + 1):
                raw_patch_martrix[k] = volume[h * stride_h:h * stride_h +
                    patch_h, w * stride_w:w * stride_w + patch_w, d *
                    stride_d:d * stride_d + patch_d]
                k += 1
    assert k == N_patches_img
    return raw_patch_martrix


<mask token>


def get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
    patch_shape, validating, testing, num_images_training):
    patch_shape_1d = patch_shape[0]
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d,
        num_mod), dtype='float32')
    y = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d), dtype
        ='uint8')
    for idx in range(len(T1_vols)):
        y_length = len(y)
        if testing:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 2))
        elif validating:
            print('Extracting Patches from Image %2d ....' % (
                num_images_training + idx + 1))
        else:
            print('Extracting Patches from Image %2d ....' % (1 + idx))
        label_patches = extract_patches(label_vols[idx], patch_shape,
            extraction_step, datype='uint8')
        if testing or validating:
            valid_idxs = np.where(np.sum(label_patches, axis=(1, 2, 3)) != -1)
        else:
            valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 
                2, 3)) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod), dtype='float32')))
        y = np.vstack((y, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d), dtype='uint8')))
        y[y_length:, :, :, :] = label_patches
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[y_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x, y


<mask token>


def preprocess_dynamic_lab(dir, seed, num_classes, extraction_step,
    patch_shape, num_images_training=2, validating=False, testing=False,
    num_images_testing=7):
    x = list(range(1, 11))
    if testing:
        print('Testing')
        index_start = num_images_training + 2
        index_end = index_start + num_images_testing
        T1_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        T2_vols = np.empty((num_images_testing, 144, 192, 256), dtype='float32'
            )
        label_vols = np.empty((num_images_testing, 144, 192, 256), dtype=
            'uint8')
    elif validating:
        print('Validating')
        index_start = num_images_training + 1
        index_end = index_start + 1
        T1_vols = np.empty((1, 144, 192, 256), dtype='float32')
        T2_vols = np.empty((1, 144, 192, 256), dtype='float32')
        label_vols = np.empty((1, 144, 192, 256), dtype='uint8')
    else:
        print('Training')
        index_start = 1
        index_end = index_start + num_images_training
        T1_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        T2_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'float32')
        label_vols = np.empty((num_images_training, 144, 192, 256), dtype=
            'uint8')
    i = 0
    for index in range(index_start, index_end):
        print(x[index - 1])
        T1_vols[i, :, :, :] = read_vol(x[index - 1], 'T1', dir)
        T2_vols[i, :, :, :] = read_vol(x[index - 1], 'T2', dir)
        label_vols[i, :, :, :] = read_vol(x[index - 1], 'label', dir)
        i = i + 1
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = (T1_vols[i] - np.min(T1_vols[i])) / (np.max(T1_vols[i]
            ) - np.min(T1_vols[i])) * 255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = (T2_vols[i] - np.min(T2_vols[i])) / (np.max(T2_vols[i]
            ) - np.min(T2_vols[i])) * 255
    T1_vols = T1_vols / 127.5 - 1.0
    T2_vols = T2_vols / 127.5 - 1.0
    x, y = get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
        patch_shape, validating=validating, testing=testing,
        num_images_training=num_images_training)
    print('Total Extracted Labelled Patches Shape:', x.shape, y.shape)
    if testing:
        return np.rollaxis(x, 4, 1), label_vols
    elif validating:
        return np.rollaxis(x, 4, 1), y, label_vols
    else:
        return np.rollaxis(x, 4, 1), y


<mask token>


def get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape, dir):
    patch_shape_1d = patch_shape[0]
    label_ref = np.empty((1, 144, 192, 256), dtype='uint8')
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod))
    label_ref = read_vol(1, 'label', dir)
    for idx in range(len(T1_vols)):
        x_length = len(x)
        print('Processing the Image %2d ....' % (idx + 11))
        label_patches = extract_patches(label_ref, patch_shape, extraction_step
            )
        valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3
            )) > 6000)
        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
            patch_shape_1d, patch_shape_1d, num_mod))))
        T1_train = extract_patches(T1_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 0] = T1_train[valid_idxs]
        T2_train = extract_patches(T2_vols[idx], patch_shape,
            extraction_step, datype='float32')
        x[x_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x


<mask token>


def preprocess_dynamic_unlab(dir, extraction_step, patch_shape,
    num_images_training_unlab):
    T1_vols = np.empty((num_images_training_unlab, 144, 192, 256), dtype=
        'float32')
    T2_vols = np.empty((num_images_training_unlab, 144, 192, 256), dtype=
        'float32')
    for case_idx in range(11, 11 + num_images_training_unlab):
        T1_vols[case_idx - 11, :, :, :] = read_vol(case_idx, 'T1', dir)
        T2_vols[case_idx - 11, :, :, :] = read_vol(case_idx, 'T2', dir)
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = (T1_vols[i] - np.min(T1_vols[i])) / (np.max(T1_vols[i]
            ) - np.min(T1_vols[i])) * 255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = (T2_vols[i] - np.min(T2_vols[i])) / (np.max(T2_vols[i]
            ) - np.min(T2_vols[i])) * 255
    T1_vols = T1_vols / 127.5 - 1.0
    T2_vols = T2_vols / 127.5 - 1.0
    x = get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape, dir)
    print('Total Extracted Unlabeled Patches Shape:', x.shape)
    return np.rollaxis(x, 4, 1)

import os
import glob
import nibabel as nib
import numpy as np
import shutil
from nipype.interfaces.ants import N4BiasFieldCorrection
from sklearn.feature_extraction.image import extract_patches as sk_extract_patches
from sklearn.utils import shuffle
import scipy.misc

num_mod = 2

def get_filename(set_name, case_idx, input_name, loc):
    pattern = '{0}/{1}/{3}/subject-{2}-{3}.nii'
    return pattern.format(loc, set_name, case_idx, input_name)

def get_set_name(case_idx):
    return 'Training' if case_idx < 11 else 'Testing'

def read_data(case_idx, input_name, loc):
    set_name = get_set_name(case_idx)
    image_path = get_filename(set_name, case_idx, input_name, loc)
    print(image_path)
    return nib.load(image_path)

def read_vol(case_idx, input_name, dir):
    image_data = read_data(case_idx, input_name, dir)
    return image_data.get_data()

def correct_bias(in_file, out_file):
    correct = N4BiasFieldCorrection()
    correct.inputs.input_image = in_file
    correct.inputs.output_image = out_file
    done = correct.run()
    return done.outputs.output_image

def normalise(case_idx, input_name, in_dir, out_dir,copy=False):
	set_name = get_set_name(case_idx)
	image_in_path = get_filename(set_name, case_idx, input_name, in_dir)
	image_out_path = get_filename(set_name, case_idx, input_name, out_dir)
	if copy:
		shutil.copy(image_in_path, image_out_path)
	else:
		correct_bias(image_in_path, image_out_path)
	print(image_in_path + " done.")

"""
To extract patches from a 3D image
"""
def extract_patches(volume, patch_shape, extraction_step,datype='float32'):
    patch_h, patch_w, patch_d = patch_shape[0], patch_shape[1], patch_shape[2]
    stride_h, stride_w, stride_d = extraction_step[0], extraction_step[1], extraction_step[2]
    img_h, img_w, img_d = volume.shape[0],volume.shape[1],volume.shape[2]
    N_patches_h = (img_h-patch_h)//stride_h+1
    N_patches_w = (img_w-patch_w)//stride_w+1
    N_patches_d = (img_d-patch_d)//stride_d+1
    N_patches_img = N_patches_h * N_patches_w * N_patches_d
    raw_patch_martrix = np.zeros((N_patches_img,patch_h,patch_w,patch_d),dtype=datype)
    k=0

    #iterator over all the patches
    for h in range((img_h-patch_h)//stride_h+1):
        for w in range((img_w-patch_w)//stride_w+1):
            for d in range((img_d-patch_d)//stride_d+1):
                raw_patch_martrix[k]=volume[h*stride_h:(h*stride_h)+patch_h,\
                                			w*stride_w:(w*stride_w)+patch_w,\
                                  				d*stride_d:(d*stride_d)+patch_d]
                k+=1
    assert(k==N_patches_img)
    return raw_patch_martrix


"""
To extract labeled patches from array of 3D labeled images
"""
def get_patches_lab(T1_vols, T2_vols, label_vols, extraction_step,
                    patch_shape,validating,testing,num_images_training):
    patch_shape_1d=patch_shape[0]
    # Extract patches from input volumes and ground truth
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod),dtype="float32")
    y = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d),dtype="uint8")
    for idx in range(len(T1_vols)) :
        y_length = len(y)
        if testing:
            print(("Extracting Patches from Image %2d ....")%(num_images_training+idx+2))
        elif validating:
            print(("Extracting Patches from Image %2d ....")%(num_images_training+idx+1))
        else:
            print(("Extracting Patches from Image %2d ....")%(1+idx))
        label_patches = extract_patches(label_vols[idx], patch_shape, extraction_step,
        														datype="uint8")

        # Select only those who are important for processing
        if testing or validating:
            valid_idxs = np.where(np.sum(label_patches, axis=(1, 2, 3)) != -1)
        else:
            valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3)) > 6000)

        # Filtering extracted patches
        label_patches = label_patches[valid_idxs]

        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
                                                patch_shape_1d, patch_shape_1d, num_mod),dtype="float32")))
        y = np.vstack((y, np.zeros((len(label_patches), patch_shape_1d,
                                                patch_shape_1d, patch_shape_1d),dtype="uint8")))

        y[y_length:, :, :, :] = label_patches

        # Sampling strategy: reject samples which labels are mostly 0 and have less than 6000 nonzero elements
        T1_train = extract_patches(T1_vols[idx], patch_shape, extraction_step,datype="float32")
        x[y_length:, :, :, :, 0] = T1_train[valid_idxs]

        # Sampling strategy: reject samples which labels are mostly 0 and have less than 6000 nonzero elements
        T2_train = extract_patches(T2_vols[idx], patch_shape, extraction_step,datype="float32")
        x[y_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x, y

"""
To preprocess the labeled training data
"""
def preprocess_dynamic_lab(dir, seed, num_classes, extraction_step,patch_shape,num_images_training=2,
                                validating=False,testing=False,num_images_testing=7):
    x = list(range(1,11))
    if testing:
        print("Testing")
        index_start = num_images_training + 2
        index_end = index_start + num_images_testing
        T1_vols = np.empty((num_images_testing, 144, 192, 256),dtype="float32")
        T2_vols = np.empty((num_images_testing, 144, 192, 256),dtype="float32")
        label_vols = np.empty((num_images_testing, 144, 192, 256),dtype="uint8")
    elif validating:
        print("Validating")
        index_start = num_images_training + 1
        index_end = index_start + 1
        T1_vols = np.empty((1, 144, 192, 256),dtype="float32")
        T2_vols = np.empty((1, 144, 192, 256),dtype="float32")
        label_vols = np.empty((1, 144, 192, 256),dtype="uint8")
    else:
        print("Training")
        index_start = 1
        index_end = index_start + num_images_training
        T1_vols = np.empty((num_images_training, 144, 192, 256),dtype="float32")
        T2_vols = np.empty((num_images_training, 144, 192, 256),dtype="float32")
        label_vols = np.empty((num_images_training, 144, 192, 256),dtype="uint8")
    i = 0
    for index in range(index_start, index_end) :
        print(x[index-1])
        T1_vols[i, :, :, :] = read_vol(x[index-1], 'T1', dir)
        T2_vols[i, :, :, :] = read_vol(x[index-1], 'T2', dir)
        label_vols[i, :, :, :] = read_vol(x[index-1], 'label', dir)
        i = i + 1
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std

    for i in range(T1_vols.shape[0]):
        T1_vols[i] = ((T1_vols[i] - np.min(T1_vols[i])) /
                                    (np.max(T1_vols[i])-np.min(T1_vols[i])))*255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = ((T2_vols[i] - np.min(T2_vols[i])) /
                                    (np.max(T2_vols[i])-np.min(T2_vols[i])))*255
    T1_vols = T1_vols/127.5 -1.
    T2_vols = T2_vols/127.5 -1.
    x,y=get_patches_lab(T1_vols,T2_vols,label_vols,extraction_step,patch_shape,validating=validating,
                                testing=testing,num_images_training=num_images_training)
    print("Total Extracted Labelled Patches Shape:",x.shape,y.shape)
    if testing:
        return np.rollaxis(x, 4, 1), label_vols
    elif validating:
        return np.rollaxis(x, 4, 1), y, label_vols
    else:
        return np.rollaxis(x, 4, 1), y

"""
To extract labeled patches from array of 3D ulabeled images
"""
def get_patches_unlab(T1_vols, T2_vols, extraction_step,patch_shape,dir):
    patch_shape_1d=patch_shape[0]
    # Extract patches from input volumes and ground truth
    label_ref= np.empty((1, 144, 192, 256),dtype="uint8")
    x = np.zeros((0, patch_shape_1d, patch_shape_1d, patch_shape_1d, num_mod))
    label_ref = read_vol(1, 'label', dir)
    for idx in range(len(T1_vols)) :

        x_length = len(x)
        print(("Processing the Image %2d ....")%(idx+11))
        label_patches = extract_patches(label_ref, patch_shape, extraction_step)

        # Select only those who are important for processing
        # Sampling strategy: reject samples which labels are mostly 0 and have less than 6000 nonzero elements
        valid_idxs = np.where(np.count_nonzero(label_patches, axis=(1, 2, 3)) > 6000)

        label_patches = label_patches[valid_idxs]
        x = np.vstack((x, np.zeros((len(label_patches), patch_shape_1d,
                                            patch_shape_1d, patch_shape_1d, num_mod))))

        T1_train = extract_patches(T1_vols[idx], patch_shape, extraction_step,datype="float32")
        x[x_length:, :, :, :, 0] = T1_train[valid_idxs]

        T2_train = extract_patches(T2_vols[idx], patch_shape, extraction_step,datype="float32")
        x[x_length:, :, :, :, 1] = T2_train[valid_idxs]
    return x

"""
To preprocess the unlabeled training data
"""
def preprocess_dynamic_unlab(dir,extraction_step,patch_shape,num_images_training_unlab):
    T1_vols = np.empty((num_images_training_unlab, 144, 192, 256),dtype="float32")
    T2_vols = np.empty((num_images_training_unlab, 144, 192, 256),dtype="float32")
    for case_idx in range(11, 11+num_images_training_unlab) :
        T1_vols[(case_idx - 11), :, :, :] = read_vol(case_idx, 'T1', dir)
        T2_vols[(case_idx - 11), :, :, :] = read_vol(case_idx, 'T2', dir)
        #print(read_vol(case_idx, 'T2', dir).shape)
    T1_mean = T1_vols.mean()
    T1_std = T1_vols.std()
    T1_vols = (T1_vols - T1_mean) / T1_std
    T2_mean = T2_vols.mean()
    T2_std = T2_vols.std()
    T2_vols = (T2_vols - T2_mean) / T2_std
    for i in range(T1_vols.shape[0]):
        T1_vols[i] = ((T1_vols[i] - np.min(T1_vols[i])) /
                                        (np.max(T1_vols[i])-np.min(T1_vols[i])))*255
    for i in range(T2_vols.shape[0]):
        T2_vols[i] = ((T2_vols[i] - np.min(T2_vols[i])) /
                                        (np.max(T2_vols[i])-np.min(T2_vols[i])))*255
    T1_vols = T1_vols/127.5 -1.
    T2_vols = T2_vols/127.5 -1.
    x=get_patches_unlab(T1_vols, T2_vols, extraction_step, patch_shape,dir)
    print("Total Extracted Unlabeled Patches Shape:",x.shape)
    return np.rollaxis(x, 4, 1)
