{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":20561,"databundleVersionId":1187889,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true},"colab":{"provenance":[]}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\n# IMPORTANT: RUN THIS CELL IN ORDER TO IMPORT YOUR KAGGLE DATA SOURCES\n# TO THE CORRECT LOCATION (/kaggle/input) IN YOUR NOTEBOOK,\n# THEN FEEL FREE TO DELETE THIS CELL.\n# NOTE: THIS NOTEBOOK ENVIRONMENT DIFFERS FROM KAGGLE'S PYTHON\n# ENVIRONMENT SO THERE MAY BE MISSING LIBRARIES USED BY YOUR\n# NOTEBOOK.\n\nimport os\nimport sys\nfrom tempfile import NamedTemporaryFile\nfrom urllib.request import urlopen\nfrom urllib.parse import unquote, urlparse\nfrom urllib.error import HTTPError\nfrom zipfile import ZipFile\nimport tarfile\nimport shutil\n\nCHUNK_SIZE = 40960\nDATA_SOURCE_MAPPING = 'covid-segmentation:https%3A%2F%2Fstorage.googleapis.com%2Fkaggle-competitions-data%2Fkaggle-v2%2F20561%2F1187889%2Fbundle%2Farchive.zip%3FX-Goog-Algorithm%3DGOOG4-RSA-SHA256%26X-Goog-Credential%3Dgcp-kaggle-com%2540kaggle-161607.iam.gserviceaccount.com%252F20240906%252Fauto%252Fstorage%252Fgoog4_request%26X-Goog-Date%3D20240906T172045Z%26X-Goog-Expires%3D259200%26X-Goog-SignedHeaders%3Dhost%26X-Goog-Signature%3D43742a5ce41d83ccfc3c8743c5c8079abe9d6551dac4966a23f0bc0a359986061563265bf12ebfe02c15b74ea516fd390ae87498e2b0eeee4d0c4130ff222cb97c8180aec7776cdaa156a56afbb9f3bf08925e57e8ef81a3037b39d68d6aaff017bb0bbe4ccf290b729140201790b7e8858beec45fb63e06ff0e01a0104b2d197c79b64fbc94fc2b8466da932dfe6f4d07651cb02bf789876345aa01ebdc11e32c56c168e4cbc2cd8fd289562cc0e8657c10e99c48d53bf333bb326e02f09e13f8e67438790fc2787521404e3b6ce197637c4398bd5dd04f4f099e48afc7a8b2e837ea63c9aa7ff785a50612b2fadb12bdaca4ed058122abc50b1e866742e2bf'\n\nKAGGLE_INPUT_PATH='/kaggle/input'\nKAGGLE_WORKING_PATH='/kaggle/working'\nKAGGLE_SYMLINK='kaggle'\n\n!umount /kaggle/input/ 2> /dev/null\nshutil.rmtree('/kaggle/input', ignore_errors=True)\nos.makedirs(KAGGLE_INPUT_PATH, 0o777, exist_ok=True)\nos.makedirs(KAGGLE_WORKING_PATH, 0o777, exist_ok=True)\n\ntry:\n os.symlink(KAGGLE_INPUT_PATH, os.path.join(\"..\", 'input'), target_is_directory=True)\nexcept FileExistsError:\n pass\ntry:\n os.symlink(KAGGLE_WORKING_PATH, os.path.join(\"..\", 'working'), target_is_directory=True)\nexcept FileExistsError:\n pass\n\nfor data_source_mapping in DATA_SOURCE_MAPPING.split(','):\n directory, download_url_encoded = data_source_mapping.split(':')\n download_url = unquote(download_url_encoded)\n filename = urlparse(download_url).path\n destination_path = os.path.join(KAGGLE_INPUT_PATH, directory)\n try:\n with urlopen(download_url) as fileres, NamedTemporaryFile() as tfile:\n total_length = fileres.headers['content-length']\n print(f'Downloading {directory}, {total_length} bytes compressed')\n dl = 0\n data = fileres.read(CHUNK_SIZE)\n while len(data) > 0:\n dl += len(data)\n tfile.write(data)\n done = int(50 * dl / int(total_length))\n sys.stdout.write(f\"\\r[{'=' * done}{' ' * (50-done)}] {dl} bytes downloaded\")\n sys.stdout.flush()\n data = fileres.read(CHUNK_SIZE)\n if filename.endswith('.zip'):\n with ZipFile(tfile) as zfile:\n zfile.extractall(destination_path)\n else:\n with tarfile.open(tfile.name) as tarfile:\n tarfile.extractall(destination_path)\n print(f'\\nDownloaded and uncompressed: {directory}')\n except HTTPError as e:\n print(f'Failed to load (likely expired) {download_url} to path {destination_path}')\n continue\n except OSError as e:\n print(f'Failed to load {download_url} to path {destination_path}')\n continue\n\nprint('Data source import complete.')\n","metadata":{"id":"DOaGPSVSWlEK","outputId":"e826e14d-213d-43f9-bf1c-bdf091d24b56"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# COVID-19 CT Images Segmentation","metadata":{"id":"T6P7wlUcWlEO"}},{"cell_type":"markdown","source":"This solution is based on [this notebook](https://www.kaggle.com/code/maedemaftouni/pytorch-baseline-for-semantic-segmentation) by Maede Maftouni.","metadata":{"id":"EfpDLzq9SNIk"}},{"cell_type":"markdown","source":"Team 6:\n* A01275287 Oskar Adolfo Villa López\n* A01736345 Luis Ángel Cruz García\n* A01369230 Dafne Fernández Hernández\n* A01368818 Joel Sánchez Olvera\n\n\n","metadata":{"id":"cMUjBt5yWlEQ"}},{"cell_type":"markdown","source":"# **Introduction**\n\nThis notebook addresses a critical problem in medical image segmentation related to COVID-19, aiming to identify and delineate specific lung areas in computed tomography (CT) scans. The segmentation problem focuses on detecting key pulmonary features, such as ground glass opacities and consolidations, which are essential indicators of COVID-19 infection. Accurate segmentation of these areas is crucial for medical diagnoses, as it helps assess the severity of the infection and provides valuable information for disease monitoring and treatment planning.\n\n**The Problem:**\n\nCOVID-19 primarily affects the lungs, and CT scans can show telltale signs of infection, like ground glass opacities and consolidations, which indicate inflammation or fluid buildup in the lungs. Accurately segmenting these regions in CT images is critical for diagnosing the severity of the infection and monitoring disease progression. This segmentation task aids in providing visual evidence of affected lung regions, improving clinical diagnosis, and supporting medical treatment decisions.\n\n**Kaggle Competition Objectives**\n\nThe Kaggle competition's objective is to develop robust deep learning models capable of accurately segmenting lung regions affected by COVID-19 in CT scans. The relevance of this task is significant in the medical field, as it enhances diagnostic tools and allows healthcare professionals to track disease progression. By providing precise visual information about the extent of lung involvement, these models assist in making informed clinical decisions and improving patient care.\n\nThe notebook employs popular segmentation models based on architectures like U-Net and FPN (Feature Pyramid Network), widely used in medical imaging tasks due to their ability to capture both fine and global contextual information. The notebook also makes use of the segmentation_models_pytorch library, which simplifies model building and training by providing pre-trained segmentation backbones.\n\n# **Data Analysis**\n\n**Exploratory Data Analysis (EDA)**\n\nThe dataset provided in this competition includes lung CT scans with corresponding segmentation masks. The CT images come from two datasets:\n\n* **Medseg Dataset:**\n \n 100 CT images from over 10 patients.\n\n Corresponding segmentation masks contain four channels:\n \n * Ground glass opacities\n * Consolidations\n * Unaffected lung tissue (Lungs Other)\n * Background.\n\n\n* **Radiopaedia Dataset:**\n\n 829 CT slices.\n\n The structure of the images and masks is similar to the Medseg dataset.\n\nEach image has dimensions of 520 x 520 pixels. The segmentation masks label regions of interest in the lungs, including infected areas (ground glass and consolidation) and normal lung tissue. Visualizing these images helps understand the data distribution and the extent of the infection in the lungs.\n\n**Data Preprocessing**\n\nBefore training the model, several data preprocessing steps are performed:\n\n* **Resizing:** All images are resized to 520 x 520 pixels to ensure consistency in input dimensions.\n* **Normalization:** Pixel values are normalized to a standard range, improving model convergence during training.\n* **Data Augmentation:** Techniques such as random rotations, flips, and intensity shifts are applied to increase the dataset's variability, helping the model generalize better to unseen data.\n\nThese steps ensure the model can handle different variations in the input data, making it more robust.\n\n# **Model**\n\n**Model Architecture**\n\nThe model chosen for this task is based on the U-Net architecture, a popular choice for medical image segmentation due to its symmetrical design and skip connections that retain spatial information.\n\nThe U-Net architecture is well-suited for this task because of its ability to accurately segment complex shapes in medical images, such as the lung regions affected by COVID-19. The use of skip connections ensures that the model does not lose important spatial information during downsampling. This makes U-Net an optimal choice for the challenge of identifying intricate lung structures.\n\nThe architecture consists of:\n\n* **Encoder:** A series of convolutional layers followed by max-pooling layers to downsample the input image and extract features.\n* **Decoder:** Upsampling layers that reconstruct the image to its original resolution while merging feature maps from the encoder via skip connections to retain fine details.\n* **Skip Connections:** Allow the model to merge high-level features with low-level spatial information, improving segmentation accuracy for small details.\n\nOther model details:\n\n* **Kernel Size:** The convolutional layers use a standard kernel size of 3x3, with stride 1 and appropriate padding to maintain spatial dimensions.\nMax Pooling: Reduces the spatial dimensions of the feature maps while retaining the most important features.\n* **Activation Functions:** The model uses ReLU for non-linearity in the convolutional layers and Sigmoid in the output layer for binary segmentation.\n* **Cross-Entropy Loss:**\nThe cross-entropy loss function is used in this model. It measures the difference between the predicted segmentation mask and the ground truth mask by penalizing incorrect classifications for each pixel. Cross-entropy loss is widely used in segmentation tasks due to its effectiveness in dealing with multiple classes and imbalanced datasets.\n\n \n\n\n\n* Where:\n\n 𝑦 is the true label for pixel 𝑖 , and 𝑝 is the predicted probability for that label.\n\n\n\n\n\n**References**\n\nRonneberger, O., Fischer, P., & Brox, T. (2015). U-Net: Convolutional Networks for Biomedical Image Segmentation. International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI).","metadata":{"id":"x3mGMcnvWlES"}},{"cell_type":"markdown","source":"## Ethical implications\nThe analyzed images were originally obtained from SIRM (Società Italiana di Radiologia Medica e Interventistica), which follows a series of regulations to protect individuals' information and allow its use for academic or scientific purposes. Among these regulations, we can find:\n\n1. D. Lgs. 101/2018: Implements the General Data Protection Regulation (GDPR) in Italy. This is crucial as it sets the rules for processing personal data, including medical images that may contain sensitive information. In this specific case, the following measures must be followed:\n * Anonymization of the images to ensure that patients cannot be identified.\n * Explicit consent from patients to use their images.\n\n2. Legge 219/2017: This law establishes the rules for informed consent and advance treatment directives.\n\n3. Legge 3/2018: Covers aspects of clinical experimentation and the reorganization of healthcare professions.\n\nAlso, it's important to emphasize that as a team, we focus not only on finding solutions that enhance our learning, but also on creating real-world applications that can help others. This, in turn, gives us the experience needed to develop more efficient solutions.\n\n\n**References**\n\n\nSocietà Italiana di Radiologia Medica e Interventistica (2024). Archivio Documenti. Normativa. https://sirm.org/documenti/\n\n","metadata":{"id":"FVTYdvxBco_K"}},{"cell_type":"markdown","source":"## Imports","metadata":{"id":"itocfzMnWlET"}},{"cell_type":"code","source":"!pip install segmentation_models_pytorch","metadata":{"id":"auyuVE6PvIKu","outputId":"7a959cc4-a603-438a-c034-f9f08c824ec8"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\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 20GB 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","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-09-04T16:20:22.450771Z","iopub.execute_input":"2024-09-04T16:20:22.451233Z","iopub.status.idle":"2024-09-04T16:20:22.468035Z","shell.execute_reply.started":"2024-09-04T16:20:22.451186Z","shell.execute_reply":"2024-09-04T16:20:22.466453Z"},"id":"PFjJFcKyWlEU","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data description","metadata":{"id":"Q9p05rGm-gCf"}},{"cell_type":"markdown","source":"The data used for training comes from two diffrent datasets: Medseg and Radiopedia. All images and mask have a size of 520 x 520. Examples of the images are provided in the **Extraction** section.\n\n\n---\n\n\nThe Medseg set is made from 100 CT images from more than 40 patients. The set contains the 100 original images and 100 masks with four channels (0 - \"ground glass\", 1 - \"consolidations\", 2 - \"lungs other\", 3 - \"background\" ). The test set is also taken from this set, and it consists of 10 images.\n\n---\n\n\nThe Radiopedia set follows the same structure of the images and masks. It contains 829 slices from segmented axial volumetric CTs.","metadata":{"id":"AVhUTJKU_CAe"}},{"cell_type":"markdown","source":"## Extraction","metadata":{"id":"kVsd-yprWlEV"}},{"cell_type":"markdown","source":"The data is extracted from the given repositories. The images are stored in numpy arrays. The images are already formatted and have their masks assigned as labels.","metadata":{"id":"a_lxnUNOWlEW"}},{"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)","metadata":{"execution":{"iopub.status.busy":"2024-09-04T16:20:53.45137Z","iopub.execute_input":"2024-09-04T16:20:53.451904Z","iopub.status.idle":"2024-09-04T16:20:55.849166Z","shell.execute_reply.started":"2024-09-04T16:20:53.45185Z","shell.execute_reply":"2024-09-04T16:20:55.847956Z"},"id":"-kJ61QTUWlEX","outputId":"f3317d9e-d880-4090-daaa-7ac13d3699e2","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A function for visalization is defined. The code prints a matriz with the original images in a grey scale in the superior row. The masks are shown in the correspondent column. The masks that are going to be predicted with the model are ground glass and consolidation, the second and third rows respectively.","metadata":{"id":"Qq4UoXheWlEY"}},{"cell_type":"markdown","source":"**Images from radiopedia are full CT volumes:**
\nClass 0 is \"ground glass\"
\nClass 1 is \"consolidations\"
\nClass 2 is \"lungs other\" – it doesn't mean that it is healthy lungs (not needed for analysis)
\nClass 3 is \"background\" – not lungs (not needed for analysis)
","metadata":{"id":"GiP5dP1oWlEZ"}},{"cell_type":"code","source":"def visualize(image_batch, mask_batch=None, pred_batch=None, num_samples=8, hot_encode=True):\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 if hot_encode:\n ax_image.imshow(image_batch[i,:,:,0], cmap='Greys')\n else:\n ax_image.imshow(image_batch[i,:,:])\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].cpu().numpy() > 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()","metadata":{"execution":{"iopub.status.busy":"2024-09-04T16:21:00.952079Z","iopub.execute_input":"2024-09-04T16:21:00.953333Z","iopub.status.idle":"2024-09-04T16:21:00.966102Z","shell.execute_reply.started":"2024-09-04T16:21:00.953278Z","shell.execute_reply":"2024-09-04T16:21:00.964366Z"},"id":"IgLofHlVWlEa","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize(images_radiopedia[30:], masks_radiopedia[30:])","metadata":{"execution":{"iopub.status.busy":"2024-09-04T16:24:02.251604Z","iopub.execute_input":"2024-09-04T16:24:02.252916Z","iopub.status.idle":"2024-09-04T16:24:05.374474Z","shell.execute_reply.started":"2024-09-04T16:24:02.252833Z","shell.execute_reply":"2024-09-04T16:24:05.373323Z"},"id":"N0rp0CeTWlEb","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, both classes that we want to predict are pretty small compared to the complete image, so we have an important class imbalance.","metadata":{"id":"SFN9wS3JWlEc"}},{"cell_type":"markdown","source":"Then, we define a function to convert the one hot encoded mask to a segmentation mask.","metadata":{"id":"gOwvHtR_vO_q"}},{"cell_type":"code","source":"def onehot_to_mask(mask, palette):\n \"\"\"\n Converts a mask (H, W, K) to (H, W, C)\n \"\"\"\n x = np.argmax(mask, axis=-1)\n colour_codes = np.array(palette)\n x = np.uint8(colour_codes[x.astype(np.uint8)])\n return x\n\npalette = [[0], [1], [2],[3]]\nmasks_radiopedia_recover = onehot_to_mask(masks_radiopedia, palette).squeeze() # shape = (H, W)\n\nmasks_medseg_recover = onehot_to_mask(masks_medseg, palette).squeeze() # shape = (H, W)\n\nprint('Hot encoded mask size: ',masks_radiopedia.shape)\nprint('Paletted mask size:',masks_medseg_recover.shape)\n\nvisualize(masks_medseg_recover[30:],hot_encode=False)","metadata":{"id":"k5lJaBlyvPiN"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we need to separate the image sets into validation and train. The test set is already provided by the competition. 24 images were selected from the medseg set for validation. The rest of medseg and all of the images from radiopedia were selected for training, with a size of ~800.","metadata":{"id":"jEnaSc7FvU5V"}},{"cell_type":"code","source":"masks_radiopedia_recover = onehot_to_mask(masks_radiopedia, palette).squeeze() # shape = (H, W)\nmasks_medseg_recover = onehot_to_mask(masks_medseg, palette).squeeze() # shape = (H, W)\n\n\nval_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_recover[train_indexes], masks_radiopedia_recover))\nval_images = images_medseg[val_indexes]\nval_masks = masks_medseg_recover[val_indexes]\n\nbatch_size = len(val_masks)\n\ndel masks_medseg_recover\ndel masks_radiopedia_recover\ndel images_radiopedia\ndel masks_radiopedia\ndel images_medseg\ndel masks_medseg","metadata":{"id":"lTGT1X6gvVu4"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Transform","metadata":{"id":"deTecJpwWlEc"}},{"cell_type":"markdown","source":"The training data set is relatively small, with less than 1000 images. To improve the variety of the set, we can apply random transformations to the images, which can reduce overfit and increase the accuracy of the model. The library Albumentations is used for this, because it can apply the same transformation to an image and its corresponding mask. The transformations applied are the following:\n* Rotation. With a limit of 360 degrees, which means that the images can be completely rotated. A probability of 90% is given.\n* Random sized crop. Crops a portion of the images with a limit of 75% of the original size, and rescales the cropped image to the original size.\n* Horizontal flip. With a probability of 50%.","metadata":{"id":"svKM8Fwovd_L"}},{"cell_type":"code","source":"import albumentations\n\nimport cv2\n\nSOURCE_SIZE = 512\nTARGET_SIZE = 256\n\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])","metadata":{"id":"dRJCJqA8vdyT"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load","metadata":{"id":"AmL4GrtxWlEd"}},{"cell_type":"markdown","source":"Now, a dataset object is defined. This object has three methods besides its constructor:\n* \\_\\_getitem\\_\\_. Returns an item from the dataset after transforming the image to a tensor. It also applies the tranformation defined previously.\n* \\_\\_len\\_\\_. Returns the length of the dataset.\n* tiles. Returns a two tensors containing an image an its mask, where each item is a patch of the original image.","metadata":{"id":"Mi-ACbMOvhzn"}},{"cell_type":"code","source":"from PIL import Image\nimport segmentation_models_pytorch as smp\nimport torch.nn.functional as F\nimport torch.nn as nn\nfrom tqdm.notebook import tqdm\nimport time\n\nfrom torchvision import transforms as T\nfrom torch.utils.data import Dataset, DataLoader\nimport torch\n\nclass 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 self.mean = [0.485]\n self.std = [0.229]\n\n def __getitem__(self, i):\n image = self.images[i]\n mask = self.masks[i]\n\n\n if self.augmentations is not None:\n sample = self.augmentations(image=image, mask=mask)\n\n image, mask = Image.fromarray(np.squeeze(sample['image'], axis=2)), sample['mask']\n\n if self.augmentations is None:\n image = Image.fromarray(image)\n\n t = T.Compose([T.ToTensor(), T.Normalize(self.mean, self.std)])\n image = t(image)\n mask = torch.from_numpy(mask).long()\n\n return image, mask\n\n def __len__(self):\n return len(self.images)\n\n def tiles(self, image, mask):\n\n img_patches = image.unfold(1, 512, 512).unfold(2, 768, 768)\n img_patches = img_patches.contiguous().view(3,-1, 512, 768)\n img_patches = img_patches.permute(1,0,2,3)\n\n mask_patches = mask.unfold(0, 512, 512).unfold(1, 768, 768)\n mask_patches = mask_patches.contiguous().view(-1, 512, 768)\n\n return img_patches, mask_patches\n\n\ntrain_dataset = Dataset(train_images, train_masks, train_augs)\nval_dataset = Dataset(val_images, val_masks, val_augs)\n\ntrain_dataloader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)\nval_dataloader = DataLoader(val_dataset, batch_size=batch_size, shuffle=False)","metadata":{"id":"38ry36BovjDG"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model","metadata":{"id":"fNBApoPhWlEd"}},{"cell_type":"markdown","source":"Before defining the model, we need some helper functions for visualization and analysis.","metadata":{"id":"62-vijj9v7ge"}},{"cell_type":"markdown","source":"**mask_to_onehot**\n\nWe will convert a segmentation mask into a one-hot encoded representation, which is a binary representation where a class is represented with a vector that has a value of 1 at the position corresponding to that class, and 0 at all other positions. We will follow the next steps in the mask_to_onehot function:\n\n1. An empty list called semantic_map is created that will serve to store the class maps for each color in the palette.\n2. For each color in the palette the mask is compared to the color using. This generates a Boolean matrix indicating where in the image pixels match that color. Then, it is checked if all the color channels match, generating a binary map (H, W) where True indicates that that pixel belongs to the class corresponding to the color of the palette.\nThis class map is added to the semantic_map list.\n3. The class maps for each color are stacked along a new dimension. This converts semantic_map to a matrix of dimensions (H, W, K), where K is the number of classes. Finally, the result is converted to a float tensor and then to a PyTorch tensor.\n4. The function returns a one-hot encoding tensor of dimensions (H, W, K), where each pixel is represented by a one-hot vector indicating which class it belongs to.","metadata":{"id":"JuiHVdSTwISY"}},{"cell_type":"code","source":"def mask_to_onehot(mask, palette):\n \"\"\"\n Converts a segmentation mask (H, W, C) to (H, W, K) where the last dim is a one\n hot encoding vector, C is usually 1 or 3, and K is the number of class.\n \"\"\"\n semantic_map = []\n for colour in palette:\n #print('colour',colour)\n equality = np.equal(mask, colour)\n #print('equality',equality)\n class_map = np.all(equality, axis=-1)\n semantic_map.append(class_map)\n semantic_map = np.stack(semantic_map, axis=-1).astype(np.float32)\n return torch.from_numpy(semantic_map)\n\ni,train_data = next(enumerate(train_dataloader))\n\nmask_hot_encoded = mask_to_onehot(torch.unsqueeze(train_data[1],-1).numpy(),palette)\n#visualize(torch.unsqueeze(torch.squeeze(train_data[0],1),-1),mask_hot_encoded)\nvisualize(train_data[0].permute(0, 2, 3,1),mask_hot_encoded)\n","metadata":{"id":"TUgWHJz9wFQv"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**pixel_accuracy**\n\nWe will compute the pixel-wise accuracy for a predicted segmentation mask by comparing it to the ground truth mask.\nPixel-wise accuracy is a metric used in image segmentation tasks, where the goal is to classify each pixel of an image into a specific category (e.g., sky, tree, road). This metric measures the percentage of pixels in an image that have been correctly classified, that is, those that match the correct label in the reference mask (ground truth).\n\nThe steps followed in the pixel_accuracy function are:\n\n1. The model's output is passed through a softmax activation function to convert it into probabilities for each class at each pixel location.\n2. The class with the highest probability is selected for each pixel using torch.argmax, creating a predicted segmentation mask.\n3. A comparison is made between the predicted mask and the ground truth mask (both now containing class indices), generating a binary matrix where 1 indicates correct predictions and 0 indicates incorrect ones.\n4. The total number of correct predictions is summed, and the result is divided by the total number of pixels in the mask to calculate the pixel-wise accuracy.\n5. The function returns the accuracy as a floating-point value.","metadata":{"id":"KmYrukY2wKMt"}},{"cell_type":"code","source":"def pixel_accuracy(output, mask):\n with torch.no_grad():\n output = torch.argmax(F.softmax(output, dim=1), dim=1)\n correct = torch.eq(output, mask).int()\n accuracy = float(correct.sum()) / float(correct.numel())\n return accuracy","metadata":{"id":"ZVFWdyarwNUB"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**mIoU**\n\nMean Intersection over Union (mIoU) is a metric used in image segmentation tasks to evaluate a model's performance. mIoU calculates the average Intersection over Union (IoU) between the predicted and true labels for each class. **IoU** for a specific class is calculated by dividing the number of pixels correctly predicted as that class (intersection) by the total number of pixels belonging to that class in either the prediction or the ground truth (union).\n\nThe steps followed in the function are:\n\n1. The model's output is passed through a softmax function to obtain probabilities for each class at each pixel.\n2. The class with the highest probability for each pixel is selected using torch.argmax, creating a predicted mask.\n3. Both the predicted mask and the ground truth mask are flattened into one-dimensional vectors for easier comparison.\n4. A list is initialized to store the IoU for each class.\n5. The function loops through each class:\n - It identifies the pixels in the predicted mask that belong to the current class and the corresponding pixels in the ground truth.\n - If there are no pixels of the current class in the ground truth, np.nan is added to the IoU list for that class.\n - Otherwise, the intersection (pixels correctly predicted as that class) and the union (pixels belonging to that class in either mask) are calculated.\n - The IoU for that class is computed by dividing the intersection by the union, and the result is added to the list.\n6. Finally, the function returns the average IoU across all classes, ignoring any nan values, which gives the overall mIoU score.","metadata":{"id":"usfeuLXwwQKC"}},{"cell_type":"code","source":"def mIoU(pred_mask, mask, smooth=1e-10, n_classes=4):\n with torch.no_grad():\n pred_mask = F.softmax(pred_mask, dim=1)\n pred_mask = torch.argmax(pred_mask, dim=1)\n pred_mask = pred_mask.contiguous().view(-1)\n mask = mask.contiguous().view(-1)\n\n iou_per_class = []\n for clas in range(0, n_classes): #loop per pixel class\n true_class = pred_mask == clas\n true_label = mask == clas\n\n if true_label.long().sum().item() == 0: #no exist label in this loop\n iou_per_class.append(np.nan)\n else:\n intersect = torch.logical_and(true_class, true_label).sum().float().item()\n union = torch.logical_or(true_class, true_label).sum().float().item()\n\n iou = (intersect + smooth) / (union +smooth)\n iou_per_class.append(iou)\n return np.nanmean(iou_per_class)","metadata":{"id":"rq3j5vNGwOcA"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**iou**\n\nThis function computes the Intersection over Union (IoU) between a predicted segmentation mask (pred) and a ground truth mask (gt). IoU is a metric that measures the overlap between the predicted and true labels for a specific class. It is calculated as the ratio of the intersection (where both masks agree on the class) to the union (where either mask shows the class).\n\nThe steps followed in the function are:\n\n1. The predicted mask (pred) and ground truth mask (gt) are squeezed to remove any unnecessary dimensions, moved to the CPU (if they're on a GPU), and converted to NumPy arrays.\n2. The function converts the predicted mask into a labeled format using the ToLabel() function, which likely assigns class labels to the prediction.\n3. An aggregated matrix (agg) is created by summing the predicted and ground truth masks. Each pixel value in this matrix represents a combination of the predicted and actual classes.\n4. The number of pixels where both pred and gt are 1 (i.e., both predict the same class) is calculated as i, representing the intersection.\n5. The union (u) is calculated as the number of pixels where either pred or gt (or both) are greater than 0, representing pixels where the class is either predicted or present in the ground truth.\n6. If there are no pixels in the union (u == 0), the IoU is set to 1 (perfect match). Otherwise, the IoU is computed as the ratio i/u (intersection over union).","metadata":{"id":"WYi5BZxUwSL7"}},{"cell_type":"code","source":"def iou(pred, gt):\n pred = pred.squeeze().cpu().data.numpy()\n pred = ToLabel(pred)\n gt = gt.squeeze().cpu().data.numpy()\n agg = pred + gt\n i = float(np.sum(agg == 2))\n u = float(np.sum(agg > 0))\n if u == 0:\n result = 1\n else:\n result = i/u\n return result","metadata":{"id":"6bF9GHOswUM1"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Expliquen la arquitectura que eligieron para su modelo y justifiquen su selección, e.g. U-NET, u otras arquitecturas de segmentación. Describir exactamente la arquitectura, e.g. número de capas, tamaño de kernel de convolución, stride, padding, max pool, etc.**","metadata":{"id":"ApXDBI_2wZIZ"}},{"cell_type":"markdown","source":"In this case, the selected model is a variant of the U-Net architecture using EfficientNet-B2 as the encoder.\n\nEfficientNet-B2 is part of the EfficientNet family, designed to achieve an optimal balance between computational efficiency and accuracy. EfficientNet uses a technique called \"compound scaling\" to simultaneously adjust the depth (number of layers), width (number of filters), and resolution (input size) of the model. It has also proven to be very effective in image classification tasks, achieving high accuracy in visual feature extraction. The architecture includes high-performance convolution blocks and advanced techniques like depthwise separable convolution that optimize precision.\n\nThe encoder in the U-Net with 5 levels of depth allows the model to capture features at different scales. Each additional level in the encoder extracts more complex and abstract features.\n\nThe U-Net architecture includes a decoder that increases the resolution of the extracted features using upsampling operations and convolution layers. This allows for the reconstruction of the segmented image from high-level features.\n\n###Characteristics\n**EfficientNet-B2 (Encoder)** \n- **Number of layers**: EfficientNet-B2 has a total of 16 convolutional layers in its encoder structure.\n- **Convolution kernel size**: EfficientNet-B2 typically uses 3x3 and 5x5 kernels in its main convolutions. \n- **Stride**: The stride varies depending on the layer. In the initial layers of the encoder, the stride is often 2 to quickly reduce the image's dimensionality. In later layers, the stride may be 1 to maintain resolution during deep feature extraction. \n- **Padding**: EfficientNet-B2 uses \"same\" padding in most of its convolutional layers, meaning the output size remains the same as the input after applying the convolution, preserving the resolution throughout the model. \n- **Max Pooling**: EfficientNet does not explicitly use traditional max pooling operations. Instead, it relies on strided convolutions to reduce the feature resolution, improving the model's computational efficiency.\n\n**U-Net (Decoder)** \n- **Number of layers**: The U-Net decoder has 5 levels, with convolutional layers that progressively increase the image resolution. \n- **Convolution kernel size**: In the U-Net decoder, 3x3 kernels are typically used in the convolution operations to accurately capture features during upsampling. \n- **Stride**: The stride in the decoder is usually 1 to maintain the resolution while reconstructing the segmented image. \n- **Padding**: \"Same\" padding is used in the decoder's convolutions to ensure that the output resolution matches the input at each layer. \n- **Max Pooling**: U-Net uses upsampling in the decoder instead of max pooling. The resolution is restored using interpolation or transposed convolutions, rather than pooling.","metadata":{"id":"ttBI5vCZ-Kr1"}},{"cell_type":"code","source":"model = smp.Unet('efficientnet-b2',in_channels=1, encoder_weights='imagenet',classes=4, activation=None, encoder_depth=5, decoder_channels=[256, 128, 64, 32, 16])\n#decoder_attention_type ='scse'","metadata":{"id":"umCs8Ztnwat8"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" **Expliquen cómo entrenaron su modelo, incluyendo:**\n\n La función de pérdida utilizada (criterion).\n El optimizador y learning rate, learning schedule (si aplica), y cualquier otro hyperparameter.\n El número de epochs y el tamaño del minibatch.\n Cualquier técnica de regularización (por ejemplo, weight decay, early stopping) que hayan empleado.","metadata":{"id":"k4TfMmmiwcDa"}},{"cell_type":"markdown","source":"**Fit Function**\nWe define fit function, this function is responsible for training and evaluating the model on a dataset over a specified number of epochs.\n\n\nLoss Function (criterion) The criterion argument represents the loss function used to train the model. Here we use a cross-entropy funcion.\n\nOptimizer, Learning Rate, and Scheduler The optimizer, passed as a parameter (optimizer), can be Adam or SGD, and adjusts the model's weights based on the gradients. The learning rate is controlled through the optimizer and can be dynamically adjusted with a scheduler, as seen in scheduler.step(). The function get_lr(optimizer) retrieves the current learning rate and stores it in the lrs list at each training step. In this case we use Adam.\n\nNumber of Epochs and Minibatch Size The number of epochs (epochs) defines how many times the model will train over the entire dataset. The minibatch size is determined using train_loader and val_loader, which divide the dataset into minibatches. The batch size is defined when creating the DataLoader.\n\nRegularization Techniques Early stopping: If the validation loss does not improve for 7 consecutive epochs (not_improve == 7), training is stopped to avoid overfitting. Weight decay: Although not explicitly mentioned, this could be included in the optimizer's configuration. Best model saving: The model is saved every time the validation loss decreases (if decrease % 5 == 0), preserving the best-performing model during training. Training Loop Steps Training Phase: The model is put in training mode with model.train(). Each batch of images and masks passes through the model, the loss is calculated, backpropagation is performed with loss.backward(), and the optimizer updates the weights using optimizer.step(). The scheduler adjusts the learning rate after each step with scheduler.step(). Validation Phase: After each epoch, the model is evaluated on the validation set using model.eval(), preventing weight updates. The mIoU metrics and pixel accuracy are calculated for both the training and validation sets, along with the validation loss. Logging and Monitoring After each epoch, the training loss, validation loss, accuracy, and mIoU scores are logged to monitor progress. If the validation loss improves, the model is saved.\n\nThis process is designed to train and evaluate a segmentation model, monitoring its performance and adjusting hyperparameters throughout training.","metadata":{"id":"hHBlkHgMkfyg"}},{"cell_type":"code","source":"def get_lr(optimizer):\n for param_group in optimizer.param_groups:\n return param_group['lr']\n\ndef fit(epochs, model, train_loader, val_loader, criterion, optimizer, scheduler, patch=False):\n #torch.cuda.empty_cache()\n train_losses = []\n test_losses = []\n val_iou = []; val_acc = []\n train_iou = []; train_acc = []\n lrs = []\n min_loss = np.inf\n decrease = 1 ; not_improve=0\n\n model.to(device)\n fit_time = time.time()\n for e in range(epochs):\n since = time.time()\n running_loss = 0\n iou_score = 0\n accuracy = 0\n #training loop\n model.train()\n for i, data in enumerate(tqdm(train_loader)):\n #training phase\n image_tiles, mask_tiles = data\n\n image = image_tiles.to(device); mask =mask_tiles.to(device);\n #forward\n output = model(image)\n\n loss = criterion(output, mask)\n #evaluation metrics\n iou_score += mIoU(output, mask)\n accuracy += pixel_accuracy(output, mask)\n #backward\n loss.backward()\n optimizer.step() #update weight\n optimizer.zero_grad() #reset gradient\n\n #step the learning rate\n lrs.append(get_lr(optimizer))\n scheduler.step()\n\n running_loss += loss.item()\n\n\n else:\n model.eval()\n test_loss = 0\n test_accuracy = 0\n val_iou_score = 0\n #validation loop\n with torch.no_grad():\n for i, data in enumerate(tqdm(val_loader)):\n image_tiles, mask_tiles = data\n\n image = image_tiles.to(device); mask =mask_tiles.to(device);\n output = model(image)\n #evaluation metrics\n val_iou_score += mIoU(output, mask)\n test_accuracy += pixel_accuracy(output, mask)\n #loss\n loss = criterion(output, mask)\n test_loss += loss.item()\n\n #calculatio mean for each batch\n train_losses.append(running_loss/len(train_loader))\n test_losses.append(test_loss/len(val_loader))\n\n\n if min_loss > (test_loss/len(val_loader)):\n print('Loss Decreasing.. {:.3f} >> {:.3f} '.format(min_loss, (test_loss/len(val_loader))))\n min_loss = (test_loss/len(val_loader))\n decrease += 1\n if decrease % 5 == 0:\n print('saving model...')\n torch.save(model, 'Unet_efficientnet_b2_mIoU-{:.3f}.pt'.format(val_iou_score/len(val_loader)))\n\n\n if (test_loss/len(val_loader)) > min_loss:\n not_improve += 1\n min_loss = (test_loss/len(val_loader))\n print(f'Loss Not Decrease for {not_improve} time')\n if not_improve == 7:\n print('Loss not decrease for 7 times, Stop Training')\n break\n\n #iou\n val_iou.append(val_iou_score/len(val_loader))\n train_iou.append(iou_score/len(train_loader))\n train_acc.append(accuracy/len(train_loader))\n val_acc.append(test_accuracy/ len(val_loader))\n print(\"Epoch:{}/{}..\".format(e+1, epochs),\n \"Train Loss: {:.3f}..\".format(running_loss/len(train_loader)),\n \"Val Loss: {:.3f}..\".format(test_loss/len(val_loader)),\n \"Train mIoU:{:.3f}..\".format(iou_score/len(train_loader)),\n \"Val mIoU: {:.3f}..\".format(val_iou_score/len(val_loader)),\n \"Train Acc:{:.3f}..\".format(accuracy/len(train_loader)),\n \"Val Acc:{:.3f}..\".format(test_accuracy/len(val_loader)),\n \"Time: {:.2f}m\".format((time.time()-since)/60))\n\n history = {'train_loss' : train_losses, 'val_loss': test_losses,\n 'train_miou' :train_iou, 'val_miou':val_iou,\n 'train_acc' :train_acc, 'val_acc':val_acc,\n 'lrs': lrs}\n print('Total time: {:.2f} m' .format((time.time()- fit_time)/60))\n return history","metadata":{"id":"4bo1eXJ7wgcP"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Use GPU","metadata":{"id":"pfgbUwgO5kEM"}},{"cell_type":"markdown","source":"A GPU device is used for training the set, as it speeds up the process, paralelizing the computations. The training was made in Kaggle, using the GPU P100 accelerator.","metadata":{"id":"kIGzOTm-BMB4"}},{"cell_type":"code","source":"device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint('Using device:', device)\nprint()","metadata":{"id":"SzmWC_at5jHY"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Find max_lr","metadata":{"id":"9mTbRYUV4QKa"}},{"cell_type":"markdown","source":"To find the optimal max learning rate, the model was trained using different values for learning rate, starting from $ 1 \\times 10^{-8} $ to 10. This technique is inspired by a paper published by Leslie N. Smith, called *Cyclical Learning Rates for Training Neural Networks*. The values for the loss at each iteration are stored, using a smoother function to reduce noise. The criterion function, as well as the optimizer are the same ones used for training.\n","metadata":{"id":"Bs82SIk2Ddl6"}},{"cell_type":"code","source":"\ndef find_lr(model, train_loader, criterion, optimizer, init_value=1e-8, final_value=10.0, beta=0.98):\n num_iter = len(train_loader)\n lr_values = []\n losses = []\n avg_loss = 0\n best_loss = float('inf')\n\n model.to(device)\n\n for i, (images, masks) in enumerate(tqdm(train_loader)):\n images = images.to(device)\n masks = masks.to(device)\n\n lr = init_value * (final_value / init_value) ** (i / num_iter)\n optimizer.param_groups[0]['lr'] = lr\n lr_values.append(lr)\n\n optimizer.zero_grad()\n outputs = model(images)\n loss = criterion(outputs, masks)\n avg_loss = beta * avg_loss + (1 - beta) * loss.item()\n smoothed_loss = avg_loss / (1 - beta ** (i + 1))\n losses.append(smoothed_loss)\n\n if smoothed_loss > 4 * best_loss or torch.isnan(loss):\n return lr_values, losses\n\n if smoothed_loss < best_loss:\n best_loss = smoothed_loss\n\n loss.backward()\n optimizer.step()\n\n return lr_values, losses\n\ndef plot_lr_finder(lr_values, losses):\n \"\"\"\n Plots the learning rate finder results.\n\n Args:\n lr_values: A list of learning rates.\n losses: A list of losses.\n \"\"\"\n plt.plot(lr_values, losses)\n plt.xscale('log')\n plt.xlabel('Learning Rate')\n plt.ylabel('Loss')\n plt.title('Learning Rate Finder')\n plt.show()\n\nmodel_lr = smp.Unet('efficientnet-b2',in_channels=1, encoder_weights='imagenet',classes=4, activation=None, encoder_depth=5, decoder_channels=[256, 128, 64, 32, 16])\n#decoder_attention_type ='scse'\nweight_decay = 1e-4\ncriterion = nn.CrossEntropyLoss()\noptimizer = torch.optim.AdamW(model_lr.parameters(), weight_decay=weight_decay)\n# Find the optimal learning rate\nlr_values, losses = find_lr(model_lr, train_dataloader, criterion, optimizer)\n\n# Plot the results\nplot_lr_finder(lr_values, losses)\n","metadata":{"id":"29N59v5b4Toy","outputId":"97b8402f-bdb6-47bb-b9b9-ee690b97e1e0"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To better visualize the data, we can use the rate of change in loss. The point in the graph with the lowest rate of change is the point were the loss decreasses the most, which would be our max learning rate.","metadata":{"id":"69CCoR1OF6vo"}},{"cell_type":"code","source":"# Calculate the derivative of the loss\ndiff = np.diff(losses)\n\n# Adjust lr_values to match the length of diff\nadjusted_lr_values = lr_values[1:]\n\n# Plot the derivative\nplt.plot(adjusted_lr_values, diff)\nplt.xscale('log')\nplt.xlabel('Learning Rate')\nplt.ylabel('Change of Loss')\nplt.title('Rate of change of Loss')\nplt.show()\n\n# Find the index of the minimum derivative\nmin_diff_index = np.argmin(diff)\n\n# Find the corresponding learning rate\nmin_diff_lr = adjusted_lr_values[min_diff_index]\n\n# Print the learning rate corresponding to the lowest derivative of the loss\nprint(\"Learning Rate corresponding to the highest negative change of the loss:\", min_diff_lr)\n","metadata":{"id":"J2Rg_QXF4913","outputId":"713040d8-d3df-4209-ba98-eb960ddec83d"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training","metadata":{"id":"CpRxKN_tDXzj"}},{"cell_type":"markdown","source":"The max learning rate was found in the steps before. The epochs were found by trial and error, 15 being a number were overfitting was minimum and the metrics were higher. The weight decay was not changed from the starting notebook.","metadata":{"id":"-ncZB886wl5y"}},{"cell_type":"code","source":"max_lr = 3e-2\nepoch = 15\nweight_decay = 1e-4\n\ncriterion = nn.CrossEntropyLoss()\noptimizer = torch.optim.AdamW(model.parameters(), lr=max_lr, weight_decay=weight_decay)\nsched = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr, epochs=epoch,\n steps_per_epoch=len(train_dataloader))\n\nhistory = fit(epoch, model, train_dataloader, val_dataloader, criterion, optimizer, sched)","metadata":{"id":"-LANV0FVwkFL"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"torch.save(model, 'Unet-efficientnet.pt')","metadata":{"id":"Pf9thB19wo8M"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After training the model, we can see the plots for loss, mean iou and accuracy","metadata":{"id":"PEQ6VlQxwuDq"}},{"cell_type":"code","source":"def plot_loss(history):\n plt.plot(history['val_loss'], label='val', marker='o')\n plt.plot( history['train_loss'], label='train', marker='o')\n plt.title('Loss per epoch'); plt.ylabel('loss');\n plt.xlabel('epoch')\n plt.legend(), plt.grid()\n plt.show()\n\ndef plot_score(history):\n plt.plot(history['train_miou'], label='train_mIoU', marker='*')\n plt.plot(history['val_miou'], label='val_mIoU', marker='*')\n plt.title('Score per epoch'); plt.ylabel('mean IoU')\n plt.xlabel('epoch')\n plt.legend(), plt.grid()\n plt.show()\n\ndef plot_acc(history):\n plt.plot(history['train_acc'], label='train_accuracy', marker='*')\n plt.plot(history['val_acc'], label='val_accuracy', marker='*')\n plt.title('Accuracy per epoch'); plt.ylabel('Accuracy')\n plt.xlabel('epoch')\n plt.legend(), plt.grid()\n plt.show()","metadata":{"id":"oCDCgY0dw6Kr"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_loss(history)\nplot_score(history)\nplot_acc(history)","metadata":{"id":"wvdvteIzw6_H"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Evaluation metrics","metadata":{"id":"ucYVpGtVvB4S"}},{"cell_type":"markdown","source":"**Loss:** The first graph shows the behavior of the loss over epochs for both training and validation. The loss seems extremely high in the first epoch of validation and then drops almost to zero. This indicates that the model quickly adjusts to the images.","metadata":{"id":"ewDiYe3QvEMn"}},{"cell_type":"markdown","source":"**mIoU (Mean Intersection over Union)**: The second graph shows the mIoU, which measures the overlap between the predicted and ground truth areas for segmentation. This metric is particularly important in image segmentation tasks because it is more robust for evaluating the overlap between the segmented and actual areas. Both training and validation sets have some variation before the 8th epoch, where both start to increase steadily.","metadata":{"id":"2h3xadFBvS8h"}},{"cell_type":"markdown","source":"**Accuracy**: The third graph displays the model’s accuracy, measuring the percentage of correctly segmented pixels. The results show a small gap between training accuracy. Although accuracy is useful, it is not ideal for segmentation tasks where the background class (non-lesion) might dominate the evaluation.","metadata":{"id":"z2pK_8uvvPH1"}},{"cell_type":"markdown","source":"# Results","metadata":{"id":"9ud-coMyWlEe"}},{"cell_type":"markdown","source":"Next, we create a function to predict a mask for an image , the fucntion sets. the model to evaluation mode to stop the training, it makes predictions using the model, it calculates the mIoU score and finally returns the predicted mask and the score.","metadata":{"id":"UAJdttrPd0bX"}},{"cell_type":"code","source":"def predict_image_mask_miou(model, image, mask, mean=[0.485], std=[0.229]):\n model.eval()\n model.to(device); image=image.to(device)\n mask = mask.to(device)\n\n\n with torch.no_grad():\n output = model(image)\n a,b,c,d = output.shape\n score = mIoU(output, mask)\n masked = torch.argmax(output, dim=1)\n masked = masked.cpu().squeeze(0)\n return masked, score, output.permute(0, 2, 3,1)\n\n\n","metadata":{"id":"hVAcEQrwY96M"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then we need another function, that works similalry to the previous one, but instead of calculating the mIoU, it checks how many individual pixels the model predicted correclty, using the metric pixel_accuracy.","metadata":{"id":"Fb8RUiYIftvt"}},{"cell_type":"code","source":"def predict_image_mask_pixel(model, image, mask, mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]):\n model.eval()\n model.to(device); image=image.to(device)\n mask = mask.to(device)\n with torch.no_grad():\n\n output = model(image)\n acc = pixel_accuracy(output, mask)\n masked = torch.argmax(output, dim=1)\n masked = masked.cpu().squeeze(0)\n return masked, acc\n","metadata":{"id":"PCXk8Qq2fo2a"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using another function, mask_to_onehot function converts a colored mask into one-hot encoding.","metadata":{"id":"Gyj2yQQrgIgg"}},{"cell_type":"code","source":"def mask_to_onehot(mask, palette):\n \"\"\"\n Converts a segmentation mask (H, W, C) to (H, W, K) where the last dim is a one\n hot encoding vector, C is usually 1 or 3, and K is the number of class.\n \"\"\"\n semantic_map = []\n for colour in palette:\n equality = np.equal(mask, colour)\n class_map = np.all(equality, axis=-1)\n semantic_map.append(class_map)\n semantic_map = np.stack(semantic_map, axis=-1).astype(np.float32)\n return torch.from_numpy(semantic_map)","metadata":{"id":"8vg-_Zu3frJN"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using the functions defined above, we fetch a batch of images and masks from the validation dataset and pass the image through the model to predict the mask","metadata":{"id":"xldmAvaQhZ9n"}},{"cell_type":"code","source":"image, mask = next(iter(val_dataloader))","metadata":{"id":"sCT30wpNi_Dl"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We then calculate the mIoU score to se how well the predicted mask overlaps with the actual mask, and convert the true mask to one-hot encoding.","metadata":{"id":"c5tmATDzi-4u"}},{"cell_type":"code","source":"pred_mask, score,output = predict_image_mask_miou(model, image, mask)\nsemantic_map = mask_to_onehot(torch.unsqueeze(mask,-1).cpu().numpy(),palette)","metadata":{"id":"sozS-u6ZZBYx"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We finally visualize our results.","metadata":{"id":"5Ryp3r9KjnIQ"}},{"cell_type":"code","source":"visualize(image, semantic_map, pred_batch=output)","metadata":{"id":"EbNNM-mmZDyS"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{"id":"-Oef0BGFiBtL"}},{"cell_type":"markdown","source":"* Yellow: Pixels where both model prediction and true masking are enabled for a specific class. This means that the model has been correct in its prediction for those pixels.\n* Red: The pixels where the model prediction is positive (that is, the model predicts the class), but the true mask is negative. This means that the model has incorrectly identified these pixels as belonging to the class when in fact they are not.\n* Green: Pixels where the true mask is positive (that is, the pixels actually belong to the class), but the model prediction is negative. This means that the model has not correctly identified these pixels that should have been classified as part of the class.\n\nWe can notice that there is a considerable amount of correct values with respect to the Ground Truth, however, we can notice that there are several red pixels, which means that values are being identified in a class when they really are not.","metadata":{"id":"QlUtH0bzoiyH"}},{"cell_type":"markdown","source":"After analyzing, we create a function to calculate the mIoU score for every image in the validation/test set. It loops over the dataset and predicts the mask for each image and calculates their score, it then stores this score in a list which provide an insight on how the model performs accross all images in the test set.","metadata":{"id":"adVGCQe0jtat"}},{"cell_type":"code","source":"def miou_score(model, test_set):\n score_iou = []\n for i, data in enumerate(tqdm(test_set)):\n img, mask = data\n pred_mask, score,output = predict_image_mask_miou(model, img, mask)\n score_iou.append(score)\n return score_iou\n\nmob_miou = miou_score(model, val_dataloader)\nmob_miou\n\n","metadata":{"id":"izVchv-MdMTe"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We clean up the training data","metadata":{"id":"va6xbFzzl7lh"}},{"cell_type":"code","source":"del train_images\ndel train_masks","metadata":{"id":"JbX06MRAmAHO"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then we create a function to normalize and predict the mask for a single image in the test set. using a softmax function to normalize the output probabilities. We then reshape the predicted mask to match the expected shape for processing.\n","metadata":{"id":"rLUtdevkl5kw"}},{"cell_type":"code","source":"def test_predict(model, image, mean=[0.485], std=[0.229]):\n t = T.Compose([T.ToTensor(), T.Normalize(mean, std)])\n image = t(image)\n model.eval()\n\n model.to(device); image=image.to(device)\n\n with torch.no_grad():\n output = model(torch.unsqueeze(image,1))\n output = nn.Softmax(dim=1)(output)\n return output.permute(0, 2, 3,1)","metadata":{"id":"G2E9KfEPdP-h"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We then process a batch of images from the test set and predict masks for each image in the batch, converting them from tensors to arrays. The predicted masks are thresholded to 0.5 and the final masks are displayed along the input images","metadata":{"id":"ivp9tp3HmZvM"}},{"cell_type":"code","source":"image_batch = np.stack([val_augs(image=img)['image'] for img in test_images_medseg], axis=0)\nprint(torch.from_numpy(image_batch).shape)\nprint(image_batch[i].shape)\n#output = test_predict(model, torch.from_numpy(image_batch).permute(0, 3, 1,2))\noutput = np.zeros((10,256,256,4))\nfor i in range(10):\n output[i] = test_predict(model, image_batch[i]).cpu().numpy()\nprint(output.shape)\ntest_masks_prediction = output > 0.5\nvisualize(image_batch, test_masks_prediction, num_samples=len(test_images_medseg))","metadata":{"id":"Z4z-VSE5dSuJ"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{"id":"Il5LxAzniOJL"}},{"cell_type":"markdown","source":"Finally, we resize the predicted masks using scipy,with the goal to enlarge the predicted masks to match the original image size or resolution.","metadata":{"id":"a1bcuiL_7flv"}},{"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","metadata":{"id":"UUoN92Vl7WFO"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then we flatten the masks predicted to a one-dmensional array, we transform the data to a DataFrame and save it as a CSV File, the content of the file will have the index and the predicted mask values for each pixel","metadata":{"id":"8siWmZ-k7w4e"}},{"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('sub.csv')","metadata":{"id":"H5WJL-ub7SOz"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Conclusions*\n\n*Regarding the obtained results, it can be observed that the model shows a good ability to correctly identify the infected regions in lung images, thus demonstrating its precision in detecting significant areas. It achieves correct predictions, where the model is able to segment both small and large areas for a reliable medical diagnosis. However, some areas with incorrectly predicted infections can also be observed. This may be due to the complexity of the images, where normal tissue can appear similar to infected areas under certain conditions. The same issue arises when predicting infected areas, where the model sometimes shows this behavior in image sections that are too small or where the infection is not very visible.\nOverall, the model proved to be effective for lung image segmentation. The identified errors highlight key areas for future improvements, such as adjusting the loss function to better balance the prediction of small infected areas and implementing post-processing techniques to eliminate false positives.*\n*In conclusion, the image results reveal several key observations where false positives can be identified. The model marks areas as infected when they are not. These incorrect predictions, represented in red, appear in various columns and different parts of the images. This suggests that the model might be over-predicting certain regions, indicating that it is detecting patterns that do not actually correspond to COVID-19-affected areas. False positives are a sign that the model may be overfitting the data, meaning it is becoming too specific to the training data and not generalizing well.*\n\n*On the other hand, there are also false negatives, reflected in the green areas of the images, where the model did not correctly detect infected areas that should have been marked. Although false negatives do not seem as significant compared to false positives, their presence suggests that the model is not capturing all affected areas well, especially those that might be small or difficult to detect.*\n\n*Despite these errors, there are many areas where the model made correct predictions, particularly in the last rows of the images, where the yellow areas indicate a good match between the prediction and the ground truth mask. This is a positive sign that the model has a solid understanding of the key features of lung infection, achieving effective segmentation in many of the images.*","metadata":{"id":"XEqP0jupq4aF"}},{"cell_type":"markdown","source":"## References","metadata":{"id":"87urGrZn-KW-"}},{"cell_type":"markdown","source":"Radiopedia source (images): https://radiopaedia.org/articles/covid-19-4?lang=us\n\nMedseg source (images): https://sirm.org/en/\n\nLr tunning : https://arxiv.org/abs/1506.01186\n\nUnet: Ronneberger, O., Fischer, P., & Brox, T. (2015). U-Net: Convolutional Networks for Biomedical Image Segmentation (arXiv:1505.04597). arXiv.\nhttps://doi.org/10.48550/arXiv.1505.04597","metadata":{"id":"h-2lJ5uIAK76"}}]}