{"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":"# About \n\nThe idea is discussed [here, method 2](https://www.kaggle.com/c/siim-covid19-detection/discussion/245323), that is to add segmented roi or the predicted segmentation maps from segmentation model as an additional channel of model inputs. However, there are three possible ways to achieve this in terms of **Spatial** and **Channel** supervised approach, such as \n\n1. **Channel Supervision**: Use Segmentaiton model - get predictin maps - add it with inputs (multi-channel, 2 or 4). The classifier input can be 2 channel inpuut or 4 channel, depends. If choose either case, you may not use pre-trained weights directly but need to add one 3 filter conv layer with consistant padding. Also, we can drop 1 channel out of RGB and add the segmentatoin maps to complete 3 channel classifier input for pre-trained weights.\n\n![one](https://user-images.githubusercontent.com/17668390/122654678-e6940480-d16e-11eb-8db7-38fd7da7851b.png)\n\n2. **Spatial Supervision (1)**: Use Boundig-Box cordinates to draw rectangle on inputs. In the data loading time, draw bbox on the sample; maybe fill up with soft color probably can do some good. \n\n\n![two](https://user-images.githubusercontent.com/17668390/122654774-7a65d080-d16f-11eb-8f9c-2ea09a1c09d8.png)\n\n---\n\n3. **Spatial Supervision (2)**: Use Cropped-Roi-Mask, blend it with the original image to produce **in-place reflected precise spatial location**. To do that, first we will create a gaussian kernel and then convolves it on the mask image. Then we add the convolved output to the original image to produce final input samples. \n\n```python\nkernel <- cv2.getGaussianKernel(ksize=ksize, sigma=sigma)\nkernel2d <- np.dot(kernel, kernel.T)\nmsk <- convolve2d(msk, kernel2d, mode='same')\nblended <- img + msk\n```\n![new](https://user-images.githubusercontent.com/17668390/122655916-17c50280-d178-11eb-9e30-64039bd035eb.png)\n\n\n---\n\nHere, we will be doing on approach 3, the **spatial supervison** by blending the cropped masks with corresponding x-ray samples to highlighten the precise location. It will be on the training time of course. ","metadata":{}},{"cell_type":"code","source":"import os, math\nimport psutil, random \n\nimport numpy as np \nimport pandas as pd \nimport matplotlib.pyplot as plt\n\nimport cv2; print(cv2.__version__)\nimport tensorflow as tf; print(tf.__version__) ","metadata":{"execution":{"iopub.status.busy":"2021-06-20T13:57:25.773782Z","iopub.execute_input":"2021-06-20T13:57:25.774234Z","iopub.status.idle":"2021-06-20T13:57:27.848071Z","shell.execute_reply.started":"2021-06-20T13:57:25.774140Z","shell.execute_reply":"2021-06-20T13:57:27.846616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"MIXED_PRECISION = True\nXLA_ACCELERATE  = False\n\nGPUS = tf.config.experimental.list_physical_devices('GPU')\nif GPUS:\n    try:\n        for GPU in GPUS:\n            tf.config.experimental.set_memory_growth(GPU, True)\n            logical_gpus = tf.config.experimental.list_logical_devices('GPU')\n            print(len(GPUS), \"Physical GPUs,\", len(logical_gpus), \"Logical GPUs\") \n    except RuntimeError as  RE:\n        print(RE)\n\nif MIXED_PRECISION:\n    policy = tf.keras.mixed_precision.experimental.Policy('mixed_float16')\n    tf.keras.mixed_precision.experimental.set_policy(policy)\n    print('Mixed precision enabled')\n\nif XLA_ACCELERATE:\n    tf.config.optimizer.set_jit(True)\n    print('Accelerated Linear Algebra enabled')","metadata":{"execution":{"iopub.status.busy":"2021-06-20T13:57:27.850215Z","iopub.execute_input":"2021-06-20T13:57:27.850719Z","iopub.status.idle":"2021-06-20T13:57:28.688074Z","shell.execute_reply.started":"2021-06-20T13:57:27.850676Z","shell.execute_reply":"2021-06-20T13:57:28.686721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Preprocess","metadata":{}},{"cell_type":"code","source":"study_df = pd.read_csv('../input/siim-covid19-detection/train_study_level.csv'); print(study_df.shape)\nstudy_df['StudyInstanceUID'] = study_df['id'].apply(lambda x: x.replace('_study', ''))\ndel study_df['id']\n\ndef hot_to_sparse(row):\n    return(row.index[row.apply(lambda x: x==1)][0])\nstudy_df['diagnosis'] = study_df.apply(lambda row:hot_to_sparse(row), axis=1)\ncls = {\n    'Typical Appearance':1,                    \n    'Negative for Pneumonia':2,                \n    'Indeterminate Appearance':3,                     \n    'Atypical Appearance':4,    \n}\nstudy_df['sparse_gt'] = study_df.diagnosis.map(cls) \n\nimage_df = pd.read_csv('../input/siim-covid19-detection/train_image_level.csv'); print(image_df.shape)\ntrain = image_df.merge(study_df, on='StudyInstanceUID')\ntrain['id'] = train['id'].apply(lambda x: x.replace('_image', ''))\ndisplay(train.head()); print(train.shape)","metadata":{"execution":{"iopub.status.busy":"2021-06-20T13:57:28.692084Z","iopub.execute_input":"2021-06-20T13:57:28.692406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ROI Segment: Cropped Bounding Box Mask Blending","metadata":{}},{"cell_type":"code","source":"import scipy.stats as st\nfrom scipy.signal import convolve2d\n\ndef blend(img, msk):\n    '''\n    https://gist.github.com/innat/56786d001048ccd4a68ad41e71869abb\n    '''\n    img = np.asarray(img, np.float32) / 255.\n    msk = np.asarray(msk, np.float32) / 255.\n\n    sigma = np.random.uniform(low=2, high=5)\n    gamma = np.random.uniform(low=1.3, high=1.3)\n    kernel = cv2.getGaussianKernel(ksize=3, sigma=sigma)\n    kernel2d = np.dot(kernel, kernel.T)\n\n    msk[...,0] = convolve2d(msk[...,0], kernel2d, mode='same')\n    blended = img + msk\n\n    if np.max(blended) > 1:\n        m = blended[blended > 1]\n        m = (np.mean(m) - 1) * gamma\n        msk = np.clip(msk - m, 0, 1)\n        blended = np.clip(msk + img, 0, 1)\n    return blended","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def vis(path1, path2, n_images, is_random=True, figsize=(16, 16)):\n    '''\n    https://gist.github.com/innat/00de7561033ba373745d425c6da7bf8c\n    '''\n    image_names = os.listdir(path1)\n    masks_names = os.listdir(path2)\n    \n    for i in range(n_images):\n        if is_random:\n            image_name = random.choice(masks_names)\n            masks_name = image_name\n        else:\n            image_name = masks_names[i]\n            masks_name = masks_names[i]\n        \n        img = cv2.resize(cv2.imread(os.path.join(path1, image_name)), (512, 512))\n        msk = cv2.resize(cv2.imread(os.path.join(path2, masks_name)), (512, 512))\n        bld = blend(img, msk)\n        \n        plt.figure(figsize=(15,15))\n        plt.subplot(131); plt.imshow(img);\n        plt.subplot(132); plt.imshow(msk);\n        plt.subplot(133); plt.imshow(bld);\n        plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_path = '../input/covid19-detection-890pxpng-study'\nTRAIN_IMG_PATH =  os.path.join(base_path, 'train/')\nTRAIN_MSK_PATH = os.path.join(base_path, 'ROI Mask/')\nvis(TRAIN_IMG_PATH, TRAIN_MSK_PATH, 10, is_random=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import GroupKFold\n\nskf = GroupKFold(n_splits=3)\nfor index, (train_index, val_index) in enumerate(skf.split(train, groups = train.id.tolist())):\n    train.loc[val_index, 'fold'] = index\n    \nprint(train.groupby(['fold', train.sparse_gt]).size())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Generator","metadata":{}},{"cell_type":"code","source":"import albumentations as A \n\n# For Validation \ndef albu_transforms_train(data_resize): \n    return A.Compose([\n        A.Resize(data_resize, data_resize),\n        A.ToFloat(), # no need if use keras.applicaiton.EfficientNets Bx\n        A.RandomBrightnessContrast(brightness_limit=0.2,\n                                   contrast_limit=0.2, p=0.75),\n        A.CoarseDropout(max_holes=10, p=0.5),\n        A.ShiftScaleRotate(shift_limit=0.1, scale_limit=0.1,\n                           rotate_limit=15, border_mode=0, p=0.65)\n\n    ], p=1.)\n\n# For Validation \ndef albu_transforms_valid(data_resize): \n    return A.Compose([\n        A.Resize(data_resize, data_resize),\n        A.ToFloat(), # no need if use keras.applicaiton.EfficientNets Bx\n        ], p=1.)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Covid19Generator(tf.keras.utils.Sequence):\n    def __init__(self, img_path, msk_path, data, batch_size, random_state, \n                 dim, shuffle=True, transform=None, is_train=False):\n        self.dim = dim\n        self.data = data\n        self.random_state = random_state\n        self.shuffle  = shuffle\n        \n        self.img_path = img_path\n        self.msk_path = msk_path\n        self.is_train = is_train\n        \n        self.augment  = transform\n        self.batch_size = batch_size\n        \n        self.list_idx = data.index.values\n        self.label = self.data[['Negative for Pneumonia', \n                                'Typical Appearance', \n                                'Indeterminate Appearance', \n                                'Atypical Appearance']] if self.is_train else np.nan\n        self.on_epoch_end()\n        \n    def __len__(self):\n        return int(np.floor(len(self.list_idx) / self.batch_size))\n    \n    def __getitem__(self, index):\n        batch_idx = self.indices[index*self.batch_size:(index+1)*self.batch_size]\n        idx = [self.list_idx[k] for k in batch_idx]\n        \n        Data = np.zeros((self.batch_size,) + self.dim + (1,), dtype=\"float32\")\n        Mask = np.zeros((self.batch_size,) + self.dim + (1,), dtype=\"float32\")\n        Target = np.zeros((self.batch_size, 4), dtype = np.float32)\n\n        for i, k in enumerate(idx):\n            # load the image file using cv2\n            image = cv2.imread(self.img_path + \n                               self.data['id'][k] + '.png', 0) \n            mask = cv2.imread(self.msk_path + \n                              self.data['id'][k] + '.png', 0)\n            \n            if mask is None:\n                # no mask.png for no bounding box image \n                mask = np.zeros_like(image)\n            \n            if self.shuffle: # do blend only for training - not for validation\n                image = self.blend(image[:, :, np.newaxis], \n                                   mask[:, :, np.newaxis]) # WARNING: blending can be expensive\n            else:\n                image = image[:, :, np.newaxis]\n                \n            res = self.augment(image=image)\n            image = res['image']\n            \n            # assign \n            if self.is_train:\n                Data[i,] = image\n                Target[i,] = self.label.iloc[k,].values \n            else:\n                Data[i,] =  image \n        \n        return Data, Target\n    \n    def on_epoch_end(self):\n        self.indices = np.arange(len(self.list_idx))\n        if self.shuffle:\n            np.random.seed(self.random_state)\n            np.random.shuffle(self.indices)\n            \n    def blend(self, img, msk):\n        img = np.asarray(img, np.float32) / 255.\n        msk = np.asarray(msk, np.float32) / 255.\n\n        sigma = np.random.uniform(low=2, high=5)\n        gamma = np.random.uniform(low=1.3, high=1.3)\n        kernel = cv2.getGaussianKernel(ksize=3, sigma=sigma)\n        kernel2d = np.dot(kernel, kernel.T)\n\n        msk[...,0] = convolve2d(msk[...,0], kernel2d, mode='same')\n        blended = img + msk\n        \n        if np.max(blended) > 1:\n            m = blended[blended > 1]\n            m = (np.mean(m) - 1) * gamma\n            msk = np.clip(msk - m, 0, 1)\n            blended = np.clip(msk + img, 0, 1)\n        return blended","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pylab import rcParams\n\n# helper function to plot sample \ndef plot_imgs(dataset_show, row, col):\n    rcParams['figure.figsize'] = 20,10\n    for i in range(row):\n        f, ax = plt.subplots(1,col)\n        for p in range(col):\n            idx = np.random.randint(0, len(dataset_show))\n            img, label = dataset_show[idx]\n            ax[p].grid(False)\n            ax[p].imshow(img[0], cmap='gray')\n            ax[p].set_title(label[0])\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TRAIN_IMG_PATH = '../input/covid19-detection-890pxpng-study/train/'\nTRAIN_MSK_PATH = '../input/covid19-detection-890pxpng-study/ROI Mask/'\n\nbatch_size = 86\nsize = 384\nfold = 1\n\ndef count_data_items(length, b_max):\n    batch_size = sorted([int(length/n) for n in range(1, length+1) \\\n                         if length % n == 0 and length/n <= b_max], reverse=True)[0]  \n    steps  = length / batch_size \n    return batch_size, steps\n\ndef fold_generator(fold):\n    # for way one - data generator\n    train_labels = train[train.fold != fold].reset_index(drop=True)\n    val_labels = train[train.fold == fold].reset_index(drop=True)\n\n    train_generator = Covid19Generator(TRAIN_IMG_PATH, TRAIN_MSK_PATH,\n                              train_labels, \n                              batch_size, 1234, (size, size),\n                              shuffle = True, is_train = True,\n                              transform = albu_transforms_train(size))\n    \n    valid_batch, valid_step = count_data_items(len(val_labels), batch_size)\n\n    val_generator = Covid19Generator(TRAIN_IMG_PATH, TRAIN_MSK_PATH,\n                              val_labels, \n                              valid_batch, 1234, (size, size),\n                              shuffle = False, is_train = True,\n                              transform = albu_transforms_valid(size))\n\n    return train_generator, val_generator, train_labels, val_labels, valid_step\n\ntrain_gen, val_gen, train_len, val_len, val_step = fold_generator(fold)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_imgs(train_gen, 10, 4)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_imgs(val_gen, 2, 4)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model \n","metadata":{}},{"cell_type":"code","source":"from tensorflow.keras import Model \nfrom tensorflow.keras import Sequential \nfrom tensorflow.keras import Input \nfrom tensorflow.keras import layers \nfrom tensorflow.keras import applications \n\nclass CovidNet(Model):\n    '''To handle different channel input for pre-trained weights \n    option 1: https://stackoverflow.com/a/67576025/9215780\n    option 2: https://stackoverflow.com/a/67540516/9215780\n    '''\n    def __init__(self):\n        super(CovidNet, self).__init__()\n        self.base = applications.ResNet50(input_shape=(size, size, 1),\n                                                      include_top=False,\n                                                      weights=None)\n        base_weights = applications.ResNet50(input_shape=(size, size, 3),\n                                                         include_top=False,\n                                                         weights='imagenet')\n\n        for i in range(3, len(self.base.layers)):\n            self.base.layers[i].set_weights(base_weights.layers[i].get_weights())\n\n        del base_weights\n        self.pool = layers.GlobalAveragePooling2D()\n        self.out = layers.Dense(units=4, activation='softmax')\n\n    def call(self, x, training=None, **kwargs):\n        x = self.base(x)\n        x = self.pool(x)\n        x = self.out(x)\n        return x\n\n\ntf.keras.backend.clear_session()\nmodel = CovidNet()\nmodel.build(input_shape=(None, size, size, 1))\nmodel.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training","metadata":{}},{"cell_type":"code","source":"from tensorflow.keras.optimizers.schedules import LearningRateSchedule, ExponentialDecay\n\nclass WarmupLearningRateSchedule(LearningRateSchedule):\n    \"\"\"Provides a variety of learning rate decay schedules with warm up.\"\"\"\n\n    def __init__(self,\n               initial_lr,\n               steps_per_epoch=None,\n               lr_decay_type='exponential',\n               decay_factor=0.97,\n               decay_epochs=2.4,\n               total_steps=None,\n               warmup_epochs=5,\n               minimal_lr=0):\n        super(WarmupLearningRateSchedule, self).__init__()\n        self.initial_lr = initial_lr\n        self.steps_per_epoch = steps_per_epoch\n        self.lr_decay_type = lr_decay_type\n        self.decay_factor = decay_factor\n        self.decay_epochs = decay_epochs\n        self.total_steps = total_steps\n        self.warmup_epochs = warmup_epochs\n        self.minimal_lr = minimal_lr\n\n    def __call__(self, step):\n        if self.lr_decay_type == 'exponential':\n            assert self.steps_per_epoch is not None\n            decay_steps = self.steps_per_epoch * self.decay_epochs\n            lr = ExponentialDecay(self.initial_lr, decay_steps, \n                                  self.decay_factor, staircase=True)(step)\n        elif self.lr_decay_type == 'cosine':\n            assert self.total_steps is not None\n            lr = 0.5 * self.initial_lr * (\n              1 + tf.cos(np.pi * tf.cast(step, tf.float32) / self.total_steps))\n        elif self.lr_decay_type == 'linear':\n            assert self.total_steps is not None\n            lr = (1.0 - tf.cast(step, tf.float32) / self.total_steps) * self.initial_lr\n        elif self.lr_decay_type == 'constant':\n            lr = self.initial_lr\n        else:\n            assert False, 'Unknown lr_decay_type : %s' % self.lr_decay_type\n\n        if self.minimal_lr:\n            lr = tf.math.maximum(lr, self.minimal_lr)\n\n        if self.warmup_epochs:\n            warmup_steps = int(self.warmup_epochs * self.steps_per_epoch)\n            warmup_lr = (\n              self.initial_lr * tf.cast(step, tf.float32) /\n              tf.cast(warmup_steps, tf.float32))\n            lr = tf.cond(step < warmup_steps, lambda: warmup_lr, lambda: lr)\n\n        return lr\n\n    def get_config(self):\n        return {\n            'initial_lr': self.initial_lr,\n            'steps_per_epoch': self.steps_per_epoch,\n            'lr_decay_type': self.lr_decay_type,\n            'decay_factor': self.decay_factor,\n            'decay_epochs': self.decay_epochs,\n            'total_steps': self.total_steps,\n            'warmup_epochs': self.warmup_epochs,\n            'minimal_lr': self.minimal_lr,\n        }","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"steps_per_epoch  = np.ceil(float(len(train_len)) / batch_size) \nvalidation_steps = val_step \nepochs = 10\n\nlr_sched = 'cosine'\nlr_base = 0.016\nlr_min=0\nlr_decay_epoch = 2.4\nlr_warmup_epoch = 5\nlr_decay_factor = 0.97\n\nscaled_lr = lr_base * (batch_size / 256.0)\nscaled_lr_min = lr_min * (batch_size / 256.0)\ntotal_steps = steps_per_epoch * epochs\n\nlearning_rate = WarmupLearningRateSchedule(\n    scaled_lr,\n    steps_per_epoch=steps_per_epoch,\n    decay_epochs=lr_decay_epoch,\n    warmup_epochs=lr_warmup_epoch,\n    decay_factor=lr_decay_factor,\n    lr_decay_type=lr_sched,\n    total_steps=total_steps,\n    minimal_lr=scaled_lr_min)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tensorflow.keras import metrics\nfrom tensorflow.keras import losses \nfrom tensorflow.keras import optimizers\n\n# bind all\nmodel.compile(\n    loss = losses.CategoricalCrossentropy(),\n    metrics = [metrics.SensitivityAtSpecificity(0.4, name='@specificity'), \n               metrics.SpecificityAtSensitivity(0.4, name='@sensitivity'),\n               metrics.AUC(curve='PR', summation_method='interpolation')],\n    optimizer = optimizers.Adam(learning_rate))\n\n\n# list of call backs \nfrom tensorflow.keras import callbacks\ncallback_list = [\n       callbacks.ModelCheckpoint(\n            filepath='model.{epoch:02d}-{val_loss:.4f}.h5', \n            save_freq='epoch', verbose=1, monitor='val_loss', \n            save_weights_only=True, save_best_only=True\n       )         \n]\n\n\n# fitter \nmodel.fit(train_gen, \n          steps_per_epoch=steps_per_epoch,\n          validation_data=val_gen, \n          validation_steps=validation_steps,\n          callbacks=callback_list, \n          workers=psutil.cpu_count(),\n          epochs=epochs)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Resources \n- [Add Attention Mechansim](https://www.kaggle.com/ipythonx/tf-keras-ranzcr-multi-attention-efficientnet): Try to add attention function to build hybrid model. \n- [Out-of-Fold Evaluation](https://www.kaggle.com/ipythonx/optimizing-metrics-out-of-fold-weights-ensemble): Compute the oof each fold and compare the match and non-match prediction in order to emphasize on the weak or minor cases. ","metadata":{}}]}