{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":10338,"databundleVersionId":862042,"sourceType":"competition"}],"dockerImageVersionId":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Preprocessing\n\n## Introduction\nIn this notebook we will preprocess the data for our classification task.<br />\nWe will train a classifier to predict whether an X-Ray of a patient shows signs of pneumonia or not based on the RSNA Pneumonia Detection Challenge (https://www.kaggle.com/c/rsna-pneumonia-detection-challenge).\n\nDatasets: https://www.kaggle.com/c/rsna-pneumonia-detection-challenge\n\nAt first we download the data from kaggle, by clicking on **\"Download All\"** and extract it afterwards.\n\nAcknowledgements:\nWang X, Peng Y, Lu L, Lu Z, Bagheri M, Summers RM. ChestX-ray8: Hospital-scale Chest X-ray Database and Benchmarks on Weakly-Supervised Classification and Localization of Common Thorax Diseases. IEEE CVPR 2017, http://openaccess.thecvf.com/content_cvpr_2017/papers/Wang_ChestX-ray8_Hospital-Scale_Chest_CVPR_2017_paper.pdf\n\nOriginal Source: https://nihcc.app.box.com/v/ChestXray-NIHCC","metadata":{}},{"cell_type":"markdown","source":"## Data Set and Libraries Imports:\n* pathlib for convenient path handling\n* pydicom for reading dicom files\n* numpy for storing the actual images\n* cv2 (opencv-python) for directly resizing the images\n* pandas to read the provided labels\n* matplotlib for visualizing some images\n* tqdm for nice progress bar\n\n**Concept:** Importing necessary libraries and obtaining a labeled dataset for pneumonia classification.","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\nimport pydicom\nimport numpy as np\nimport cv2\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:55:47.700415Z","iopub.execute_input":"2023-12-10T04:55:47.700909Z","iopub.status.idle":"2023-12-10T04:55:48.803776Z","shell.execute_reply.started":"2023-12-10T04:55:47.700854Z","shell.execute_reply":"2023-12-10T04:55:48.802307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Preprocessing:\n\nAt first, we read the csv file containing the labels\n\nLoad label data from a CSV file, including patient IDs, pneumonia coordinates, and a binary target variable.","metadata":{}},{"cell_type":"code","source":"labels = pd.read_csv(\"/kaggle/input/rsna-pneumonia-detection-challenge/stage_2_train_labels.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:55:48.805924Z","iopub.execute_input":"2023-12-10T04:55:48.806727Z","iopub.status.idle":"2023-12-10T04:55:48.882636Z","shell.execute_reply.started":"2023-12-10T04:55:48.806691Z","shell.execute_reply":"2023-12-10T04:55:48.881684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First 6 entries of the labels dataframe:","metadata":{}},{"cell_type":"code","source":"labels.head(6)","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:55:48.883767Z","iopub.execute_input":"2023-12-10T04:55:48.884166Z","iopub.status.idle":"2023-12-10T04:55:48.919158Z","shell.execute_reply.started":"2023-12-10T04:55:48.884126Z","shell.execute_reply":"2023-12-10T04:55:48.917930Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note that subjects may occur multiple times in the dataset because different pneumonia spots are handled indivually. For our classification task, we can remove those duplicates as we are only interested in the binary label.\n\nOur goal is to train a classifier and not a detector, so we can remove those duplicates.\n\nRemove duplicate entries with the same patient ID.","metadata":{}},{"cell_type":"code","source":"# Remove duplicate entries\nlabels = labels.drop_duplicates(\"patientId\")","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:55:48.922295Z","iopub.execute_input":"2023-12-10T04:55:48.923068Z","iopub.status.idle":"2023-12-10T04:55:48.943322Z","shell.execute_reply.started":"2023-12-10T04:55:48.923021Z","shell.execute_reply":"2023-12-10T04:55:48.941623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labels.head()","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:55:48.946758Z","iopub.execute_input":"2023-12-10T04:55:48.947937Z","iopub.status.idle":"2023-12-10T04:55:48.969033Z","shell.execute_reply.started":"2023-12-10T04:55:48.947847Z","shell.execute_reply":"2023-12-10T04:55:48.967271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Concept:** Organizing and preparing data for further analysis by handling duplicates and setting up file paths.\n\nLet's define the path to the dicom files and also the path were we want to store our processed npy files.\n\nDefine root and save paths for train images and processed files.","metadata":{}},{"cell_type":"code","source":"ROOT_PATH = Path(\"/kaggle/input/rsna-pneumonia-detection-challenge/stage_2_train_images\")\nSAVE_PATH = Path(\"Processed/\")","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:55:48.971382Z","iopub.execute_input":"2023-12-10T04:55:48.971796Z","iopub.status.idle":"2023-12-10T04:55:48.982502Z","shell.execute_reply.started":"2023-12-10T04:55:48.971761Z","shell.execute_reply":"2023-12-10T04:55:48.981097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualization:\n\n- Create a 3x3 grid of X-ray images with corresponding labels.\n- Highlight positive labels (indicating pneumonia) in the central images of the middle and lower rows.\n\n**Concept:** Visualizing a subset of the data to understand the characteristics of X-ray images and their labels.\n\nLet's look at some example images:\n\nTo do so we create 3x3 subplots to visualize 9 x-ray images and their corresponding labels.","metadata":{}},{"cell_type":"code","source":"fig, axis = plt.subplots(3, 3, figsize=(9, 9))\nc = 0\nfor i in range(3):\n    for j in range(3):\n        patient_id = labels.patientId.iloc[c]\n        dcm_path = ROOT_PATH/patient_id\n        dcm_path = dcm_path.with_suffix(\".dcm\")\n        dcm = pydicom.read_file(dcm_path).pixel_array\n        \n        label = labels[\"Target\"].iloc[c]\n        \n        axis[i][j].imshow(dcm, cmap=\"bone\")\n        axis[i][j].set_title(label)\n        c+=1","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:55:48.984654Z","iopub.execute_input":"2023-12-10T04:55:48.986195Z","iopub.status.idle":"2023-12-10T04:55:52.306605Z","shell.execute_reply.started":"2023-12-10T04:55:48.986137Z","shell.execute_reply":"2023-12-10T04:55:52.305359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Image Loading and Preprocessing:\n\n- Loop through each patient ID, standardize and resize images, and store them in the train or validation directory.\n- Divide the pixel array by 255 for standardization.\n- Resize images to 224x224 pixels to reduce computational complexity.\n- Calculate mean and standard deviation during image processing for normalization.\n\n**Concept:** Preparing images for machine learning by standardizing, resizing, and normalizing, which aids in model training and improves convergence.\n\n### Dicom Reading & Effective storage\n\nIn order to efficiently handle our data in the Dataloader, we convert the X-Ray images stored in the DICOM format to numpy arrays. Afterwards we compute the overall mean and standard deviation of the pixels of the whole dataset, for the purpose of normalization.\nThen the created numpy images are stored in two separate folders according to their binary label:\n* 0: All X-Rays which do not show signs of pneumonia\n* 1: All X-Rays which show signs of pneumonia","metadata":{}},{"cell_type":"markdown","source":"To do so, we iterate over the patient ids and concat the patient ID with the ROOT_PATH.\n\nWe then directly save the standardized and resized files into the corresponding directory (0 for healthy, 1 for pneumonia).\nThis allows to take advantage of the ready-to-use torchvision **DatasetFolder** for simple file reading\n\n\nWe standardize all images by the maximum pixel value in the provided dataset, 255.\nAll images are resized to 224x224.\n\nTo compute dataset mean and standard deviation, we compute the sum of the pixel values as well as the sum of the squared pixel values for each subject.\nThis allows to compute the overall mean and standard deviation without keeping the whole dataset in memory.\n","metadata":{}},{"cell_type":"code","source":"sums = 0\nsums_squared = 0\n\nfor c, patient_id in enumerate(tqdm(labels.patientId)):\n    dcm_path = ROOT_PATH/patient_id  # Create the path to the dcm file\n    dcm_path = dcm_path.with_suffix(\".dcm\")  # And add the .dcm suffix\n    \n    # Read the dicom file with pydicom and standardize the array\n    dcm = pydicom.read_file(dcm_path).pixel_array / 255  \n        \n    # Resize the image as 1024x1024 is way to large to be handeled by Deep Learning models at the moment\n    # Let's use a shape of 224x224\n    # In order to use less space when storing the image we convert it to float16\n    dcm_array = cv2.resize(dcm, (224, 224)).astype(np.float16)\n    \n    # Retrieve the corresponding label\n    label = labels.Target.iloc[c]\n    \n    # 4/5 train split, 1/5 val split\n    train_or_val = \"train\" if c < 24000 else \"val\" \n        \n    current_save_path = SAVE_PATH/train_or_val/str(label) # Define save path and create if necessary\n    current_save_path.mkdir(parents=True, exist_ok=True)\n    np.save(current_save_path/patient_id, dcm_array)  # Save the array in the corresponding directory\n    \n    normalizer = dcm_array.shape[0] * dcm_array.shape[1]  # Normalize sum of image\n    if train_or_val == \"train\":  # Only use train data to compute dataset statistics\n        sums += np.sum(dcm_array) / normalizer\n        sums_squared += (np.power(dcm_array, 2).sum()) / normalizer\n","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2023-12-10T04:55:52.308625Z","iopub.execute_input":"2023-12-10T04:55:52.309141Z","iopub.status.idle":"2023-12-10T04:59:33.051732Z","shell.execute_reply.started":"2023-12-10T04:55:52.309098Z","shell.execute_reply":"2023-12-10T04:59:33.047508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Mean and Standard Deviation:\n\n- Convert sum and sum squared variables into mean and standard deviation.\n- These computed values will be used for normalizing the dataset in subsequent steps.\n\n**Concept:** Calculating mean and standard deviation is a common step in normalization, ensuring that input features have similar scales, which helps in training more robust and faster models.","metadata":{}},{"cell_type":"code","source":"mean = sums / 24000\nstd = np.sqrt(sums_squared / 24000 - (mean**2))","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.053472Z","iopub.status.idle":"2023-12-10T04:59:33.054070Z","shell.execute_reply.started":"2023-12-10T04:59:33.053793Z","shell.execute_reply":"2023-12-10T04:59:33.053814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will use mean and std later in the dataloader to normalize our data","metadata":{}},{"cell_type":"code","source":"print(f\"Mean of Dataset: {mean}, STD: {std}\")","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.056272Z","iopub.status.idle":"2023-12-10T04:59:33.057516Z","shell.execute_reply.started":"2023-12-10T04:59:33.057151Z","shell.execute_reply":"2023-12-10T04:59:33.057188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Introduction\nIn this notebook we will create and train the classifier on the chest x-ray dataset to classify whether an image shows signs of pneumonia or not","metadata":{}},{"cell_type":"markdown","source":"## Imports:\n\n* torch and torchvision for model and dataloader creation\n* transforms from torchvision for Data Augmentation and Normalization\n* torchmetrics for easy metric computation\n* pytorch lightning for efficient and easy training implementation\n* ModelCheckpoint and TensorboardLogger for checkpoint saving and logging\n* tqdm for progress par when validating the model\n* numpy for all kinds of stuff :) accessing stored x-ray images\n* matplotlib for visualizing some images\n\n","metadata":{}},{"cell_type":"code","source":"import torch\nimport torchvision\nfrom torchvision import transforms\nimport torchmetrics\nimport pytorch_lightning as pl\nfrom pytorch_lightning.callbacks import ModelCheckpoint\nfrom pytorch_lightning.loggers import TensorBoardLogger\nfrom tqdm.notebook import tqdm\nimport numpy as np\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.059363Z","iopub.status.idle":"2023-12-10T04:59:33.059990Z","shell.execute_reply.started":"2023-12-10T04:59:33.059674Z","shell.execute_reply":"2023-12-10T04:59:33.059703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First we create our dataset.\nWe can leverage the **DatasetFolder** from torchvision: It allows to simply pass a root directory and return return a dataset object with access to all files within the directory and the directory name as class label. <br />\nWe only need to define a loader function, **load_file**, which defines how the files shall be loaded.\nThis is very comfortable as we only have to load our previously stored numpy files.\nAdditionally, we need to define a list of file extensions (just \"npy\" in our case).\n\nFinally we can pass a transformation sequence for Data Augmentation and Normalization.\n\nWe use:\n* RandomResizedCrops which applies a random crop of the image and resizes it to the original image size (224x224)\n* Random Rotations between -5 and 5 degrees\n* Random Translation (max 5%)\n* Random Scaling (0.9-1.1 of original image size)","metadata":{}},{"cell_type":"code","source":"def load_file(path):\n    return np.load(path).astype(np.float32)","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.061770Z","iopub.status.idle":"2023-12-10T04:59:33.062421Z","shell.execute_reply.started":"2023-12-10T04:59:33.062178Z","shell.execute_reply":"2023-12-10T04:59:33.062203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_transforms = transforms.Compose([\n                                    transforms.ToTensor(),  # Convert numpy array to tensor\n                                    transforms.Normalize(0.49, 0.248),  # Use mean and std from preprocessing notebook\n                                    transforms.RandomAffine( # Data Augmentation\n                                        degrees=(-5, 5), translate=(0, 0.05), scale=(0.9, 1.1)),\n                                        transforms.RandomResizedCrop((224, 224), scale=(0.35, 1), antialias=True)\n\n])\n\nval_transforms = transforms.Compose([\n                                    transforms.ToTensor(),  # Convert numpy array to tensor\n                                    transforms.Normalize([0.49], [0.248]),  # Use mean and std from preprocessing notebook\n])\n\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.064221Z","iopub.status.idle":"2023-12-10T04:59:33.064694Z","shell.execute_reply.started":"2023-12-10T04:59:33.064482Z","shell.execute_reply":"2023-12-10T04:59:33.064503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we create the train and val dataset and the corresponding data loaders.\n\nPlease adapt batch size and num_workers according to your hardware ressources.","metadata":{}},{"cell_type":"code","source":"train_dataset = torchvision.datasets.DatasetFolder(\n    \"Processed/train/\",\n    loader=load_file, extensions=\"npy\", transform=train_transforms)\n\nval_dataset = torchvision.datasets.DatasetFolder(\n    \"Processed/val/\",\n    loader=load_file, extensions=\"npy\", transform=val_transforms)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.065649Z","iopub.status.idle":"2023-12-10T04:59:33.066096Z","shell.execute_reply.started":"2023-12-10T04:59:33.065890Z","shell.execute_reply":"2023-12-10T04:59:33.065911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's inspect some augmented train images","metadata":{}},{"cell_type":"code","source":"fig, axis = plt.subplots(2, 2, figsize=(9, 9))\nfor i in range(2):\n    for j in range(2):\n        random_index = np.random.randint(0, 20000)\n        x_ray, label = train_dataset[random_index]\n        axis[i][j].imshow(x_ray[0], cmap=\"bone\")\n        axis[i][j].set_title(f\"Label:{label}\")\n","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.068639Z","iopub.status.idle":"2023-12-10T04:59:33.069119Z","shell.execute_reply.started":"2023-12-10T04:59:33.068909Z","shell.execute_reply":"2023-12-10T04:59:33.068931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batch_size = 64#TODO\nnum_workers = 4# TODO\n\ntrain_loader = torch.utils.data.DataLoader(train_dataset, batch_size=batch_size, num_workers=num_workers, shuffle=True)\nval_loader = torch.utils.data.DataLoader(val_dataset, batch_size=batch_size, num_workers=num_workers, shuffle=False, persistent_workers=True)\n\nprint(f\"There are {len(train_dataset)} train images and {len(val_dataset)} val images\")","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.071073Z","iopub.status.idle":"2023-12-10T04:59:33.071535Z","shell.execute_reply.started":"2023-12-10T04:59:33.071321Z","shell.execute_reply":"2023-12-10T04:59:33.071341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"The classes are imbalanced: There are more images without signs of pneumonia than with pneumonia.\nThere are multiple ways to deal with imbalanced datasets:\n* Weighted Loss\n* Oversampling\n* Doing nothing :)\n\nIn this example, we will simply do nothing as this often yields the best results.\nBuf feel free to play around with a weighted loss. A template to define a customized weighted loss function is provided below.\n\nOversampling will be shown in a later lecture.","metadata":{}},{"cell_type":"code","source":"np.unique(train_dataset.targets, return_counts=True), np.unique(val_dataset.targets, return_counts=True)","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.073149Z","iopub.status.idle":"2023-12-10T04:59:33.073577Z","shell.execute_reply.started":"2023-12-10T04:59:33.073365Z","shell.execute_reply":"2023-12-10T04:59:33.073384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model Creation in pytorch lightning\n\nEach pytorch lightning model is defined by at least an initialization method, a **forward** function which defines the forward pass/prediction, a **training_step** which yields the loss and **configure_optimizers** to specify the optimization algorithm.\n\nAdditionally, we can use a **training_epoch_end** callback to compute overall dataset statistics and metrics such as accuracy.\n\nSubsequently, we define the **validation_step**. The validation step performs more or less the same steps as the training step, however, on the validation data. In this case, pytorch lightning doesn't update the weights.\nAgain, we can use **validation_epoch_end** to compute overall dataset metrics.\n\nNo loops or manual weight updates are needed!<br />\nAdditionally, pl also handles device management.  Just pass the number of GPUS when creating the trainer.","metadata":{}},{"cell_type":"markdown","source":"**Now it is time to create the model** - We will use the ResNet18 network architecture.\n\nAs most of the torchvision models, the original ResNet expects a three channel input in **conv1**. <br />\nHowever, our X-Ray image data has only one channel.\nThus we need to change the in_channel parameter from 3 to 1.\n\nAdditionally, we will change the last fully connected layer to have only one output as we have a binary class label.","metadata":{}},{"cell_type":"markdown","source":"### Optimizer and Loss\nWe use the **Adam** Optimizer with a learning rate of 0.0001 and the **BinaryCrossEntropy** Loss function.<br />\n(In fact we use **BCEWithLogitsLoss** which directly accepts the raw unprocessed predicted values and computes the sigmoid activation function before applying Cross Entropy).\nFeel free to pass a weight different from 1 to the Pneumonia model in order to use the weighted loss function.\n","metadata":{}},{"cell_type":"code","source":"# torchvision.models.resnet18()","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.074853Z","iopub.status.idle":"2023-12-10T04:59:33.075270Z","shell.execute_reply.started":"2023-12-10T04:59:33.075080Z","shell.execute_reply":"2023-12-10T04:59:33.075099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class PneumoniaModel(pl.LightningModule):\n    def __init__(self, weight=1):\n        super().__init__()\n        \n        self.model = torchvision.models.resnet18()\n        # change conv1 from 3 to 1 input channels\n        self.model.conv1 = torch.nn.Conv2d(1, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)\n        # change out_feature of the last fully connected layer (called fc in resnet18) from 1000 to 1\n        self.model.fc = torch.nn.Linear(in_features=512, out_features=1)\n        \n        self.optimizer = torch.optim.Adam(self.model.parameters(), lr=1e-4)\n        self.loss_fn = torch.nn.BCEWithLogitsLoss(pos_weight=torch.tensor([weight]))\n        \n        # simple accuracy computation\n        self.train_acc = torchmetrics.Accuracy(task=\"binary\")\n        self.val_acc = torchmetrics.Accuracy(task=\"binary\")\n\n    def forward(self, data):\n        pred = self.model(data)\n        return pred\n    \n    def training_step(self, batch, batch_idx):\n        x_ray, label = batch\n        label = label.float()  # Convert label to float (just needed for loss computation)\n        pred = self(x_ray)[:,0]  # Prediction: Make sure prediction and label have same shape\n        loss = self.loss_fn(pred, label)  # Compute the loss\n        \n        # Log loss and batch accuracy\n        self.log(\"Train Loss\", loss)\n        self.log(\"Step Train Acc\", self.train_acc(torch.sigmoid(pred), label.int()))\n        return loss\n    \n    \n    def on_train_epoch_end(self):\n        # After one epoch compute the whole train_data accuracy\n        self.log(\"Train Acc\", self.train_acc.compute())\n        \n        \n    def validation_step(self, batch, batch_idx):\n        # Same steps as in the training_step\n        x_ray, label = batch\n        label = label.float()\n        pred = self(x_ray)[:,0]  # make sure prediction and label have same shape\n\n        loss = self.loss_fn(pred, label)\n        \n        # Log validation metrics\n        self.log(\"Val Loss\", loss)\n        self.log(\"Step Val Acc\", self.val_acc(torch.sigmoid(pred), label.int()))\n        return loss\n    \n    def on_validation_epoch_end(self):\n        self.log(\"Val Acc\", self.val_acc.compute())\n    \n    def configure_optimizers(self):\n        #Caution! You always need to return a list here (just pack your optimizer into one :))\n        return [self.optimizer]\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.076723Z","iopub.status.idle":"2023-12-10T04:59:33.077133Z","shell.execute_reply.started":"2023-12-10T04:59:33.076943Z","shell.execute_reply":"2023-12-10T04:59:33.076962Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = PneumoniaModel()  # Instantiate the model","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.078431Z","iopub.status.idle":"2023-12-10T04:59:33.078831Z","shell.execute_reply.started":"2023-12-10T04:59:33.078643Z","shell.execute_reply":"2023-12-10T04:59:33.078662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We create a checkpoint callback which only stores the 10 best models based on the validation accuracy.","metadata":{}},{"cell_type":"code","source":"# Create the checkpoint callback\ncheckpoint_callback = ModelCheckpoint(\n    monitor='Val Acc',\n    save_top_k=10,\n    mode='max')","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.080546Z","iopub.status.idle":"2023-12-10T04:59:33.081022Z","shell.execute_reply.started":"2023-12-10T04:59:33.080791Z","shell.execute_reply":"2023-12-10T04:59:33.080811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Trainer documentation: https://pytorch-lightning.readthedocs.io/en/latest/common/trainer.html","metadata":{}},{"cell_type":"code","source":"# Create the trainer\n# Change the gpus parameter to the number of available gpus on your system. Use 0 for CPU training\n\ngpus = 1 #TODO\n# gpus = 0 #TODO\ntrainer = pl.Trainer(logger=TensorBoardLogger(save_dir=\"./logs\"), log_every_n_steps=1,\n                    callbacks=checkpoint_callback,\n                    max_epochs=35)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.082615Z","iopub.status.idle":"2023-12-10T04:59:33.083036Z","shell.execute_reply.started":"2023-12-10T04:59:33.082825Z","shell.execute_reply":"2023-12-10T04:59:33.082844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainer.fit(model, train_loader, val_loader)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.084412Z","iopub.status.idle":"2023-12-10T04:59:33.084830Z","shell.execute_reply.started":"2023-12-10T04:59:33.084632Z","shell.execute_reply":"2023-12-10T04:59:33.084650Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Evaluation\n\nLet's evaluate our model!\nAt first, we load the latest checkpoint and send the model to the GPU, if possible.","metadata":{}},{"cell_type":"code","source":"device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n\n# Use strict=False, otherwise we would want to match the pos_weight which is not necessary\nmodel = PneumoniaModel.load_from_checkpoint(\"weights/weights_1.ckpt\")\nmodel.eval()\nmodel.to(device);","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compute prediction on the complete validation set and store predictions and labels","metadata":{}},{"cell_type":"code","source":"preds = []\nlabels = []\n\nwith torch.no_grad():\n    for data, label in tqdm(val_dataset):\n        data = data.to(device).float().unsqueeze(0)\n        pred = torch.sigmoid(model(data)[0].cpu())\n        preds.append(pred)\n        labels.append(label)\npreds = torch.tensor(preds)\nlabels = torch.tensor(labels).int()\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compute metrics:\nWe can see that the overall result is already decent with our simple model.<br />\nHowever, we suffer from a large amount of False Negatives due to the data imbalance.<br />\nThis is of particular importance in to avoid in medical imaging as missing findings might be fatal.<br />\nFeel free to try what happens if you increase or decrease the weight in the loss.\n\nAn alternative to retraining with a weighted loss is to reduce the classification threshold from 0.5 to e.g 0.25. It produces way less false negatives but increases the number of False positives. <br />\nThis is called the precision-recall tradeoff","metadata":{}},{"cell_type":"code","source":"acc = torchmetrics.Accuracy()(preds, labels)\nprecision = torchmetrics.Precision()(preds, labels)\nrecall = torchmetrics.Recall()(preds, labels)\ncm = torchmetrics.ConfusionMatrix(num_classes=2)(preds, labels)\ncm_threshed = torchmetrics.ConfusionMatrix(num_classes=2, threshold=0.25)(preds, labels)\n\nprint(f\"Val Accuracy: {acc}\")\nprint(f\"Val Precision: {precision}\")\nprint(f\"Val Recall: {recall}\")\nprint(f\"Confusion Matrix:\\n {cm}\")\nprint(f\"Confusion Matrix 2:\\n {cm_threshed}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axis = plt.subplots(3, 3, figsize=(9, 9))\n\nfor i in range(3):\n    for j in range(3):\n        rnd_idx = np.random.randint(0, len(preds))\n        axis[i][j].imshow(val_dataset[rnd_idx][0][0], cmap=\"bone\")\n        axis[i][j].set_title(f\"Pred:{int(preds[rnd_idx] > 0.5)}, Label:{labels[rnd_idx]}\")\n        axis[i][j].axis(\"off\")\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Congratulation! You made it to the end of the notebook and to the end of the classification lecture.","metadata":{}},{"cell_type":"markdown","source":"# Interpretability\n\n## Visualize the decision of the classifier\nIn this notebook we will learn how to use Class Acvitation Maps (CAM) (https://arxiv.org/abs/1512.04150). <br />\nA method to visualize image regions most important for the decision of the classifier.","metadata":{}},{"cell_type":"markdown","source":"## Imports\n\n* torch for tensor manipulation\n* torchvision for resnet18\n* transforms for Normalization\n* pytorch lightning for model creation\n* numpy for data loading\n* matplotlib for plotting","metadata":{}},{"cell_type":"code","source":"%matplotlib notebook\nimport torch\nimport torchvision\nfrom torchvision import transforms\nimport pytorch_lightning as pl\nimport numpy as np\nimport matplotlib.pyplot as plt","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The dataset of the previous notebook is used.","metadata":{}},{"cell_type":"code","source":"def load_file(path):\n    return np.load(path).astype(np.float32)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-10T04:59:33.086073Z","iopub.status.idle":"2023-12-10T04:59:33.086472Z","shell.execute_reply.started":"2023-12-10T04:59:33.086278Z","shell.execute_reply":"2023-12-10T04:59:33.086297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_transforms = transforms.Compose([\n                                transforms.ToTensor(),\n                                transforms.Normalize(0.49, 0.248),\n\n])\n\nval_dataset = torchvision.datasets.DatasetFolder(\"Processed/val/\", loader=load_file, extensions=\"npy\", transform=val_transforms)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The key idea of CAM is to multiply the output of the last convolutional layer (BasicBlock 1 of layer 4) $A_k$ (consisting of k channels) with the parameters $w$ of the subsequent fully connected layer to compute an activation map $M$:\n$$ M = \\sum_k w_kA_k$$\n\nTo do so, we need to access this particular output of the trained resnet18.<br />\nLet's recap the resnet architecture:","metadata":{}},{"cell_type":"code","source":"temp_model = torchvision.models.resnet18()\ntemp_model","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can convert the network to a generator using the **children()** function.<br />\nThis means that we can use the list function to convert it into a list!<br />\nThe convolutional part of the network comprises all layers up to the AdaptiveAvgPool2d layer.","metadata":{}},{"cell_type":"code","source":"list(temp_model.children())[:-2]  # get all layers up to avgpool","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Using **Sequential** from pytorch, we convert the list of layers back to a Sequential Model.\n","metadata":{}},{"cell_type":"code","source":"torch.nn.Sequential(*list(temp_model.children())[:-2])","metadata":{"execution":{"iopub.status.busy":"2023-12-10T05:00:41.709243Z","iopub.execute_input":"2023-12-10T05:00:41.709727Z","iopub.status.idle":"2023-12-10T05:00:41.749416Z","shell.execute_reply.started":"2023-12-10T05:00:41.709692Z","shell.execute_reply":"2023-12-10T05:00:41.747490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we are ready to go.<br />\nWe add an additional output to the forward function of our pneumonia model, to return the feature maps of the last convolutional layer ($A$)\n\nWe extract the feature map in the forward pass, followed by global average pooling and flattening.\nFinally we use the fully connected layer to compute the final class prediction.","metadata":{}},{"cell_type":"code","source":"class PneumoniaModel(pl.LightningModule):\n    def __init__(self):\n        super().__init__()\n        \n        self.model = torchvision.models.resnet18()\n        # Change conv1 from 3 to 1 input channels\n        self.model.conv1 = torch.nn.Conv2d(1, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)\n        # Change out_feature of the last fully connected layer (called fc in resnet18) from 1000 to 1\n        self.model.fc = torch.nn.Linear(in_features=512, out_features=1)\n        \n        # Extract the feature map\n        self.feature_map = torch.nn.Sequential(*list(self.model.children())[:-2])    \n    def forward(self, data):\n        \n        # Compute feature map\n        feature_map = self.feature_map(data)\n        # Use Adaptive Average Pooling as in the original model\n        avg_pool_output = torch.nn.functional.adaptive_avg_pool2d(input=feature_map, output_size=(1, 1))\n        print(avg_pool_output.shape)\n        # Flatten the output into a 512 element vector\n        avg_pool_output_flattened = torch.flatten(avg_pool_output)\n        print(avg_pool_output_flattened.shape)\n        # Compute prediction\n        pred = self.model.fc(avg_pool_output_flattened)\n        return pred, feature_map\n    \n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def cam(model, img):\n    with torch.no_grad():\n        pred, features = model(img.unsqueeze(0))\n    features = features.reshape((512, 49))\n    weight_params = list(model.model.fc.parameters())[0]\n    weight = weight_params[0].detach()\n    \n    \n    cam = torch.matmul(weight, features)\n    cam_img = cam.reshape(7, 7).cpu()\n    return cam_img, torch.sigmoid(pred)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use strict to prevent pytorch from loading weights for self.feature_map\nmodel = PneumoniaModel.load_from_checkpoint(\"weights/weights_3.ckpt\", strict=False)\nmodel.eval();","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## CAM\nNow we can define the CAM function by using the formula from above:","metadata":{}},{"cell_type":"code","source":"def cam(model, img):\n    \"\"\"\n    Compute class activation map according to cam algorithm\n    \"\"\"\n    with torch.no_grad():\n        pred, features = model(img.unsqueeze(0))\n    b, c, h, w = features.shape\n\n    # We reshape the 512x7x7 feature tensor into a 512x49 tensor in order to simplify the multiplication\n    features = features.reshape((c, h*w))\n    \n    # Get only the weights, not the bias\n    weight_params = list(model.model.fc.parameters())[0] \n    \n    # Remove gradient information from weight parameters to enable numpy conversion\n    weight = weight_params[0].detach()\n    print(weight.shape)\n    # Compute multiplication between weight and features with the formula from above.\n    # We use matmul because it directly multiplies each filter with the weights\n    # and then computes the sum. This yields a vector of 49 (7x7 elements)\n    cam = torch.matmul(weight, features)\n    print(features.shape)\n    \n    ### The following loop performs the same operations in a less optimized way\n    #cam = torch.zeros((7 * 7))\n    #for i in range(len(cam)):\n    #    cam[i] = torch.sum(weight*features[:,i])\n    ##################################################################\n    \n    # Normalize and standardize the class activation map (Not always necessary, thus not shown in the lecture)\n    cam = cam - torch.min(cam)\n    cam_img = cam / torch.max(cam)\n    # Reshape the class activation map to 512x7x7 and move the tensor back to CPU\n    cam_img = cam_img.reshape(h, w).cpu()\n\n    return cam_img, torch.sigmoid(pred)\n\ndef visualize(img, heatmap, pred):\n    \"\"\"\n    Visualization function for class activation maps\n    \"\"\"\n    img = img[0]\n    # Resize the activation map of size 7x7 to the original image size (224x224)\n    heatmap = transforms.functional.resize(heatmap.unsqueeze(0), (img.shape[0], img.shape[1]))[0]\n    \n    # Create a figure\n    fig, axis = plt.subplots(1, 2)\n    \n    axis[0].imshow(img, cmap=\"bone\")\n    # Overlay the original image with the upscaled class activation map\n    axis[1].imshow(img, cmap=\"bone\")\n    axis[1].imshow(heatmap, alpha=0.5, cmap=\"jet\")\n    plt.title(f\"Pneumonia: {(pred > 0.5).item()}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize(img, cam, pred):\n    img = img[0]\n    cam = transforms.functional.resize(cam.unsqueeze(0), (224, 224))[0]\n    \n    fig, axis = plt.subplots(1, 2)\n    axis[0].imshow(img, cmap=\"bone\")\n    axis[1].imshow(img, cmap=\"bone\")\n    axis[1].imshow(cam, alpha=0.5, cmap=\"jet\")\n    plt.title(pred)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Demo Time","metadata":{}},{"cell_type":"code","source":"img = val_dataset[-6][0]  # Select a subject\nactivation_map, pred = cam(model, img)  # Compute the Class activation map given the subject","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize(img, activation_map, pred)  # Visualize CAM","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The heatmap focuses on the area which shows signs of pneumonia, isn't that cool?","metadata":{}},{"cell_type":"markdown","source":"Awesome, you made it! <br />\nYou can try out this method for all neural networks with one restriction: The classic CAM algorithm requrires a specific network architecture:<br />\nLast Convolution $\\rightarrow$ Global Average Pool $\\rightarrow$ 1 FC Layer.\n\nFor other network architectures, you either need to adapt the architecture, followed by retraining our use different methods such as GradCAM or ScoreCAM (https://arxiv.org/abs/1610.02391), (https://arxiv.org/abs/1910.01279)","metadata":{}}]}