{"cells":[{"metadata":{},"cell_type":"markdown","source":"# High-res samples into multi-input CNN\nThis notebook will go through the steps required to take some high-resolution samples from an image and feed them into a multi-input CNN. The model will never see the entire image and only three 256x256 pixels samples. This approach assumes that there is more value to learn from high-resolution local areas in an image than a resized version of the entire image. To do so, we will use a custom generator which will prepare the images \"live\" as we build the batches during training.\nIn previous versions of this notebook, I tried to stack the three samples together and feed as one input into a ResNet50 but the model never learnt.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Version notes:\n* **V7**: Stacked images randonly sampled and feed them into a ResNet50 - *score:0.0*\n* **V8**: Randomly samples 3 images at full resolution and feed them into a multi-input CNN - *score:0.15*\n* **V9**: Apply to the model 3 times over each image at prediction time - *score:0.32*\n* **V11**: Take the mean from the predicted isup score over 5 runs instead of the max over 3 runs - *score:0.34*\n* **V12**: Now train the input branches first, before training a model with 3 input branches, with the input branches mostly frozen. \n* **V13**: It appeared that the model training was sometimes unstable as seen in V12's log. This version brings so minor changes to try to fix that.\n* **V14**: Reduce the number of epochs in order to run in less than 6 hours and submit the notebook.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport rasterio\nfrom sklearn.utils import shuffle\nimport openslide\n\nimport os\nimport sys\nfrom shutil import copyfile, move\nfrom tqdm import tqdm\nimport h5py\nimport random\nfrom random import randint\n\nimport tensorflow as tf\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom tensorflow.keras import layers\nfrom tensorflow.keras import Input\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import InputLayer, Input\nfrom tensorflow.keras.layers import Conv2D, Dense, Flatten, Dropout, Activation\nfrom tensorflow.keras.layers import BatchNormalization, Reshape, MaxPooling2D, GlobalAveragePooling2D\nfrom tensorflow.keras.callbacks import Callback\nfrom tensorflow.keras.callbacks import ModelCheckpoint, ReduceLROnPlateau, EarlyStopping\nfrom tensorflow.keras.applications import ResNet50, VGG16\nfrom keras.losses import mean_squared_error\nimport keras as K\nfrom sklearn.metrics import cohen_kappa_score","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(\"../input/prostate-cancer-grade-assessment/train.csv\")\nimage_path = \"../input/prostate-cancer-grade-assessment/train_images/\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(train_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Generate paths to the input files and a list with corresponding groundtruth","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"image_size = 256\ntraining_sample_percentage = 0.9\ntraining_item_count = int(len(train_df)*training_sample_percentage)\ntrain_df[\"image_path\"] = [image_path+image_id+\".tiff\" for image_id in train_df[\"image_id\"]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#remove all image file that don't have a mask file\nindex_to_drop = []\nfor idx, row in train_df.iterrows():\n    mask_path = row.image_path.replace(\"train_images\",\"train_label_masks\").replace(\".tiff\",\"_mask.tiff\")\n\n    if not os.path.isfile(mask_path):\n        index_to_drop.append(idx)\n\ntrain_df.drop(index_to_drop,0,inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = openslide.OpenSlide(train_df.iloc[0].image_path)\nprint(example.dimensions)\nclipped_example = example.read_region((5000, 5000), 0, (256, 256))\nplt.imshow(clipped_example)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We're quickly shuffling the data and split the dataframe in two in order to keep a validation set.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = shuffle(train_df)\nvalidation_df = train_df[training_item_count:]\ntrain_df = train_df[:training_item_count]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Image preparation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_single_sample(image_path,image_size=256,training=False,display=False):\n    '''\n    Return a single 256x256 sample\n    with possibility of returning a gleason score using the masks\n    '''\n    \n    image = openslide.OpenSlide(image_path)\n    \n    mask_path = image_path.replace(\"train_images\",\"train_label_masks\").replace(\".tiff\",\"_mask.tiff\")\n    mask = openslide.OpenSlide(mask_path)\n    \n    stacked_image = []\n    groundtruth_per_image = []\n    \n    maximum_iteration = 0\n    selected_sample = False\n    while not selected_sample:\n        sampling_start_x = randint(image_size,image.dimensions[0]-image_size)\n        sampling_start_y = randint(image_size,image.dimensions[1]-image_size)\n\n        clipped_sample = image.read_region((sampling_start_x, sampling_start_y), 0, (256, 256))\n        clipped_array = np.asarray(clipped_sample)\n        \n        #check that the sample is not empty\n        #and use the standard deviation to make sure\n        #there is something happening in the sample\n        if (not np.all(clipped_array==255) and np.std(clipped_array)>20) or maximum_iteration>200:\n            if display:\n                plt.imshow(clipped_sample)\n                plt.show()\n                \n            sampled_image = clipped_array[:,:,:3]\n            \n            if training:\n                clipped_mask = mask.read_region((sampling_start_x, sampling_start_y), 0, (256, 256))\n                groundtruth_per_image.append(np.mean(np.asarray(clipped_mask)[:,:,0]))\n            \n            selected_sample = True\n        maximum_iteration+=1\n    \n    if training: \n        return np.array(sampled_image), np.array(groundtruth_per_image)\n    else:\n        return np.array(sampled_image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_random_samples(image_path,image_size=256,display=False):\n    '''\n    Load an image and select random areas.\n    Return a list of 3 images from areas where there is data.\n    '''\n    \n    image = openslide.OpenSlide(image_path)\n    stacked_image = []\n    \n    selected_samples = 0\n    maximum_iteration = 0\n    while selected_samples<3:\n        sampling_start_x = randint(image_size,image.dimensions[0]-image_size)\n        sampling_start_y = randint(image_size,image.dimensions[1]-image_size)\n\n        clipped_sample = image.read_region((sampling_start_x, sampling_start_y), 0, (256, 256))\n        clipped_array = np.asarray(clipped_sample)\n        \n        #check that the sample is not empty\n        #and use the standard deviation to make sure\n        #there is something happening in the sample\n        if (not np.all(clipped_array==255) and np.std(clipped_array)>20) or maximum_iteration>200:\n            if display:\n                plt.imshow(clipped_sample)\n                plt.show()\n\n            stacked_image.append(clipped_array[:,:,:3])\n            selected_samples+=1\n        maximum_iteration+=1\n    return np.array(stacked_image)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Just to double-check. We know generate images with 3 images with 3 channels each.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"get_random_samples(train_df.iloc[0].image_path).shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can check that the function works by displaying the random samples that will be used as inputs.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = get_random_samples(train_df.iloc[0].image_path, display=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"output = get_single_sample(train_df.iloc[0].image_path, display=True, training=True)\nprint(output[1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We create a custom generator that will randomly take samples from each full-scale image. It ensures there is actually some kind of content in the returned sample and give the associated groundtruth.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def custom_single_image_generator(image_path_list, batch_size=16):\n    '''\n    return an image and a corresponding gleason score from the mask\n    '''\n    \n    while True:\n        for start in range(0, len(image_path_list), batch_size):\n            X_batch = []\n            Y_batch = []\n            end = min(start + batch_size, training_item_count)\n\n            image_info_list = [get_single_sample(image_path, training=True) for image_path in image_path_list[start:end]]\n            X_batch = np.array([image_info[0]/255. for image_info in image_info_list])\n            Y_batch = np.array([image_info[1] for image_info in image_info_list])\n            \n            yield X_batch, Y_batch","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train our future input branch","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We build a fairly \"basic\" CNN that will be trained using random 256x256pixel samples and the gleason score. The gleason score will be calculcated from the corresponding masks for each image.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"num_channel = 3\nimage_shape = (image_size, image_size, num_channel)\n\ndef branch(input_image):\n    x = Conv2D(128, (3, 3))(input_image)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = Conv2D(128, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = Conv2D(128, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = MaxPooling2D(pool_size=(2, 2))(x)\n    \n    x = Conv2D(64, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = Conv2D(64, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = Conv2D(64, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = MaxPooling2D(pool_size=(2, 2))(x)\n    \n    x = Conv2D(32, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = Conv2D(32, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = Conv2D(32, (3, 3))(x)\n    x = BatchNormalization()(x)\n    x = Activation('relu')(x)\n    x = GlobalAveragePooling2D()(x)\n    \n    x = layers.Dense(256)(x)\n    x = Activation('relu')(x)\n    \n    return layers.Dropout(0.3)(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input_image = Input(shape=image_shape)\ncore_branch = branch(input_image)\noutput = Dense(1, activation='linear')(core_branch)\n\nbranch_model = Model(input_image,output)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"branch_model.compile(loss=\"mse\",\n                      optimizer=tf.keras.optimizers.RMSprop(learning_rate=0.001))\ncallbacks = [ReduceLROnPlateau(monitor='val_loss', patience=1, verbose=1, factor=0.5),\n             EarlyStopping(monitor='val_loss', patience=3),\n             ModelCheckpoint(filepath='best_branch.h5', monitor='val_loss', save_best_only=True)]\n\nbatch_size = 16\n\nhistory = branch_model.fit_generator(custom_single_image_generator(train_df[\"image_path\"], batch_size=batch_size),\n                        steps_per_epoch = int(len(train_df)/batch_size),\n                        validation_data=custom_single_image_generator(validation_df[\"image_path\"], batch_size=batch_size),\n                        validation_steps= int(len(validation_df)/batch_size),\n                        epochs=2,\n                        callbacks=callbacks)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Build the core model","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Below, we have our custom generator which will use the function to generate original-scale 256x256 images from random samples.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def custom_generator(image_path_list, groundtruth_list, batch_size=16):\n    num_classes=6\n    while True:\n        for start in range(0, len(image_path_list), batch_size):\n            X_batch = []\n            Y_batch = []\n            end = min(start + batch_size, training_item_count)\n            \n            X_batch = np.array([get_random_samples(image_path)/255. for image_path in image_path_list[start:end]])\n            input_image1 = X_batch[:,0,:,:]\n            input_image2 = X_batch[:,1,:,:]\n            input_image3 = X_batch[:,2,:,:]\n            \n            Y_batch = tf.keras.utils.to_categorical(np.array(groundtruth_list[start:end]),num_classes) \n            \n            yield [input_image1,input_image2,input_image3], Y_batch","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, time to build the model. We have an `input_branch` function that will generate the 3 input branches with our pretrained weights, before merging them together. The 3 input branches have their weight mostly frozen as they have already been trained using the gleason score. The architecture itself will need some additional tuning later. The activation layer for our output is set to `softmax` with 6 units, representing the gradient from 0 to 5. It could also be set to `linear` if we approached this as a regression problem.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def input_branch(input_image):\n    '''\n    Generate a new input branch using our previous weights\n    \n    '''\n    input_image = Input(shape=image_shape)\n    core_branch = branch(input_image)\n    output = Dense(1, activation='linear')(core_branch)\n    branch_model = Model(input_image,output)\n    branch_model.load_weights(\"../working/best_branch.h5\")\n        \n    new_branch = Model(inputs=branch_model.input, outputs=branch_model.layers[-2].output)\n    \n    for layer in new_branch.layers[:-3]:\n        layer.trainable = False\n    \n    return new_branch","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input_image1 = Input(shape=image_shape)\ninput_image2 = Input(shape=image_shape)\ninput_image3 = Input(shape=image_shape)\n\nfirst_branch = branch(input_image1)\nsecond_branch = branch(input_image2)\nthird_branch = branch(input_image3)\n\nmerge = layers.Concatenate(axis=-1)([first_branch,second_branch,third_branch])\ndense = layers.Dense(256)(merge)\ndropout = layers.Dropout(0.3)(dense)\noutput = Dense(6, activation='softmax')(dropout)\n\nmodel = Model([input_image1,input_image2,input_image3],output)\nmodel.compile(loss='categorical_crossentropy',\n              optimizer=tf.keras.optimizers.Adam(learning_rate=0.0001))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Time to use our custom generator which loads 3 samples into the branches. We need to use it both for our training data and validation data. Be careful if you decide to add data augmentation into your custom generator as your performance on the validation set will be biased. In this case, add a parameter to your custom generator to activate the augmentation.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"callbacks = [ReduceLROnPlateau(monitor='val_loss', patience=1, verbose=1, factor=0.5),\n             EarlyStopping(monitor='val_loss', patience=3),\n             ModelCheckpoint(filepath='best_model.h5', monitor='val_loss', save_best_only=True)]\n\nbatch_size = 16\nhistory = model.fit_generator(custom_generator(train_df[\"image_path\"], train_df[\"isup_grade\"], batch_size=batch_size),\n                        steps_per_epoch = int(len(train_df)/batch_size),\n                        validation_data=custom_generator(validation_df[\"image_path\"],np.array(validation_df[\"isup_grade\"]), batch_size=batch_size),\n                        validation_steps=int(len(validation_df)/batch_size),\n                        epochs=3,\n                        callbacks=callbacks)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(history.history['loss'])\nplt.plot(history.history['val_loss'])\nplt.title('Loss over epochs')\nplt.ylabel('Loss')\nplt.xlabel('Epoch')\nplt.legend(['Train', 'Validation'], loc='best')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Make predictions using our newly trained model","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model.load_weights(\"best_model.h5\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Time to predict on new data! We use the function to retrieve random samples from the original full-size images and simply feed it to the model. The `passes` parameter will decide how many times we attempt to predict on a given image. Due to the random sampling, we will do 3 passes on each image to increase our chances of finding the cancerous zones and keep the mean of all grades. In future version, we will try to optimise the sampling strategy.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def predict_submission(df, path, passes=1):\n    \n    df[\"image_path\"] = [path+image_id+\".tiff\" for image_id in df[\"image_id\"]]\n    df[\"isup_grade\"] = 0\n    \n    for idx, row in df.iterrows():\n        prediction_per_pass = []\n        for i in range(passes):\n            model_input = np.array([get_random_samples(row.image_path)/255.])\n            input_image1 = model_input[:,0,:,:]\n            input_image2 = model_input[:,1,:,:]\n            input_image3 = model_input[:,2,:,:]\n\n            prediction = model.predict([input_image1,input_image2,input_image3])\n            prediction_per_pass.append(np.argmax(prediction))\n            \n        df.at[idx,\"isup_grade\"] = np.mean(prediction_per_pass)\n    df = df.drop('image_path', 1)\n    return df[[\"image_id\",\"isup_grade\"]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One last check to ensure that the `predict_submission` function works as intended. We use 100 images from the training set.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntest_from_training_df = pd.read_csv(\"../input/prostate-cancer-grade-assessment/train.csv\")[:20]\npredict_submission(test_from_training_df, image_path, passes=5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To finish, we add a condition that will check that the test images are available when running the notebook in `submission` mode. If it finds the folder, it will load `test.csv` and apply predictions before saving the submission file.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_path = \"../input/prostate-cancer-grade-assessment/test_images/\"\nsubmission_df = pd.read_csv(\"../input/prostate-cancer-grade-assessment/sample_submission.csv\")\n\nif os.path.exists(test_path):\n    test_df = pd.read_csv(\"../input/prostate-cancer-grade-assessment/test.csv\")\n    submission_df = predict_submission(test_df, test_path, passes=5)\n\nsubmission_df.to_csv('submission.csv', index=False)\nsubmission_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### If you found this notebook helpful, please give it an upvote. It is always greatly appreciated!","execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}