{"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":"# Begin","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom PIL import Image\nimport tensorflow as tf\nfrom random import seed\nimport matplotlib.pyplot as plt\nfrom skimage import color\nfrom tqdm.notebook import tqdm\nimport time\nimport os\nimport gc\n\ntqdm.pandas()\ngc.enable()\n\nt_size = 1024 #Size of the image to rescale\ntime_limit = 7 * 3600 #time limit in seconds\nstart_time = time.time()\n\nrandom_seed = 541\nseed(random_seed)\nnp.random.seed(random_seed)\ntf.random.set_seed(random_seed)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-08-09T00:13:53.357085Z","iopub.execute_input":"2022-08-09T00:13:53.357951Z","iopub.status.idle":"2022-08-09T00:13:58.259349Z","shell.execute_reply.started":"2022-08-09T00:13:53.357845Z","shell.execute_reply":"2022-08-09T00:13:58.258257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Supporting Functions\n\nThe following functions are used repeatedly across the whole notebook, so we will define them here.","metadata":{}},{"cell_type":"code","source":"def rle_to_mask(rle, width, target_size=None):\n    if target_size == None:\n        target_size = width\n\n    rle = np.array(list(map(int, rle.split())))\n    label = np.zeros((width*width))\n    \n    for start, end in zip(rle[::2], rle[1::2]):\n        label[start:start+end] = 1\n        \n    #Convert label to image\n    label = Image.fromarray(label.reshape(width, width))\n    #Resize label\n    label = label.resize((target_size, target_size))\n    label = np.array(label).astype(float)\n    #rescale label\n    label = np.round((label - label.min())/(label.max() - label.min()))\n    \n    return label.T\n\ndef mask_to_rle(mask, orig_dim=160):\n    #Rescale image to original size\n    size = int(len(mask.flatten())**.5)\n    n = Image.fromarray(mask.reshape((size, size))*255.0)\n    n = n.resize((orig_dim, orig_dim))\n    n = np.array(n).astype(np.float32)\n    #Get pixels to flatten\n    pixels = n.T.flatten()\n    #Round the pixels using the half of the range of pixel value\n    pixels = (pixels-min(pixels) > ((max(pixels)-min(pixels))/2)).astype(int)\n    pixels = np.nan_to_num(pixels) #incase of zero-div-error\n    \n    pixels = np.concatenate([[0], pixels, [0]])\n    runs = np.where(pixels[1:] != pixels[:-1])[0]\n    runs[1::2] -= runs[::2]\n    \n    return ' '.join(str(x) for x in runs)\n\ndef show_masked_img(img, mask, title=''):\n    #rescalse mask to 0-1 range regardless of min and max value\n    mask = (mask - mask.min())/(mask.max()-mask.min())\n    mask = np.nan_to_num(mask) #Adress NaN incase of zerodiv error\n    #Round up the mask\n    mask = np.round(mask)\n    mask = mask.reshape(img.shape[:2])\n    \n    fig, ax = plt.subplots(1, 3, figsize=(9, 3))\n    fig.suptitle(title, fontsize=16)\n    \n    ax[0].imshow(mask); ax[0].set_title('Mask')\n    ax[1].imshow(img); ax[1].set_title('Image')\n    ax[2].imshow(color.label2rgb(mask, img,\n                               bg_label=0, bg_color=(1.,1.,1.), alpha=0.25))\n    ax[2].set_title('Masked Image')\n    plt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-09T00:13:58.264959Z","iopub.execute_input":"2022-08-09T00:13:58.266005Z","iopub.status.idle":"2022-08-09T00:13:58.284410Z","shell.execute_reply.started":"2022-08-09T00:13:58.265966Z","shell.execute_reply":"2022-08-09T00:13:58.283602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv('../input/hubmap-organ-segmentation/train.csv')\ndisplay(train_df)","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:13:58.289536Z","iopub.execute_input":"2022-08-09T00:13:58.291797Z","iopub.status.idle":"2022-08-09T00:13:58.633603Z","shell.execute_reply.started":"2022-08-09T00:13:58.291751Z","shell.execute_reply":"2022-08-09T00:13:58.632414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"All length is equal to all width:\", np.all(train_df.img_height == train_df.img_width))","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:13:58.636117Z","iopub.execute_input":"2022-08-09T00:13:58.636388Z","iopub.status.idle":"2022-08-09T00:13:58.658858Z","shell.execute_reply.started":"2022-08-09T00:13:58.636363Z","shell.execute_reply":"2022-08-09T00:13:58.657254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"With the above snippet, we can confirm that all input images are square shaped. After careful observation, only the the `organ` and one `side` measure of the original image is useful for further training. Next, we confirm the helper functions if they work correctly.","metadata":{}},{"cell_type":"code","source":"rle = train_df.rle.iloc[0]\nl = train_df.img_width.iloc[0]\n# Test converting the rle back and forth and check for loss of data\nmask = rle_to_mask(rle, l)\nnrle = mask_to_rle(mask.flatten(), l)\nprint(\"Original is same as Reconstructed:\", rle == nrle)","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:13:58.660591Z","iopub.execute_input":"2022-08-09T00:13:58.661325Z","iopub.status.idle":"2022-08-09T00:14:01.099741Z","shell.execute_reply.started":"2022-08-09T00:13:58.661282Z","shell.execute_reply":"2022-08-09T00:14:01.098620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Loader","metadata":{}},{"cell_type":"code","source":"class ImageDataGenerator(tf.keras.utils.Sequence):\n    def __init__(self, df, batch_size=32, train=False):\n        self.df = df.reset_index(drop=True)\n        self.dim = t_size\n        self.train = train\n        if self.train: self.batch_size = batch_size // 4\n        else: self.batch_size = batch_size\n        self.pref = 'train' if train else 'test'\n    \n    def __len__(self):\n        return np.ceil(len(self.df) / self.batch_size).astype(int)\n    \n    def on_epoch_end(self):\n        if self.train: #Reshuffle train on end of epoch\n            self.df = self.df.sample(frac=1.0).reset_index(drop=True)\n            \n    def __getitem__(self, idx):\n        batch_x = self.df.iloc[idx*self.batch_size:(idx+1)*self.batch_size].id.values\n        \n        if not self.train:\n            X = np.zeros((batch_x.shape[0], self.dim, self.dim, 3))\n            \n            for i in range(batch_x.shape[0]):\n                image = Image.open(f\"../input/hubmap-organ-segmentation/{self.pref}_images/{batch_x[i]}.tiff\")\n                image = image.resize((self.dim, self.dim))\n                image = np.array(image) / 255.\n                X[i,] = image\n                \n            return X\n                \n        else:\n            batch_y = self.df.iloc[idx*self.batch_size:(idx+1)*self.batch_size].rle.values\n            batch_w = self.df.iloc[idx*self.batch_size:(idx+1)*self.batch_size].img_width.values\n            \n            X = np.zeros((batch_x.shape[0]*4, self.dim, self.dim, 3))\n            Y = np.zeros((batch_x.shape[0]*4, self.dim, self.dim, 1))\n            \n            for i in range(batch_x.shape[0]):\n                image = Image.open(f\"../input/hubmap-organ-segmentation/{self.pref}_images/{batch_x[i]}.tiff\")\n                image = image.resize((self.dim, self.dim))\n                image = np.array(image) / 255.\n                rle = rle_to_mask(batch_y[i], batch_w[i], self.dim)\n                rle = rle.reshape((self.dim, self.dim, 1))\n                \n                for n, (h, v) in enumerate([(0, 0), (0, 1), (1, 0), (1, 1)]):\n                    X[i*4 + n, :, :, :] = self.flip_image(self.swapaxes(image), h, v)[:, :, :]\n                    Y[i*4 + n, :, :, :] = self.flip_image(rle, h, v)[:, :, :]\n                    \n            return X, Y#.reshape(Y.shape[:-1])\n                \n    \n    def flip_image(self, image, horizontal=True, vertical=True):\n        n_image = image.copy()\n        \n        if horizontal:\n            n_image = n_image[:, ::-1, :]\n        if vertical:\n            n_image = n_image[::-1, :, :]\n        return n_image\n    \n    def swapaxes(self, o_image):\n        \"\"\"Randomly swap RGB channels to add variation/noise to train data.\"\"\"\n        swap_point_1 = np.random.choice([0, 1, 2, 3])\n        swap_point_2 = np.random.choice([c for c in [0, 1, 2, 3] if c != swap_point_1])\n            \n            \n        image = o_image.copy()\n        if 0 in [swap_point_1, swap_point_2]:\n            image[:, :, 0], image[:, :, 1] = image[:, :, 1], image[:, :, 0]\n        if 1 in [swap_point_1, swap_point_2]:\n            image[:, :, 0], image[:, :, 2] = image[:, :, 2], image[:, :, 0]\n        if 2 in [swap_point_1, swap_point_2]:\n            image[:, :, 1], image[:, :, 2] = image[:, :, 2], image[:, :, 1]\n        \n        return image","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:14:01.102149Z","iopub.execute_input":"2022-08-09T00:14:01.102804Z","iopub.status.idle":"2022-08-09T00:14:02.003573Z","shell.execute_reply.started":"2022-08-09T00:14:01.102762Z","shell.execute_reply":"2022-08-09T00:14:02.002621Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model\n\nWe will be creating a U-Net model using Keras Layers.","metadata":{}},{"cell_type":"code","source":"from tensorflow.keras.layers import Input, Conv2D, MaxPool2D, Dropout\nfrom tensorflow.keras.layers import Conv2DTranspose, concatenate\nfrom tensorflow.keras import Model\n\ndef double_conv_block(x, n_filters, step=3):\n    o = Conv2D(n_filters, step, padding='same', activation='leaky_relu', kernel_initializer='he_normal')(x)\n    o = Conv2D(n_filters, step, padding='same', activation='leaky_relu', kernel_initializer='he_normal')(o)\n    return o\n\ndef downsample_block(x, n_filters, step=3):\n    f = double_conv_block(x, n_filters, step)\n    p = MaxPool2D(2)(f)\n    p = Dropout(0.25)(p)\n    \n    return f, p\n\ndef upsample_block(x, conv_f, n_filters, step=3):\n    o = Conv2DTranspose(n_filters, step, 2, padding='same')(x)\n    o = concatenate([o, conv_f])\n    o = Dropout(0.25)(o)\n    o = double_conv_block(o, n_filters, step)\n    \n    return o\n\ndef create_model(in_shape, n=32):\n    s = 3\n    i = Input(shape=in_shape)\n    f1, p1 = downsample_block(i, n, s)\n    f2, p2 = downsample_block(p1, n*2, s)\n    f3, p3 = downsample_block(p2, n*2**2, s)\n    f4, p4 = downsample_block(p3, n*2**3, s)\n    \n    bridge = double_conv_block(p4, n*2**4, s)\n    \n    u4 = upsample_block(bridge, f4, n*2**3, s)\n    u3 = upsample_block(u4, f3, n*2**2, s)\n    u2 = upsample_block(u3, f2, n*2, s)\n    u1 = upsample_block(u2, f1, n, s)\n    \n    o = Conv2D(2, 1, padding='same', activation='softmax')(u1)\n    unet = Model(inputs=i, outputs=o)\n    return unet\n\nmodel = create_model((t_size, t_size, 3), 32)\nmodel.summary()\ntf.keras.utils.plot_model(model, show_shapes=True, rankdir='TB')","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:14:21.677282Z","iopub.execute_input":"2022-08-09T00:14:21.677845Z","iopub.status.idle":"2022-08-09T00:14:26.007852Z","shell.execute_reply.started":"2022-08-09T00:14:21.677806Z","shell.execute_reply":"2022-08-09T00:14:26.006811Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training","metadata":{"_kg_hide-input":true}},{"cell_type":"code","source":"from IPython.display import clear_output\nfrom tensorflow.keras import backend as K\n\nclass PlotLearning(tf.keras.callbacks.Callback):\n    \"\"\"\n    Callback to plot the learning curves of the model during training.\n    \"\"\"\n    def __init__(self, org_name='', **kwargs):\n        super().__init__(**kwargs)\n        self.org_name = org_name\n        \n    def on_train_begin(self, logs={}):\n        self.metrics = {}\n        for metric in logs:\n            self.metrics[metric] = []\n        \n    def on_epoch_end(self, epoch, logs={}):\n        # Storing metrics\n        for metric in logs:\n            if metric in self.metrics:\n                self.metrics[metric].append(logs.get(metric))\n            else:\n                self.metrics[metric] = [logs.get(metric)]\n        \n        # Plotting\n        metrics = [x for x in logs if 'val' not in x]\n        \n        f, axs = plt.subplots(1, len(metrics), figsize=(15,5))\n        clear_output(wait=True)\n        f.suptitle(f'Organ: {self.org_name}')\n        \n        for i, metric in enumerate(metrics):\n            axs[i].plot(range(1, epoch + 2), \n                        self.metrics[metric], \n                        label=metric)\n            if logs['val_' + metric]:\n                axs[i].plot(range(1, epoch + 2), \n                            self.metrics['val_' + metric], \n                            label='val_' + metric)\n                \n            axs[i].legend()\n            axs[i].grid()\n\n        plt.tight_layout()\n        plt.show()\n\ndef dice_coef(y_true, y_pred, smooth=1):\n    # flatten\n    y_true_f = K.flatten(K.cast(y_true, np.float32))\n    #y_pred_f = K.flatten(K.cast(K.argmax(y_pred, axis=-1), np.float32))\n    y_pred_f = K.flatten(K.cast(K.argmax(y_pred, axis=-1), np.float32))\n    # compute numerator and denominator\n    intersection = K.sum(y_true_f * y_pred_f)\n    union = K.sum(y_true_f) + K.sum(y_pred_f)\n    # apply dice formula\n    dice = K.mean((2.*intersection + smooth)/(union + smooth))\n    return dice\n\ndef weighted_loss(y_true, y_pred):\n    # Calculate the base loss\n    ce = K.sparse_categorical_crossentropy(y_true, y_pred)\n    # Apply the weights\n    one_weight = 1.0\n    zero_weight = 1e-2\n    weight_vector = y_true * one_weight + (1. - y_true) * zero_weight\n    weight_vector = K.squeeze(weight_vector, axis=-1)\n    weighted_ce = weight_vector * ce\n\n    # Return the mean error\n    return K.mean(weighted_ce)","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:14:26.010310Z","iopub.execute_input":"2022-08-09T00:14:26.010926Z","iopub.status.idle":"2022-08-09T00:14:26.030916Z","shell.execute_reply.started":"2022-08-09T00:14:26.010877Z","shell.execute_reply":"2022-08-09T00:14:26.029723Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split as tts\nfrom tensorflow.keras.callbacks import ModelCheckpoint\n\n\n#We will be training a model for each organ\norgans = train_df.organ.unique()\nepochs = 50\nn_kernels = 16\nbatch_size = 4\nbasic_lr = 3e-4\nweights_path = \"../input/hubmap-hpa-gpu-train-1024px-model\"","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:14:39.159592Z","iopub.execute_input":"2022-08-09T00:14:39.161854Z","iopub.status.idle":"2022-08-09T00:14:39.844079Z","shell.execute_reply.started":"2022-08-09T00:14:39.161816Z","shell.execute_reply":"2022-08-09T00:14:39.842950Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction","metadata":{}},{"cell_type":"code","source":"test_df = pd.read_csv('../input/hubmap-organ-segmentation/test.csv')\nsub = {'id':[], 'rle':[]}\n\nfor organ in organs:\n    print(\"Predicting on organ:\", organ)\n\n    test_X = test_df[test_df.organ == organ].reset_index(drop=True)\n    if len(test_X) == 0: continue #Skip organs without item in test\n        \n    model = create_model((t_size, t_size, 3), n_kernels)\n    model.load_weights(f\"{weights_path}/{organ}_model.hdf5\")\n    \n    test_loader = ImageDataGenerator(test_X, batch_size)\n    preds = model.predict(test_loader)\n    rle = [mask_to_rle(m, d) for m,d in zip(np.argmax(preds, axis=-1), test_X.img_width)]\n    sub['id'] += test_X.id.values.tolist()\n    sub['rle'] += rle\n    \nsub = pd.DataFrame(sub)\ntest_df = test_df.merge(sub, on='id')\nsub = test_df[['id', 'rle']].copy()\nsub.to_csv('submission.csv', index=False)\nsub","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:14:39.845420Z","iopub.execute_input":"2022-08-09T00:14:39.845775Z","iopub.status.idle":"2022-08-09T00:14:42.181885Z","shell.execute_reply.started":"2022-08-09T00:14:39.845740Z","shell.execute_reply":"2022-08-09T00:14:42.180701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vx = test_loader[0]\npvy = np.argmax(model.predict(vx), axis=-1)\nfor n, (a, b) in enumerate(zip(vx, pvy)):\n    show_masked_img(a, b, f\"Predicted Mask {n}\")\n    print(\"-\"*72)\n    if n == 2:\n        break","metadata":{"execution":{"iopub.status.busy":"2022-08-09T00:14:42.183374Z","iopub.execute_input":"2022-08-09T00:14:42.184372Z","iopub.status.idle":"2022-08-09T00:14:44.026054Z","shell.execute_reply.started":"2022-08-09T00:14:42.184330Z","shell.execute_reply":"2022-08-09T00:14:44.024950Z"},"trusted":true},"execution_count":null,"outputs":[]}]}