{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Setup","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom PIL import Image\nimport matplotlib.pyplot as plt\nfrom skimage import io\nimport glob, os, json, cv2, gc, shutil\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras import layers\nfrom tensorflow.keras import backend as K\n\nIMG_SIZE = 1024\nINPUT_PATH = \"../input/hubmap-kidney-segmentation\"\nPRINT_PLOTS = False\nFULL_RUN = True\n\ndef create_folder(folder):\n    if not os.path.exists(folder):\n        os.makedirs(folder)\n        \ncreate_folder(\"./plots\")\ncreate_folder(\"./train\")\ncreate_folder(\"./test\")\n\ndf_train = pd.read_csv(os.path.join(INPUT_PATH, 'train.csv'))\ndisplay(df_train)\n\ndf_image_info = pd.read_csv(os.path.join(INPUT_PATH,'HuBMAP-20-dataset_information.csv'))\ndisplay(df_image_info)","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:51.72996Z","iopub.execute_input":"2021-06-08T15:34:51.730407Z","iopub.status.idle":"2021-06-08T15:34:58.655479Z","shell.execute_reply.started":"2021-06-08T15:34:51.730292Z","shell.execute_reply":"2021-06-08T15:34:58.654499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Image Functions","metadata":{}},{"cell_type":"code","source":"def read_tiff(image_path):\n    image = io.imread(image_path)\n    image = np.squeeze(image) # some images have unnecessary axes with shape 1 --> remove\n    if image.shape[0] == 3: # some images have color as first axis -> swap axes\n        image = image.swapaxes(0,1)\n        image = image.swapaxes(1,2)\n    return image\n\ndef read_mask(image, encoded_mask):\n    mask = rle_decode(encoded_mask, (image.shape[1], image.shape[0])) # with inverted axes\n    mask = mask.swapaxes(0,1) # swap back axes\n    mask = np.expand_dims(mask, -1) # add one axis to have same shape as images\n    return mask\n\ndef delete_directory_contents(dir):\n    for file in os.scandir(dir):\n        os.remove(file.path)\n        \ndef plot_masked_image(image, mask, name):\n    plt.imshow(image, interpolation='none')\n    plt.imshow(mask, cmap='jet', alpha=0.3, interpolation='none')\n    \n    plt.savefig(f\"./plots/{name}.png\", dpi = 1000)\n    plt.show()\n    \ndef slice_images(image_id, image, mask=[], folder=\"\"):\n    print('Slicing Image ' + image_id + ' ...')\n\n    possible_slices_x = image.shape[0] // IMG_SIZE\n    possible_slices_y = image.shape[1] // IMG_SIZE\n\n    for x in range(possible_slices_x):\n        for y in range(possible_slices_y):\n            image_slice = image[x * IMG_SIZE : (x+1) * IMG_SIZE, y * IMG_SIZE : (y+1) * IMG_SIZE]\n            \n            #if np.any(image_slice) and not (image_slice > 200).all(): # only process non-black and non-gray images --> no background images\n\n            if not len(mask) == 0:\n                mask_slice = mask[x * IMG_SIZE : (x+1) * IMG_SIZE, y * IMG_SIZE : (y+1) * IMG_SIZE] * 255\n                if 255 in mask_slice:\n                    cv2.imwrite(f\"./{folder}/{image_id}-imgslice.{x}.{y}.jpg\", image_slice)\n                    cv2.imwrite(f\"./{folder}/{image_id}-maskslice.{x}.{y}.png\", mask_slice.astype(int))\n            else:\n                cv2.imwrite(f\"./{folder}/{image_id}-imgslice.{x}.{y}.jpg\", image_slice)\n\n# ref.: https://www.kaggle.com/stainsby/fast-tested-rle\ndef rle_decode(mask_rle, shape):\n    '''\n    mask_rle: run-length as string formated (start length)\n    shape: (height,width) of array to return \n    Returns numpy array, 1 - mask, 0 - background\n\n    '''\n    s = mask_rle.split()\n    starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]\n    starts -= 1\n    ends = starts + lengths\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    for lo, hi in zip(starts, ends):\n        img[lo:hi] = 1\n    return img.reshape(shape)\n\n## ref.: https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\ndef rle_encode_less_memory(img):\n    pixels = img.T.flatten()\n    \n    # This simplified method requires first and last pixel to be zero\n    pixels[0] = 0\n    pixels[-1] = 0\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 2\n    runs[1::2] -= runs[::2]\n    \n    return ' '.join(str(x) for x in runs)","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:58.658787Z","iopub.execute_input":"2021-06-08T15:34:58.659052Z","iopub.status.idle":"2021-06-08T15:34:58.67646Z","shell.execute_reply.started":"2021-06-08T15:34:58.659026Z","shell.execute_reply":"2021-06-08T15:34:58.674707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Slice Training Images to Uniform Size","metadata":{}},{"cell_type":"code","source":"def slice_training_images(df_train):\n    if not FULL_RUN:\n        df_train = df_train.iloc[0:1, :]  # only use one training image for quicker debug runs\n    else:\n        df_train = df_train.iloc[1:, :]\n    for index, train_sample in df_train.iterrows():\n        image_id = train_sample['id']\n        encoded_mask = train_sample['encoding']\n\n        image_path = os.path.join(INPUT_PATH, f\"train/{image_id}.tiff\")\n        image = read_tiff(image_path)\n        mask = read_mask(image, encoded_mask)\n\n        if PRINT_PLOTS:\n            image_to_plot = cv2.resize(image, (0, 0), fx=0.25, fy=0.25)\n            mask_to_plot = cv2.resize((mask * 255).astype('float32'), (0, 0), fx=0.25, fy=0.25)\n            plot_masked_image(image_to_plot, mask_to_plot, image_id)\n\n        slice_images(image_id, image, mask, \"train\")","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:58.67839Z","iopub.execute_input":"2021-06-08T15:34:58.678762Z","iopub.status.idle":"2021-06-08T15:34:58.690238Z","shell.execute_reply.started":"2021-06-08T15:34:58.678725Z","shell.execute_reply":"2021-06-08T15:34:58.689327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loss and Metrics","metadata":{}},{"cell_type":"code","source":"# ref.: https://gist.github.com/CarloSegat/1a2816676c48607dac9dda38afe4f3d9\ndef weighted_binary_crossentropy(y_true, y_pred, weight1=5, weight0=1):\n    y_true = K.clip(y_true, K.epsilon(), 1-K.epsilon())\n    y_pred = K.clip(y_pred, K.epsilon(), 1-K.epsilon())\n    logloss = -(y_true * K.log(y_pred) * weight1 + (1 - y_true) * K.log(1 - y_pred) * weight0 )\n    return K.mean( logloss, axis=-1)","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:58.69177Z","iopub.execute_input":"2021-06-08T15:34:58.692244Z","iopub.status.idle":"2021-06-08T15:34:58.704052Z","shell.execute_reply.started":"2021-06-08T15:34:58.692205Z","shell.execute_reply":"2021-06-08T15:34:58.703071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ref.: https://github.com/keras-team/keras/issues/3611\ndef dice_coef(y_true, y_pred, smooth=1):\n    intersection = K.sum(y_true * y_pred, axis=[1,2,3])\n    union = K.sum(y_true, axis=[1,2,3]) + K.sum(y_pred, axis=[1,2,3])\n    return K.mean( (2. * intersection + smooth) / (union + smooth), axis=0)","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:58.706879Z","iopub.execute_input":"2021-06-08T15:34:58.707379Z","iopub.status.idle":"2021-06-08T15:34:58.715329Z","shell.execute_reply.started":"2021-06-08T15:34:58.707339Z","shell.execute_reply":"2021-06-08T15:34:58.714614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# U-Net Model for our Kidney Images\nReference: https://arxiv.org/abs/1505.04597","metadata":{}},{"cell_type":"code","source":"def unet():\n    input = layers.Input(shape=(IMG_SIZE, IMG_SIZE, 3 ))\n\n    contraction1 = layers.Conv2D(32, 3, activation='relu', padding='same')(input)\n    contraction1 = layers.Conv2D(32, 3, activation='relu', padding='same')(contraction1)\n    pooling1 = layers.MaxPooling2D(2)(contraction1)\n\n    contraction2 = layers.Conv2D(64, 3, activation='relu', padding='same')(pooling1)\n    contraction2 = layers.Conv2D(64, 3, activation='relu', padding='same')(contraction2)\n    pooling2 = layers.MaxPooling2D(2)(contraction2)\n\n    contraction3 = layers.Conv2D(128, 3, activation='relu', padding='same')(pooling2)\n    contraction3 = layers.Conv2D(128, 3, activation='relu', padding='same')(contraction3)\n    pooling3 = layers.MaxPooling2D(2)(contraction3)\n\n    contraction4 = layers.Conv2D(256, 3, activation='relu', padding='same')(pooling3)\n    contraction4 = layers.Conv2D(256, 3, activation='relu', padding='same')(contraction4)\n    pooling4 = layers.MaxPooling2D(2)(contraction4)\n    \n    contraction5 = layers.Conv2D(512, 3, activation='relu', padding='same')(pooling4)\n    contraction5 = layers.Conv2D(256, 3, activation='relu', padding='same')(contraction5)\n\n    up_sampling1 = layers.UpSampling2D(2)(contraction5)\n    expansion1 = layers.concatenate([up_sampling1, contraction4], axis=-1)\n    expansion1 = layers.Conv2D(512, 3, activation='relu', padding='same')(expansion1)\n    expansion1 = layers.Conv2D(128, 3, activation='relu', padding='same')(expansion1)\n    \n    up_sampling2 = layers.UpSampling2D(2)(expansion1)\n    expansion2 = layers.concatenate([up_sampling2, contraction3], axis=-1)\n    expansion2 = layers.Conv2D(256, 3, activation='relu', padding='same')(expansion2)\n    expansion2 = layers.Conv2D(64, 3, activation='relu', padding='same')(expansion2)\n    \n    up_sampling3 = layers.UpSampling2D(2)(expansion2)\n    expansion3 = layers.concatenate([up_sampling3, contraction2], axis=-1)\n    expansion3 = layers.Conv2D(128, 3, activation='relu', padding='same')(expansion3)\n    expansion3 = layers.Conv2D(32, 3, activation='relu', padding='same')(expansion3)\n    \n    up_sampling4 = layers.UpSampling2D(2)(expansion3)\n    expansion4 = layers.concatenate([up_sampling4, contraction1], axis=-1)\n    expansion4 = layers.Conv2D(64, 3, activation='relu', padding='same')(expansion4)\n    expansion4 = layers.Conv2D(32, 3, activation='relu', padding='same')(expansion4)\n    output = layers.Conv2D(1, 1, activation='sigmoid', name='output', padding='same')(expansion4)\n    \n    learning_rate = 5e-4\n    optimizer = tf.keras.optimizers.Adam(learning_rate)\n    model = keras.Model(inputs=[input], outputs=[output])\n    losses = {'output': weighted_binary_crossentropy}\n    metrics = {'output': [\"acc\", dice_coef]}\n    model.compile(optimizer=optimizer, loss = losses, metrics=metrics)\n    model.summary()\n    \n    if PRINT_PLOTS:\n        keras.utils.plot_model(model, to_file=\"./model.png\", show_shapes=True)\n\n    return model","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:58.716927Z","iopub.execute_input":"2021-06-08T15:34:58.717363Z","iopub.status.idle":"2021-06-08T15:34:58.735681Z","shell.execute_reply.started":"2021-06-08T15:34:58.717321Z","shell.execute_reply":"2021-06-08T15:34:58.734495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mapping Functions to Read Sliced Images and Convert them to Tensors","metadata":{}},{"cell_type":"code","source":"def image_to_tensor(image_path, is_png):\n    image = tf.io.read_file(image_path)\n    if is_png:\n        image = tf.image.decode_png(image)\n    else:\n        image = tf.image.decode_jpeg(image)\n        \n    image = tf.image.convert_image_dtype(image, tf.float32)\n    \n    return image\n\n\ndef prepare_images_for_unet(image_path, mask_path):\n    image = image_to_tensor(image_path, is_png=False)\n    mask = image_to_tensor(mask_path, is_png=True)\n\n    return image, mask\n\n\ndef prepare_image_for_unet(image_path):\n    image = image_to_tensor(image_path, is_png=False)\n\n    return image","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:58.737208Z","iopub.execute_input":"2021-06-08T15:34:58.737628Z","iopub.status.idle":"2021-06-08T15:34:58.748351Z","shell.execute_reply.started":"2021-06-08T15:34:58.73759Z","shell.execute_reply":"2021-06-08T15:34:58.747462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Training","metadata":{}},{"cell_type":"code","source":"def train_model():\n    image_slice_paths = sorted(glob.glob('../working/train/*.jpg'))\n    mask_slice_paths = sorted(glob.glob('../working/train/*.png'))\n    \n    train_dataset = tf.data.Dataset.from_tensor_slices((image_slice_paths, mask_slice_paths)).shuffle(len(image_slice_paths)).map(prepare_images_for_unet).batch(16).prefetch(tf.data.experimental.AUTOTUNE)\n    print(len(image_slice_paths))\n    if PRINT_PLOTS:\n        for image_index in range(110, 130):\n            image_slice_path = image_slice_paths[image_index]\n            mask_slice_path = mask_slice_paths[image_index]\n\n            image = cv2.imread(image_slice_path)\n            mask = cv2.imread(mask_slice_path)\n\n            plot_masked_image(image, mask, f\"slice-{image_index}\")\n            del image\n            del mask\n    train_dataset = tf.data.Dataset.from_tensor_slices((image_slice_paths, mask_slice_paths)).shuffle(\n        len(image_slice_paths)).map(prepare_images_for_unet).batch(1).prefetch(tf.data.experimental.AUTOTUNE)\n    epochs = 30\n    verbose = 0\n    model = unet()\n    model.fit(train_dataset,\n              epochs=epochs,\n              verbose=verbose)\n    model.save('./model.h5')\n    del train_dataset\n    del image_slice_paths\n    del mask_slice_paths","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:58.74983Z","iopub.execute_input":"2021-06-08T15:34:58.750241Z","iopub.status.idle":"2021-06-08T15:34:58.761071Z","shell.execute_reply.started":"2021-06-08T15:34:58.750199Z","shell.execute_reply":"2021-06-08T15:34:58.760313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slice_training_images(df_train)","metadata":{"execution":{"iopub.status.busy":"2021-06-08T15:34:59.747583Z","iopub.execute_input":"2021-06-08T15:34:59.747929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_model()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cleanup","metadata":{}},{"cell_type":"code","source":"shutil.rmtree('./train')\nshutil.rmtree('./test')\nshutil.rmtree('./plots')","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}