{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!pip uninstall keras -y\n!pip install git+https://github.com/qubvel/segmentation_models\n!git clone https://github.com/SlinkoIgor/ImageDataAugmentor.git","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"prefix = '/kaggle/input/covid-segmentation/'\n\nimages_radiopedia = np.load(os.path.join(prefix, 'images_radiopedia.npy')).astype(np.float32)\nmasks_radiopedia = np.load(os.path.join(prefix, 'masks_radiopedia.npy')).astype(np.int8)\nimages_medseg = np.load(os.path.join(prefix, 'images_medseg.npy')).astype(np.float32)\nmasks_medseg = np.load(os.path.join(prefix, 'masks_medseg.npy')).astype(np.int8)\n\ntest_images_medseg = np.load(os.path.join(prefix, 'test_images_medseg.npy')).astype(np.float32)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np\n\n\ndef visualize(image_batch, mask_batch=None, pred_batch=None, num_samples=8):\n    num_classes = mask_batch.shape[-1] if mask_batch is not None else 0\n    fix, ax = plt.subplots(num_classes + 1, num_samples, figsize=(num_samples * 2, (num_classes + 1) * 2))\n\n    for i in range(num_samples):\n        ax_image = ax[0, i] if num_classes > 0 else ax[i]\n        ax_image.imshow(image_batch[i,:,:,0], cmap='Greys')\n        ax_image.set_xticks([]) \n        ax_image.set_yticks([])\n        \n        if mask_batch is not None:\n            for j in range(num_classes):\n                if pred_batch is None:\n                    mask_to_show = mask_batch[i,:,:,j]\n                else:\n                    mask_to_show = np.zeros(shape=(*mask_batch.shape[1:-1], 3)) \n                    mask_to_show[..., 0] = pred_batch[i,:,:,j] > 0.5\n                    mask_to_show[..., 1] = mask_batch[i,:,:,j]\n                ax[j + 1, i].imshow(mask_to_show, vmin=0, vmax=1)\n                ax[j + 1, i].set_xticks([]) \n                ax[j + 1, i].set_yticks([]) \n\n    plt.tight_layout()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Images from radiopedia are full CT volumes:\nClass 0 is \"ground glass\"<br>\nClass 1 is \"consolidations\"<br>\nClass 2 is \"lungs other\" – it doesn't mean that it is healthy lungs (you don't need to predict this class)<br>\nClass 3 is \"background\" – not lungs (you don't need to predict this class)<br>"},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize(images_radiopedia[30:], masks_radiopedia[30:])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Images from medseg are individual axial slices:\nClasses are same as in radiopedia part"},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize(images_medseg, masks_medseg)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Test images from medseg are individual axial slices:\nYou should make predictions for class 0 \"ground glass\" and class 1 \"consolidation\""},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize(test_images_medseg)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_hists(images1, images2=None):\n    plt.hist(images1.ravel(), bins=100, density=True, color='b', alpha=1 if images2 is None else 0.5)\n    if images2 is not None:\n        plt.hist(images2.ravel(), bins=100, density=True, alpha=0.5, color='orange')\n    plt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Plot images hists:\nHU (Hounsfield scale) of radiopedia data (blue) vs medseg data (orange):"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_hists(images_radiopedia, images_medseg)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"HU (Hounsfield scale) of test medseg data (blue) vs medseg data (orange):"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_hists(test_images_medseg, images_medseg)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Preprocess images:"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"def preprocess_images(images_arr, mean_std=None):\n    images_arr[images_arr > 500] = 500\n    images_arr[images_arr < -1500] = -1500\n    min_perc, max_perc = np.percentile(images_arr, 5), np.percentile(images_arr, 95)\n    images_arr_valid = images_arr[(images_arr > min_perc) & (images_arr < max_perc)]\n    mean, std = (images_arr_valid.mean(), images_arr_valid.std()) if mean_std is None else mean_std\n    images_arr = (images_arr - mean) / std\n    print(f'mean {mean}, std {std}')\n    return images_arr, (mean, std)\n\nimages_radiopedia, mean_std = preprocess_images(images_radiopedia)\nimages_medseg, _ = preprocess_images(images_medseg, mean_std)\ntest_images_medseg, _ = preprocess_images(test_images_medseg, mean_std)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Normalized values of radiopedia data (blue) vs medseg data (orange):"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_hists(images_radiopedia, images_medseg)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Normalized values of test medseg data (blue) vs medseg data (orange):"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_hists(test_images_medseg, images_medseg)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Split train / val:"},{"metadata":{"trusted":true},"cell_type":"code","source":"val_indexes, train_indexes = list(range(24)), list(range(24, 100))\n\ntrain_images = np.concatenate((images_medseg[train_indexes], images_radiopedia))\ntrain_masks = np.concatenate((masks_medseg[train_indexes], masks_radiopedia))\nval_images = images_medseg[val_indexes]\nval_masks = masks_medseg[val_indexes]\n\nbatch_size = len(val_masks)\n\ndel images_radiopedia\ndel masks_radiopedia\ndel images_medseg\ndel masks_medseg","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data generator and augmentations:"},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow\n\nimport albumentations\nimport cv2\n\n#for PSPNet\nSOURCE_SIZE = 510\nTARGET_SIZE = 384\n\ntrain_augs = albumentations.Compose([\n    albumentations.Rotate(limit=360, p=0.9, border_mode=cv2.BORDER_REPLICATE),\n    albumentations.RandomSizedCrop((int(SOURCE_SIZE * 0.75), SOURCE_SIZE), \n                                   TARGET_SIZE, \n                                   TARGET_SIZE, \n                                   interpolation=cv2.INTER_NEAREST),\n    albumentations.HorizontalFlip(p=0.5),\n\n])\n\nval_augs = albumentations.Compose([\n    albumentations.Resize(TARGET_SIZE, TARGET_SIZE, interpolation=cv2.INTER_NEAREST)\n])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class Dataset:   \n    def __init__(\n            self, \n            images, \n            masks,\n            augmentations=None\n    ):\n        self.images = images\n        self.masks = masks\n        self.augmentations = augmentations\n    \n    def __getitem__(self, i):\n        image = self.images[i]\n        mask = self.masks[i]\n        \n        if self.augmentations:\n            sample = self.augmentations(image=image, mask=mask)\n            image, mask = sample['image'], sample['mask']\n        return image, mask\n        \n    def __len__(self):\n        return len(self.images)\n    \n    \nclass Dataloder(tensorflow.keras.utils.Sequence):\n    \"\"\"Load data from dataset and form batches\n    \n    Args:\n        dataset: instance of Dataset class for image loading and preprocessing.\n        batch_size: Integet number of images in batch.\n        shuffle: Boolean, if `True` shuffle image indexes each epoch.\n    \"\"\"\n    \n    def __init__(self, dataset, batch_size=1, shuffle=False):\n        self.dataset = dataset\n        self.batch_size = batch_size\n        self.shuffle = shuffle\n        self.indexes = np.arange(len(dataset))\n\n        self.on_epoch_end()\n\n    def __getitem__(self, i):\n        \n        # collect batch data\n        start = i * self.batch_size\n        stop = (i + 1) * self.batch_size\n        images = []\n        masks = []\n        for j in range(start, stop):\n            image, mask = self.dataset[self.indexes[j]]\n            images.append(image)\n            masks.append(mask)\n        \n        images = np.stack(images, axis=0)\n        masks = np.stack(masks, axis=0).astype(np.float32)\n        \n        return (images, masks)\n    \n    def __len__(self):\n        \"\"\"Denotes the number of batches per epoch\"\"\"\n        return len(self.indexes) // self.batch_size\n    \n    def on_epoch_end(self):\n        \"\"\"Callback function to shuffle indexes each epoch\"\"\"\n        if self.shuffle:\n            self.indexes = np.random.permutation(self.indexes)\n            \ntrain_dataset = Dataset(train_images, train_masks, train_augs)\nval_dataset = Dataset(val_images, val_masks, val_augs)\n\ntrain_dataloader = Dataloder(train_dataset, batch_size=batch_size, shuffle=True)\nval_dataloader = Dataloder(val_dataset, batch_size=batch_size, shuffle=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"assert train_dataloader[0][0].shape == (batch_size, TARGET_SIZE, TARGET_SIZE, 1)\nassert train_dataloader[0][1].shape == (batch_size, TARGET_SIZE, TARGET_SIZE, 4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize(*next(iter(train_dataloader)))\nvisualize(*next(iter(val_dataloader)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Metrics:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def fscore_glass(y_true, y_pred):\n    return sm.metrics.f1_score(y_true[..., 0:1], \n                               y_pred[..., 0:1])\n    \ndef fscore_consolidation(y_true, y_pred):\n    return sm.metrics.f1_score(y_true[..., 1:2], \n                               y_pred[..., 1:2])\n\ndef fscore_lungs_other(y_true, y_pred):\n    return sm.metrics.f1_score(y_true[..., 2:3], \n                               y_pred[..., 2:3])\n\ndef fscore_glass_and_consolidation(y_true, y_pred):\n    return sm.metrics.f1_score(y_true[..., :2], \n                               y_pred[..., :2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nfrom segmentation_models import PSPNet\nimport segmentation_models as sm\n\nfrom tensorflow.keras.layers import Input, Conv2D\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.callbacks import ModelCheckpoint\nnp.random.seed(0)\n\nbase_model = PSPNet(backbone_name='efficientnetb0',\n                  encoder_weights='imagenet',\n                  classes=4, \n                  activation='softmax')\n\nmodel = Sequential([Input(shape=(TARGET_SIZE, TARGET_SIZE, 1)),\n                    Conv2D(3, (1, 1)),  # map N channels data to 3 channels\n                    base_model])\n    \nmodel.compile(Adam(learning_rate=0.001, amsgrad=True),\n              loss=sm.losses.categorical_crossentropy,\n              metrics=[fscore_glass, fscore_consolidation, fscore_lungs_other, fscore_glass_and_consolidation])\n\ncheckpoint_callback = ModelCheckpoint('best_model',\n                                      monitor='fscore_glass_and_consolidation',\n                                      mode='max',\n                                      save_best_only=True)\n\nmodel.fit(\n    train_dataloader,\n    steps_per_epoch=len(train_dataloader) * 6,\n    epochs=10,\n    validation_data=val_dataloader,\n    validation_steps=len(val_dataloader),\n    callbacks=[checkpoint_callback],\n    workers=4)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load best model and visualize predicions on val:"},{"metadata":{"trusted":true},"cell_type":"code","source":"del train_images\ndel train_masks","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = tensorflow.keras.models.load_model('best_model/',\n                                           compile=False,\n                                           custom_objects={\n                                                'categorical_crossentropy': sm.losses.categorical_crossentropy,\n                                                'fscore_consolidation': fscore_consolidation, \n                                                'fscore_glass': fscore_glass, \n                                                'fscore_lungs_other': fscore_lungs_other,\n                                                'fscore_glass_and_consolidation': fscore_glass_and_consolidation})\n\nmodel.compile(Adam(learning_rate=0.001, amsgrad=True),\n              loss=sm.losses.jaccard_loss)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input = val_dataloader[0]\nimage_batch, mask_batch = input\n\npreds = model.predict_on_batch(image_batch)\nvisualize(image_batch, mask_batch, pred_batch=preds)\n\n# yellow is TP, red is FP, green is FN.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Test preds:"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_batch = np.stack([val_augs(image=img)['image'] for img in test_images_medseg], axis=0)\ntest_preds = model.predict_on_batch(image_batch)\ntest_masks_prediction = test_preds > 0.5\nvisualize(image_batch, test_masks_prediction, num_samples=len(test_images_medseg))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Resize prediction to original size:"},{"metadata":{"trusted":true},"cell_type":"code","source":"import scipy\ntest_masks_prediction_original_size = scipy.ndimage.zoom(test_masks_prediction[..., :-2], (1, 2, 2, 1), order=0)\ntest_masks_prediction_original_size.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\n\npd.DataFrame(\n             data=np.stack((np.arange(len(test_masks_prediction_original_size.ravel())), \n                            test_masks_prediction_original_size.ravel().astype(int)),\n                            axis=-1), \n             columns=['Id', 'Predicted'])\\\n.set_index('Id').to_csv('submission.csv')","execution_count":null,"outputs":[]}],"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}