<jupyter_start><jupyter_text>Monomodal Affine Registration using Convolutional Neural Networks  In this notebook a Convolutional Neural Network (CNN) is going to be trained to perform an algorithm called image registration. The image registration algorithm consists in, given a reference image, align another image according to the reference. To align means to find a geometrical transformation which makes the other image more similar with the reference.The archtecture of the CNN was proposed by SLOAN et al. in 2018 [1] and is focused in register monomodal images.<jupyter_code>from keras.layers import Input, Concatenate, Conv2D, concatenate, Flatten, Dense, Dropout
from keras.models import Model
from keras.models import load_model
from keras.optimizers import Adadelta
from keras.utils import Sequence
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import SimpleITK as sitk
import os<jupyter_output><empty_output><jupyter_text>Creating some useful functions We are going to use these functions later in this notebook.<jupyter_code>def myshow(img, title=None, margin=0.05, dpi=80, size=(1, 3)):
    nda = sitk.GetArrayViewFromImage(img)
    spacing = img.GetSpacing()
        
    ysize = nda.shape[0]
    xsize = nda.shape[1]
      
    figsize = (1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi

    fig = plt.figure(title, figsize=figsize, dpi=dpi)
    ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])
    
    extent = (0, xsize*spacing[1], 0, ysize*spacing[0])
    
    t = ax.imshow(nda,
            extent=extent,
            interpolation='hamming',
            cmap='gray',
            origin='lower')
    
    if(title):
        plt.title(title)
        
    plt.show()

    extent = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
    fig.savefig(title + '.png', bbox_inches=extent)

#This function receives an image (Image class of sitk) and does a transformation that rotates and translates the image
def affine_rotate_translate(image, x_translation, y_translation, degrees):
    center = np.floor(np.array(image.GetSize())/2)
    transform = sitk.Euler2DTransform(center, np.deg2rad(degrees), (x_translation, y_translation))
    
    reference_image = image
    interpolator = sitk.sitkCosineWindowedSinc
    default_value = 0
    rotated_translated_image = sitk.Resample(image, reference_image, transform, interpolator, default_value)
    
    return rotated_translated_image

#This function returns the minimum and maximum of an image object
def get_max_min(image):
    max_min = sitk.MinimumMaximumImageFilter()
    max_min.Execute(image)
    
    max_value = max_min.GetMaximum()
    min_value = max_min.GetMinimum()
    
    return np.array([min_value, max_value])
#This class assists the training in the batch form
class My_Custom_Generator(Sequence) :
  
  def __init__(self, image_filenames, labels, batch_size) :
    self.image_filenames = image_filenames
    self.labels = labels
    self.batch_size = batch_size
    
    
  def __len__(self) :
    return (np.ceil(len(self.image_filenames) / float(self.batch_size))).astype(np.int)
  
  
  def __getitem__(self, idx) :
    batch_x = self.image_filenames[idx * self.batch_size : (idx+1) * self.batch_size]
    batch_y = self.labels[idx * self.batch_size : (idx+1) * self.batch_size]
    
    return [np.expand_dims(np.array([sitk.GetArrayFromImage(sitk.ReadImage(str(file_name) + '/reference.mha')) 
                     for file_name in batch_x]), axis=3), np.expand_dims(np.array([sitk.GetArrayFromImage(sitk.ReadImage(str(file_name) + 
                     '/moving.mha')) for file_name in batch_x]), axis=3)], np.array(batch_y)<jupyter_output><empty_output><jupyter_text>Knowing the MRI Database The database used in this project is the OASIS-1: Cross-sectional MRI Data in Young, Middle Aged, Nondemented and Demented Older Adults (https://www.oasis-brains.org/data). This database is composed of 416 subjects and each of these ones have 3 or 4 images. or 20 of the nondemented subjects, images from a subsequent scan session after a short delay (less than 90 days) are also included as a means of assessing acquisition reliability. For more information see: https://www.oasis-brains.org/files/oasis_cross-sectional_facts.pdf<jupyter_code>volume_names = pd.read_csv('Data/oasis_cross-sectional.csv')

volume_names = volume_names.iloc[0:416, :]
volume_names['ID']<jupyter_output><empty_output><jupyter_text>Creating the Database for Registration The data used in the training step will be obtained by selecting 121 random slices in each of the 416 volumes and applying a rigid transformation to them, which consists in three degrees of freedom: An translation in X and Y, and a rotation around the center of the image. The selected slice is chosen in order to avoid slices that are composed mostly by air. After that, a Gaussian Noise is added to the trasformed image.At the end, the dataset is composed by 50336 pairs of image. 30201 will be used for training, 10067 for validation and the rest 10068 for testing.<jupyter_code>'''#Selecting 121 slices in each of the 416 volumes, because we want 50'000 pairs Reference and Moving
path_names = [] #save the path of the images
parameters = [] #Saving the parameters of the transformation in a sequence X_translation Y_translation Rotation_degree

for name in volume_names['ID']:
    volume = sitk.ReadImage('Data/Volumes/' + name + '/RAW/' + name + '_mpr-1_anon.hdr', sitk.sitkUInt16)

    for i in range(121):
        print("Volume " + name + ", Slice " + str(i))
        slc = np.random.randint(0, 127)
        reference = volume[:, :, slc]

        #Verifying if the selected slice is primarily composed  by air
        while np.sum(sitk.GetArrayFromImage(reference)) < 1e6:
            slc = np.random.randint(0, 127)
            reference = volume[:, :, slc]
            
        #Applying the transformation
        x = np.random.uniform(-30, 30)
        y = np.random.uniform(-30, 30)
        degrees = np.random.uniform(-15, 15)
        transformed = affine_rotate_translate(reference, x, y, degrees)
        moving = sitk.AdditiveGaussianNoise(transformed, standardDeviation=0.01)
        
        #Saving the Reference, Moving and transformation parameters
        path_folder = "Data/Register_Data_Base/" + name + "/Case_" + str(i)
        if not os.path.exists(path_folder):
            os.makedirs(path_folder)
        
        path_names.append(path_folder)
        reference_path = path_folder + "/reference.mha"
        moving_path = path_folder + "/moving.mha"
        
        sitk.WriteImage(reference, reference_path)
        sitk.WriteImage(moving, moving_path)
        
        #Saving the parameters
        parameters.append([-x, -y, -degrees])
        
parameters_df = pd.DataFrame(data=parameters, index=path_names, columns=['x', 'y', 'degrees'], dtype=np.float64, )
parameters_df.to_csv(os.getcwd() + "/Data/Register_Data_Base/transformation_parameters.csv", encoding='utf-8')'''<jupyter_output><empty_output><jupyter_text>Creating the Convolutional Neural Network The CNN archtecture is presented below. It basically consists in two 2D convolutional layers concatenated which selects some features necessary for the registration. After these layers, the output is inserted as input of a series of three dense layers. The output of the whole NN is three parameters, the X, Y for the translation and d degrees for the rotation.<jupyter_code>'''#Be careful: Input pixel range - 0 -- 1; concatenation approach, flatten approach
reference_input = Input(shape=(256, 256, 1), name='Reference')
moving_input = Input(shape=(256, 256, 1), name='Moving')

shared_conv1 = Conv2D(1, (3, 3), activation='relu', name='Shared_Conv1')
filtered_reference1 = shared_conv1(reference_input)
filtered_moving1 = shared_conv1(moving_input)

shared_conv2 = Conv2D(1, (3, 3), activation='relu', name='Shared_Conv2')
filtered_reference2 = shared_conv2(filtered_reference1)
filtered_moving2 = shared_conv2(filtered_moving1)

merged_vector = concatenate([filtered_reference2, filtered_moving2], axis=-1)
flatten_layer = Flatten()
flatten_vector = flatten_layer(merged_vector)

drop_layer = Dropout(0.5)
droped_merged_vector = drop_layer(flatten_vector)

densed_output1 = Dense(64, activation='relu', name='Dense_64_1')(droped_merged_vector)
densed_output2 = Dense(64, activation='relu', name='Dense_64_2')(densed_output1)
parameters = Dense(3, activation='linear', name='Dense_3')(densed_output2)

model = Model(inputs=[reference_input, moving_input], outputs=parameters)
model.compile(optimizer=Adadelta(1, 0.95), loss='mean_squared_error')

model.summary()'''<jupyter_output><empty_output><jupyter_text>Because the training set is very large, the training step will be done in batches. So, in order to do that, we create some keras sequences that will feed the CNN with the batches. The validation and test sets will also be inserted in the algorithm in the batch approach due to their sizes.<jupyter_code>#Reading the filenames
data_filenames = pd.read_csv('Data/Register_Data_Base/transformation_parameters.csv')

#Percentage of the whole data that corresponds to each set (training, validation and test)
training_rate = 0.6
validation_rate = 0.2
test_rate = 1 - training_rate - validation_rate

training_indice = int(training_rate*data_filenames.shape[0])
validation_rate = int(training_indice + validation_rate*data_filenames.shape[0])
test_rate = int(validation_rate + test_rate*data_filenames.shape[0])

X_train_filenames = data_filenames.iloc[0:training_indice, 0].to_numpy()
y_train = data_filenames.iloc[0:training_indice, 1:].to_numpy()

X_val_filenames = data_filenames.iloc[training_indice:validation_rate, 0].to_numpy()
y_val = data_filenames.iloc[training_indice:validation_rate, 1:].to_numpy()

X_test_filenames = data_filenames.iloc[validation_rate:test_rate, 0].to_numpy()
y_test = data_filenames.iloc[validation_rate:test_rate, 1:].to_numpy()

#Creating the instances for the batch training
batch_size = 32

my_training_batch_generator = My_Custom_Generator(X_train_filenames, y_train, batch_size)
my_validation_batch_generator = My_Custom_Generator(X_val_filenames, y_val, batch_size)
my_test_batch_generator = My_Custom_Generator(X_test_filenames, y_test, batch_size)<jupyter_output><empty_output><jupyter_text>Training the CNN<jupyter_code>'''hist = model.fit_generator(generator=my_training_batch_generator,
                           steps_per_epoch = int(len(X_train_filenames) // batch_size),
                           epochs = 30,
                           verbose = 1,
                           validation_data = my_validation_batch_generator,
                           validation_steps = int(len(X_val_filenames) // batch_size))'''
'''plt.plot(hist.history['loss'], label='Training Error')
plt.plot(hist.history['val_loss'], label='Validation Error')
plt.xlabel('Epochs')
plt.ylabel('MSE')
plt.title('MSE x Epochs')

plt.legend(loc='upper right')
plt.savefig('MSE_Epochs_Evolution.png')'''<jupyter_output><empty_output><jupyter_text>Saving the Neural Network Architecture and Parameters<jupyter_code>#model.save('CNN_2D_Registration_Monomodal_without_padding_intermediate.h5')<jupyter_output><empty_output><jupyter_text>Loading the Model<jupyter_code>model = load_model('CNN_2D_Registration_Monomodal_without_padding_intermediate.h5')<jupyter_output><empty_output><jupyter_text>Analyzing the CNN in the Test Data<jupyter_code>#Why the NN asks for a 4D image?

#Reading the test examples
test_reference = sitk.ReadImage('Data/Register_Data_Base/OAS1_0368_MR1/Case_96/reference.mha')
test_moving = sitk.ReadImage('Data/Register_Data_Base/OAS1_0368_MR1/Case_96/moving.mha')

test_reference_array = np.expand_dims(np.expand_dims(sitk.GetArrayFromImage(test_reference), axis=3), axis=0)
test_moving_array = np.expand_dims(np.expand_dims(sitk.GetArrayFromImage(test_moving), axis=3), axis=0)

parameters = model.predict(x=my_test_batch_generator)
parameters = parameters.astype('float')
Error = parameters - y_test
mse = (np.square(Error)).mean()
print('Mean Square Error in the Test Data: ' + str(mse))
X_mean = Error[0].mean()
X_sd = Error[0].std()

Y_mean = Error[1].mean()
Y_sd = Error[1].std()

d_mean = Error[2].mean()
d_sd = Error[2].std()

print('X error: ' + str(X_mean) + ' +- ' + str(X_sd))
print('Y error: ' + str(Y_mean) + ' +- ' + str(Y_sd))
print('d error: ' + str(d_mean) + ' +- ' + str(d_sd))
fig, axs = plt.subplots(3, 1, figsize=(4, 12), sharey='none')

axs[0].scatter(y_test[:, 0], parameters[:, 0])
axs[0].set_xlabel('X Target')
axs[0].set_ylabel('X Predicted')

axs[1].scatter(y_test[:, 1], parameters[:, 1])
axs[1].set_ylabel('Y Predicted')
axs[1].set_xlabel('Y Target')

axs[2].scatter(y_test[:, 2], parameters[:, 2])
axs[2].set_xlabel('d Target')
axs[2].set_ylabel('d Predicted')

fig.tight_layout()

fig.savefig('Parameters_distribution.png')
reference_test = sitk.ReadImage('Data/Register_Data_Base/OAS1_0439_MR1/Case_62/reference.mha')
moving_test = sitk.ReadImage('Data/Register_Data_Base/OAS1_0439_MR1/Case_62/moving.mha')

reference_array = np.expand_dims(np.expand_dims(sitk.GetArrayFromImage(reference_test), axis = 3), axis = 0)
moving_array = np.expand_dims(np.expand_dims(sitk.GetArrayFromImage(moving_test), axis = 3), axis = 0)

predicted_parameters = model.predict({'Reference': reference_array, 'Moving': moving_array})
predicted_parameters = predicted_parameters.astype(float)
predicted_image = affine_rotate_translate(moving_test, predicted_parameters[0][0], predicted_parameters[0][1], predicted_parameters[0][2])

myshow(reference_test, 'Reference')
myshow(moving_test, 'Moving')
myshow(predicted_image, 'Predicted')

print(X_test_filenames[7831])
print('Predicted Parameters: ' + str(predicted_parameters))
print('Original Parameters: ' + str(y_test[7831]))<jupyter_output>/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:4: DeprecationWarning: Both axis > a.ndim and axis < -a.ndim - 1 are deprecated and will raise an AxisError in the future.
  after removing the cwd from sys.path.
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:5: DeprecationWarning: Both axis > a.ndim and axis < -a.ndim - 1 are deprecated and will raise an AxisError in the future.
  """
{'fraction_non_alphanumeric': 0.10796211855489302, 'fraction_numerical': 0.02076464398456682, 'mean_word_length': 4.932584269662922, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 42, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '562100', 'license': 'MIT', 'n_tokens_mistral': 4875, 'n_tokens_neox': 4597, 'n_words': 1134, 'path': 'Project/Projeto_IA006.ipynb', 'repo_name': 'felippe-trigueiro/Machine-Learning', 'size': '14255'}