{"cells":[{"metadata":{},"cell_type":"markdown","source":"# AIML Online Capstone -Pneumonia Detection Challenge\n\n## What is Pneumonia?\n\n**Pneumonia** is an infection in one or both lungs. Bacteria, viruses, and fungi cause it. The infection causes inflammation in the air sacs in your lungs, which are called alveoli.\n\nThe alveoli fill with fluid or pus, making it difficult to breathe.Pneumonia is a lung infection that can range from mild to so severe that you have to go to the hospital.\n![](https://www.mayoclinic.org/-/media/kcms/gbs/patient-consumer/images/2016/05/18/13/02/ww5r032t-8col-jpg.jpg)\n\nPneumonia accounts for over 15% of all deaths of children under 5 years old internationally. In 2017, 920,000 children under the age of 5 died from the disease. It requires review of a chest radiograph (CXR) by highly trained specialists and confirmation through clinical history, vital signs and laboratory exams. Pneumonia usually manifests as an area or areas of increased opacity on CXR. However, the diagnosis of pneumonia on CXR is complicated because of a number of other conditions in the lungs such as fluid overload (pulmonary edema), bleeding, volume loss (atelectasis or collapse), lung cancer, or post-radiation or surgical changes. Outside of the lungs, fluid in the pleural space (pleural effusion) also appears as increased opacity on CXR. When available, comparison of CXRs of the patient taken at different time points and correlation with clinical symptoms and history are helpful in making the diagnosis.\n\nCXRs are the most commonly performed diagnostic imaging study. A number of factors such as positioning of the patient and depth of inspiration can alter the appearance of the CXR, complicating interpretation further. In addition, clinicians are faced with reading high volumes of images every shift.\n\n## Pneumonia Detection\n\nNow to detect Pneumonia, we need to detect **Inflammation** of the lungs. In this project, you’re challenged to build an algorithm to detect a visual signal for pneumonia in medical images. Specifically, your algorithm needs to automatically locate lung opacities on chest radiographs.\n\n## How Is Pneumonia Diagnosed?\nSometimes pneumonia can be difficult to diagnose because the symptoms are so variable, and are often very similar to those seen in a cold or influenza. To diagnose pneumonia, and to try to identify the germ that is causing the illness, your doctor will ask questions about your medical history, do a physical exam, and run some tests.\n\n### Medical history\nYour doctor will ask you questions about your signs and symptoms, and how and when they began. To help figure out if your infection is caused by bacteria, viruses or fungi, you may be asked some questions about possible exposures, such as:\n\n### Physical exam\nYour doctor will listen to your lungs with a stethoscope. If you have pneumonia, your lungs may make crackling, bubbling, and rumbling sounds when you inhale.\n\n### Diagnostic Tests\nIf your doctor suspects you may have pneumonia, they will probably recommend some tests to confirm the diagnosis and learn more about your infection. These may include:\n\n1 Blood tests to confirm the infection and to try to identify the germ that is causing your illness.\n\n2) Chest X-ray to look for the location and extent of inflammation in your lungs.\n\n3) Pulse oximetry to measure the oxygen level in your blood. Pneumonia can prevent your lungs from moving enough oxygen into your bloodstream.\n\n4) Sputum test on a sample of mucus (sputum) taken after a deep cough, to look for the source of the infection. If you are considered a high-risk patient because of your age and overall health, or if you are hospitalized, the doctors may want to do some additional tests, including:\n\n5) CT scan of the chest to get a better view of the lungs and look for abscesses or other complications.\n\n6) Arterial blood gas test, to measure the amount of oxygen in a blood sample taken from an artery, usually in      your wrist. This is more accurate than the simpler pulse oximetry.\n\n7) Pleural fluid culture, which removes a small amount of fluid from around tissues that surround the lung, to      analyze and identify bacteria causing the pneumonia.\n\n8) Bronchoscopy, a procedure used to look into the lungs' airways. If you are hospitalized and your treatment is    not working well, doctors may want to see whether something else is affecting your airways, such as a            blockage. They may also take fluid samples or a biopsy of lung tissue.\n\n\n## Business Domain Value\n\nAutomating Pneumonia screening in chest radiographs, providing affected area details through bounding box. \n\nAssist physicians to make better clinical decisions or even replace human judgement in certain functional areas of healthcare (eg, radiology).\n\nGuided by relevant clinical questions, powerful AI techniques can unlock clinically relevant information hidden in the massive amount of data, which in turn can assist clinical decision making.\n\n## Image DataType\n\nMedical images are stored in a special format called DICOM files (`*.dcm`). They contain a combination of header metadata as well as underlying raw image arrays for pixel data.\n\nDataset link: https://www.kaggle.com/c/rsna-pneumonia-detection-challenge/data\n\n\n## Prediction Output\n\nIn this project, we have to predict whether pneumonia exists in a given image. This is done by predicting bounding boxes around areas of the lung. Samples without bounding boxes are negative and contain no definitive evidence of pneumonia. Samples with bounding boxes indicate evidence of pneumonia.\n\nWhen making predictions, the model should predict as many bounding boxes as necessary, in the format: `confidence x-min y-min width height`\n\nThere will be only ONE predicted row per image. This row may include multiple bounding boxes.\n\nA properly formatted row may look like any of the following.\n\nFor patientIds with no predicted pneumonia / bounding boxes: `0004cfab-14fd-4e49-80ba-63a80b6bddd6,`\n\nFor patientIds with a single predicted bounding box: `0004cfab-14fd-4e49-80ba-63a80b6bddd6,0.5 0 0 100 100`\n\nFor patientIds with multiple predicted bounding boxes: `0004cfab-14fd-4e49-80ba-63a80b6bddd6,0.5 0 0 100 100 0.5 0 0 100 100,` etc.\n\nThe general format is as follows:\n\n`patientId,{confidence x-min y-min width height},{confidence x-min y-min width height}, etc.`\n\n## Dataset File Description\n\n`stage_2_train_labels.csv` - the training set. Contains `patientId`s and bounding box / target information.\n\n`stage_2_detailed_class_info.csv` - provides detailed information about the type of positive or negative class for each image.\n\n## Data Fields\n\n- `patientId_` - Each patientId corresponds to a unique image.\n- `x_` - the upper-left x coordinate of the bounding box.\n- `y_` - the upper-left y coordinate of the bounding box.\n- `width_` - the width of the bounding box.\n- `height_` - the height of the bounding box.\n- `Target_` - the binary Target, indicating whether this sample has evidence of pneumonia.\n\n\n## Lung Opacity\n\nTissues with sparse material, such as lungs which are full of air, do not absorb the X-rays and appear black in the image. Dense tissues such as bones absorb X-rays and appear white in the image.\n\nWhile we are theoretically detecting “lung opacities”, there are lung opacities that are not pneumonia related.\n\nIn the data, some of these are labeled “Not Normal No Lung Opacity”.This extra third class indicates that while pneumonia was determined not to be present, there was nonetheless some type of abnormality on the image and often times this finding may mimic the appearance of true pneumonia.\n\nIt's important to note that the various shades of gray in the chest X-Ray refer to the following:\n\n- **Black** = Air\n- **White** = Bone\n- **Grey** = Tissue or Fluid\n\nIn a normal image (shown above) we see the lungs as black, but they have different projections on them - mainly the rib cage bones, main airways, blood vessels and the heart.\n\nIn case of pneumonia, a **haziness** (also referred to as **consolidation**) is present in the chest x-ray image. \n\nImages with no lung opacity and no pneumonia are images where the patient can have rounded hazy boundaries or masses (probably because of lung nodules or masses which can be because of cancer).\n\nThere are other exceptional cases as well where there can be no lung opacity but no pneumonia either. Some of these cases include **pneumonectomy** (lung removed by surgery), **enlarged heart**, **pleural effusion**, etc. \n\n**Reference**: https://www.kaggle.com/zahaviguy/what-are-lung-opacities","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Data Visualization\n\n**Reference**: https://www.kaggle.com/peterchang77/exploratory-data-analysis\n\nLet's start the project with visualizing the data we have. Since the images are in DICOM format, we will use the Python module - `pydicom` for reading these images. We will also use `pylab` module for Data Visualization purposes, `pandas` for EDA, `numpy` for numerical processing and `glob` for directory listing based operations.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import glob, pylab \nimport pandas as pd\nimport pydicom\nimport numpy as np\nimport numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import Rectangle\nimport seaborn as sns\nimport gc\nimport os\nimport cv2\nimport warnings\nwarnings.simplefilter(action = 'ignore')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0f368672adb58cc7bb57ecc95a04487251bb9c2c"},"cell_type":"markdown","source":"## Directory Structure\n\nThe challenge data is organized in several files and folders.","execution_count":null},{"metadata":{"trusted":true,"_uuid":"1b6c13edf09c7d80145173edad847c34db1ba5f2"},"cell_type":"code","source":"!ls ../input","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"88fb25d6864f224a0995a28140de532f80cb3e6b"},"cell_type":"markdown","source":"The several key items in this folder:\n* `stage_2_train_labels.csv`: CSV file containing training set patientIds and  labels (including bounding boxes)\n* `stage_2_detailed_class_info.csv`: CSV file containing detailed labels (explored further below)\n* `stage_2_train_images/`:  directory containing training set raw image (DICOM) files\n* `stage_2_test_images/` : directory containing testing set image (DICOM) files\n\nLet's go ahead and take a look at the first labels CSV file first:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"b4c54974f191b477e46370c79254f75657ed2a85","scrolled":true},"cell_type":"code","source":"df = pd.read_csv('../input/stage_2_train_labels.csv')\nprint(df.iloc[0])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5435a871306c834a2282ae867a86861f9b80a610"},"cell_type":"markdown","source":"As you can see, each row in the CSV file contains a `patientId` (one unique value per patient), a target (either 0 or 1 for absence or presence of pneumonia, respectively) and the corresponding abnormality bounding box defined by the upper-left hand corner (x, y) coordinate and its corresponding width and height. In this particular case, the patient does *not* have pneumonia and so the corresponding bounding box information is set to `NaN`. See an example case with pnuemonia here:","execution_count":null},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"98c60a2ddf5bba070c0908d9bb705c0c5976ac7a"},"cell_type":"code","source":"print(df.iloc[4])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7ad152504ec74be4ac1992ca669da7ad3ce3aac2"},"cell_type":"markdown","source":"One important thing to keep in mind is that a given `patientId` may have **multiple** boxes if more than one area of pneumonia is detected.","execution_count":null},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"# Overview of DICOM files and medical images\n\nMedical images are stored in a special format known as DICOM files (`*.dcm`). They contain a combination of header metadata as well as underlying raw image arrays for pixel data. In Python, one popular library to access and manipulate DICOM files is the `pydicom` module. To use the `pydicom` library, first find the DICOM file for a given `patientId` by simply looking for the matching file in the `stage_2_train_images/` folder, and the use the `pydicom.read_file()` method to load the data:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"5f2c15162a0d1390624b42ef94d4f9e260be56ac"},"cell_type":"code","source":"patientId = df['patientId'][0]\ndcm_file = '../input/stage_2_train_images/%s.dcm' % patientId\ndcm_data = pydicom.read_file(dcm_file)\n\nprint(dcm_data)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2cfbe9eb43f4e4922c42481739046943767f765c"},"cell_type":"markdown","source":"Most of the standard headers containing patient identifable information have been anonymized (removed) so we are left with a relatively sparse set of metadata. The primary field we will be accessing is the underlying pixel data as follows:","execution_count":null},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"953ccb9e176398bcad8dad0241dbc1e194fa8a8e"},"cell_type":"code","source":"im = dcm_data.pixel_array\nprint(type(im))\nprint(im.dtype)\nprint(im.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ed51c8db0a041023f7fb19ecd66b08bb630a1485"},"cell_type":"markdown","source":"## Considerations\n\nAs we can see here, the pixel array data is stored as a Numpy array, a powerful numeric Python library for handling and manipulating matrix data (among other things). In addition, it is apparent here that the original radiographs have been preprocessed for us as follows:\n\n* The relatively high dynamic range, high bit-depth original images have been rescaled to 8-bit encoding (256 grayscales). For the radiologists out there, this means that the images have been windowed and leveled already. In clinical practice, manipulating the image bit-depth is typically done manually by a radiologist to highlight certain disease processes. To visually assess the quality of the automated bit-depth downscaling and for considerations on potentially improving this baseline, consider consultation with a radiologist physician.\n\n* The relativley large original image matrices (typically acquired at >2000 x 2000) have been resized to the data-science friendly shape of 1024 x 1024. For the purposes of this challenge, the diagnosis of most pneumonia cases can typically be made at this resolution. To visually assess the feasibility of diagnosis at this resolution, and to determine the optimal resolution for pneumonia detection (oftentimes can be done at a resolution *even smaller* than 1024 x 1024), consider consultation with a radiogist physician.\n\n## Visualizing An Example\n\nTo take a look at this first DICOM image, let's use the `pylab.imshow()` method:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"59ac176a83e8fc1f11531a216723791577948fbe"},"cell_type":"code","source":"pylab.imshow(im, cmap=pylab.cm.gist_gray)\npylab.axis('off')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"78c0e24742ce4f61c153f363e66ac5d7e5efbb0e"},"cell_type":"markdown","source":"# Data Overview\n\nAs alluded to above, any given patient may potentially have many boxes if there are several different suspicious areas of pneumonia. To collapse the current CSV file dataframe into a dictionary with unique entries, consider the following method:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"626ee0aa4038a988bac8e8c800c4bf4270af1b76"},"cell_type":"code","source":"def parse_data(df):\n    \"\"\"\n    Method to read a CSV file (Pandas dataframe) and parse the \n    data into the following nested dictionary:\n\n      parsed = {\n        \n        'patientId-00': {\n            'dicom': path/to/dicom/file,\n            'label': either 0 or 1 for normal or pnuemonia, \n            'boxes': list of box(es)\n        },\n        'patientId-01': {\n            'dicom': path/to/dicom/file,\n            'label': either 0 or 1 for normal or pnuemonia, \n            'boxes': list of box(es)\n        }, ...\n\n      }\n\n    \"\"\"\n    # --- Define lambda to extract coords in list [y, x, height, width]\n    extract_box = lambda row: [row['y'], row['x'], row['height'], row['width']]\n\n    parsed = {}\n    for n, row in df.iterrows():\n        # --- Initialize patient entry into parsed \n        pid = row['patientId']\n        if pid not in parsed:\n            parsed[pid] = {\n                'dicom': '../input/stage_2_train_images/%s.dcm' % pid,\n                'label': row['Target'],\n                'boxes': []}\n\n        # --- Add box if opacity is present\n        if parsed[pid]['label'] == 1:\n            parsed[pid]['boxes'].append(extract_box(row))\n\n    return parsed","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cbb1a9cf9eb93df25896f47f80257aba561c0d00"},"cell_type":"markdown","source":"Let's use the method here:","execution_count":null},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"fbf2c2f60dcc661c249989eb063382ad7304aa10"},"cell_type":"code","source":"parsed = parse_data(df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7b851e847ce442d690a8dc55937070235c7a16f0"},"cell_type":"markdown","source":"As we saw above, patient `00436515-870c-4b36-a041-de91049b9ab4` has pnuemonia so lets check our new `parsed` dict here to see the patients corresponding bounding boxes:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"459c07f8fadd3f30aff6ff26bdd4b4c8fe4f79f5"},"cell_type":"code","source":"print(parsed['00436515-870c-4b36-a041-de91049b9ab4'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d54c891b2477a2c022b8b4e3ffd73422eaf862e1"},"cell_type":"markdown","source":"# Visualizing Boxes\n\nIn order to overlay color boxes on the original grayscale DICOM files, consider using the following  methods (below, the main method `draw()` requires the method `overlay_box()`):","execution_count":null},{"metadata":{"trusted":true,"_uuid":"6be7b744e42ac7204ea661c0c7175881b043c5c6"},"cell_type":"code","source":"def draw(data):\n    \"\"\"\n    Method to draw single patient with bounding box(es) if present \n\n    \"\"\"\n    # --- Open DICOM file\n    d = pydicom.read_file(data['dicom'])\n    im = d.pixel_array\n\n    # --- Convert from single-channel grayscale to 3-channel RGB\n    im = np.stack([im] * 3, axis=2)\n\n    # --- Add boxes with random color if present\n    for box in data['boxes']:\n        rgb = np.floor(np.random.rand(3) * 256).astype('int')\n        im = overlay_box(im=im, box=box, rgb=rgb, stroke=6)\n\n    pylab.imshow(im, cmap=pylab.cm.gist_gray)\n    pylab.axis('off')\n\ndef overlay_box(im, box, rgb, stroke=1):\n    \"\"\"\n    Method to overlay single box on image\n\n    \"\"\"\n    # --- Convert coordinates to integers\n    box = [int(b) for b in box]\n    \n    # --- Extract coordinates\n    y1, x1, height, width = box\n    y2 = y1 + height\n    x2 = x1 + width\n\n    im[y1:y1 + stroke, x1:x2] = rgb\n    im[y2:y2 + stroke, x1:x2] = rgb\n    im[y1:y2, x1:x1 + stroke] = rgb\n    im[y1:y2, x2:x2 + stroke] = rgb\n\n    return im","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d10aef5768f6147aaa4c5c8dec7444fb90f0842d"},"cell_type":"markdown","source":"As we saw above, patient `00436515-870c-4b36-a041-de91049b9ab4` has pnuemonia so let's take a look at the overlaid bounding boxes:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"ee8f69746904ae88da12f631b6c9e0b1e214c4fa","scrolled":true},"cell_type":"code","source":"draw(parsed['00436515-870c-4b36-a041-de91049b9ab4'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"eb00fb542999b5ffe6b32ee38321830f9b090d26"},"cell_type":"markdown","source":"## Exploring Detailed Labels\n\nIn this challenge, the primary endpoint will be the detection of bounding boxes consisting of a binary classification---e.g. the presence or absence of pneumonia. However, in addition to the binary classification, each bounding box *without* pneumonia is further categorized into *normal* or *no lung opacity / not normal*. This extra third class indicates that while pneumonia was determined not to be present, there was nonetheless some type of abnormality on the image---and oftentimes this finding may mimic the appearance of true pneumonia. Keep in mind that this extra class is provided as supplemental information to help improve algorithm accuracy if needed; generation of this separate class **will not** be a formal metric used to evaluate performance in this competition.\n\nAs above, we saw that the first patient in the CSV file did not have pneumonia. Let's look at the detailed label information for this patient:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"a639baad366af644bdd8efe63712dd47cc0fe9ed"},"cell_type":"code","source":"df_detailed = pd.read_csv('../input/stage_2_detailed_class_info.csv')\nprint(df_detailed.iloc[0])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8bb1f27ca89d8f05142e7ac462ccbde61d4228bd"},"cell_type":"markdown","source":"As we see here, the patient does not have pneumonia however *does* have another imaging abnormality present. Let's take a closer look:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"ce3843dbf0f47cda1299d05fa8791148c1c2e806"},"cell_type":"code","source":"patientId = df_detailed['patientId'][0]\ndraw(parsed[patientId])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b54a1136f966df74c15872a073fc333a8dd37294"},"cell_type":"markdown","source":"While the image displayed inline within the notebook is small, it is evident that the patient has several well circumscribed nodular densities in the left lung (right side of image). This can be because of lung cancer masses.","execution_count":null},{"metadata":{"_uuid":"f52f87aafb11d94916d984e19089c28e25b95eed"},"cell_type":"markdown","source":"## Label Summary\n\nFinally, let us take a closer look at the distribution of labels in the dataset. To do so we will first parse the detailed label information:","execution_count":null},{"metadata":{"trusted":true,"_uuid":"44cdd1bff2692927dabd6d08f231be9aae4d31c9"},"cell_type":"code","source":"summary = {}\nfor n, row in df_detailed.iterrows():\n    if row['class'] not in summary:\n        summary[row['class']] = 0\n    summary[row['class']] += 1\n    \nprint(summary)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7e1575f979fe05c0853ba19bf28b56099fd38d7a"},"cell_type":"markdown","source":"As we can see, there is a relatively even split between the three classes, with nearly 2/3rd of the data comprising of no pneumonia (either completely *normal* or *no lung opacity / not normal*). Compared to most medical imaging datasets, where the prevalence of disease is quite low, this dataset has been significantly enriched with pathology.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Before we move on with further EDA on the Metadata, let's have a look at some more images from the three categories to get a better understanding of the dataset.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def show_dicom_image(data_df):\n        img_data = list(data_df.T.to_dict().values())\n        f, ax = plt.subplots(2,2, figsize=(16,18))\n        for i,data_row in enumerate(img_data):\n            pid = data_row['patientId']\n            dcm_file = '../input/stage_2_train_images/%s.dcm' % pid\n            dcm_data = pydicom.read_file(dcm_file)                    \n            ax[i//2, i%2].imshow(dcm_data.pixel_array, cmap=plt.cm.bone)\n            ax[i//2, i%2].set_title('ID: {}\\n Age: {} Sex: {}'.format(\n                data_row['patientId'],dcm_data.PatientAge, dcm_data.PatientSex))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will start off with the images of patients with pnuemonia.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_orig = df.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.concat([df_orig,df_detailed[\"class\"]],axis=1,sort=False)\nshow_dicom_image(df[df['Target']==1].sample(n=4))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, let's have a look at chest x-rays of patients who do not have pneumonia but don't have a normal chest x-ray either.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"show_dicom_image(df[ (df['Target']==0) & (df['class']=='No Lung Opacity / Not Normal')].sample(n=4))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, let's visualize a few images from `\"Normal\"` class as well.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"show_dicom_image(df[ (df['Target']==0) & (df['class']=='Normal')].sample(n=4))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can also display the positive pnuemonia chest x-rays with bounding boxes.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def show_dicome_with_boundingbox(data_df):\n    img_data = list(data_df.T.to_dict().values())\n    f, ax = plt.subplots(2,2, figsize=(16,18))\n    for i,data_row in enumerate(img_data):\n        pid = data_row['patientId']\n        dcm_file = '../input/stage_2_train_images/%s.dcm' % pid\n        dcm_data = pydicom.read_file(dcm_file)                    \n        ax[i//2, i%2].imshow(dcm_data.pixel_array, cmap=plt.cm.bone)\n        ax[i//2, i%2].set_title('ID: {}\\n Age: {} Sex: {}'.format(\n                data_row['patientId'],dcm_data.PatientAge, dcm_data.PatientSex))\n        rows = data_df[data_df['patientId']==data_row['patientId']]\n        box_data = list(rows.T.to_dict().values())        \n        for j, row in enumerate(box_data):            \n            x,y,width,height = row['x'], row['y'],row['width'],row['height']\n            rectangle = Rectangle(xy=(x,y),width=width, height=height, color=\"red\",alpha = 0.1)\n            ax[i//2, i%2].add_patch(rectangle)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_dicome_with_boundingbox(df[df['Target']==1].sample(n=4))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we have a basic idea regarding the labels and the images, we can explore the metadata present in the dataset to see whether some information can be extracted directly from them. This will be extremely important since images will require deep learning based solutions which will take large amount of resources and time for coming up with efficient results, whereas analysing and extracting useful information from metadata is extremely fast and easy to do.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# EDA on Metadata\n\nSince the `target` information is directly given in the metadata, we can approach the problem using 2 approaches:\n\n1. Predicting `target` based on metadata and finding the bounding boxes for the images with `target=1` using the deep neural network model.\n2. Predicting `target` and bounding boxes based on DNN model.\n\nLet's try the first approach first. For this, we will start off with the EDA on metadata.\n\n**Reference**: https://www.kaggle.com/aantonova/practical-eda-on-numerical-data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# For data visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import gc","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load detailed class information\ndetailed_class_info = pd.read_csv('../input/stage_2_detailed_class_info.csv')\n# Load training dataset labels\ntrain_labels = pd.read_csv('../input/stage_2_train_labels.csv')\n\n# Merge the above data information into one dataframe\ndf = pd.merge(left = detailed_class_info, right = train_labels, how = 'left', on = 'patientId')\n\n# Remove the original dataframes since they don't hold any useful information now\ndel detailed_class_info, train_labels\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Display information about the merged dataframe\ndf.info(null_counts = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, there are 37.6k patients' information. Out of these patients, only for 16.9k patients the bounding box coordinates are available meaning that there was no Pneumonia detected for the rest of the patients. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# First 5 rows of the dataframe\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's first start with data cleaning. The first step is to remove the duplicates, if any.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df.drop_duplicates(inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.info()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are now left with data for 30.2k unique patients. Out of these only 9.5k patients have pneumonia.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Next, we know that it's possible to have multiple bounding boxes in an image corresponding to a patient. Let's see this.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df['patientId'].value_counts().head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's a see sample of this data for the first patient id - `32408669-c137-4e8d-bd62-fe8345b40e73`","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['patientId'] == '32408669-c137-4e8d-bd62-fe8345b40e73']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, this patient's chest x-ray image has 4 bounding boxes.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Let's see the distribution of number of bounding boxes in patients' x-ray images.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df['patientId'].value_counts().value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We know that a patient without Pneumonia will have only one row (0 bounding boxes). Let's see how many such cases are there.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['Target'] == 0]['patientId'].value_counts().value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see here, most of the patients have no pneumonia. But this also means that there are some other cases as well. Let's explore the same.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['Target'] == 1]['patientId'].value_counts().value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Notice how that even for patients with pneumonia, there are 2614 patients which have only one row.\n\nNext, let's have a look at the distribution of the target class.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Distribution of `class`","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(x = 'class', hue = 'Target', data = df)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Notice here that most of the images are from the class `No Lung Opacity/Not Normal` or `Normal class`.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Let's have a look at the `target` distribution for these classes. Again, \n\n- `Lung Opacity` - 1\n- `Normal Class` - 0\n- `No Lung Opacity/Not Normal` - 0","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['class'] == 'Lung Opacity']['Target'].value_counts(dropna = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['class'] == 'No Lung Opacity / Not Normal']['Target'].value_counts(dropna = False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['class'] == 'Normal']['Target'].value_counts(dropna = False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature Engineering\n\nNext, using the metadata we have, we will create some new features regarding the bounding box.\n\n1. `x_2`, `y_2` - Corner point coordinates opposite to `x`, `y`\n2. `area` - Area of the bounding box\n3. `x_center`, `y_center` - Center point of the bounding box","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_areas = df.dropna()[['x', 'y', 'width', 'height']].copy()\ndf_areas['x_2'] = df_areas['x'] + df_areas['width']\ndf_areas['y_2'] = df_areas['y'] + df_areas['height']\ndf_areas['x_center'] = df_areas['x'] + df_areas['width'] / 2\ndf_areas['y_center'] = df_areas['y'] + df_areas['height'] / 2\ndf_areas['area'] = df_areas['width'] * df_areas['height']\n\ndf_areas.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we have some new features, let's see if there is any relationship between these features using a `jointplot`.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def createJointplot(df, x, y):\n    sns.jointplot(x = x, y = y, data = df, kind = 'hex', gridsize = 20)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"createJointplot(df_areas,'x','y')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"createJointplot(df_areas,'x_center','y_center')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"createJointplot(df_areas,'x_2','y_2')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, there is no significant correlation between the points (x,y) and (x_2, y_2). Let's also study the correlation between width and height.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"createJointplot(df_areas,'width','height')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see from the above joint plot, width and height have high correlation, so while building a model, we can choose to keep only one of the 2 for increasing the simplicity of the model.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Outlier Analysis\n\nIn this section, we will use boxplots to see whether there are any outliers in the features.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"n_columns = 3\nn_rows = 3\n_, axes = plt.subplots(n_rows, n_columns, figsize=(8 * n_columns, 5 * n_rows))\nfor i, c in enumerate(df_areas.columns):\n    sns.boxplot(y = c, data = df_areas, ax = axes[i // n_columns, i % n_columns])\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are some outliers in the `width` and `height` features. Let's find the rows where these values exist.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_areas[df_areas['width'] > 500]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pid_width = list(df[df['width'] > 500]['patientId'].values)\ndf[df['patientId'].isin(pid_width)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_areas[df_areas['height'] > 900].shape[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pid_height = list(df[df['height'] > 900]['patientId'].values)\ndf[df['patientId'].isin(pid_height)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's drop all these rows for the 2 patients. This will remove the extreme outliers.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df[~df['patientId'].isin(pid_width + pid_height)]\ndf.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are now left with around 30k patients.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Metadata Cleanup\n\nFinally, let's clean up the entire metadata to keep only the relevant columns. We will use the cleaned up data to see whether there is any correlation between the `target` column and other columns.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_meta = df.drop('class', axis = 1).copy()\ndcm_columns = None\n\nfor n, pid in enumerate(df_meta['patientId'].unique()):\n    if n%1000==0:\n        print(n,len(df_meta['patientId'].unique()))\n    dcm_file = '../input/stage_2_train_images/%s.dcm' % pid\n    dcm_data = pydicom.read_file(dcm_file)\n    \n    if not dcm_columns:\n        dcm_columns = dcm_data.dir()\n        dcm_columns.remove('PixelSpacing')\n        dcm_columns.remove('PixelData')\n    \n    for col in dcm_columns:\n        if not (col in df_meta.columns):\n            df_meta[col] = np.nan\n        index = df_meta[df_meta['patientId'] == pid].index\n        df_meta.loc[index, col] = dcm_data.data_element(col).value\n        \n    del dcm_data\n    \ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_meta.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"to_drop = df_meta.nunique()\nto_drop = to_drop[(to_drop <= 1) | (to_drop == to_drop['patientId'])].index\nto_drop = to_drop.drop('patientId')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_meta.drop(to_drop, axis = 1, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_meta.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# try:\n#     df_meta.drop('ReferringPhysicianName', axis = 1, inplace = True)\n# except:\n#     print(\"Referring Physician Name not found\")\ndf_meta['PatientAge'] = df_meta['PatientAge'].astype(int)\ndf_meta['SeriesDescription'] = df_meta['SeriesDescription'].map({'view: AP': 'AP', 'view: PA': 'PA'})\n\ndf_meta.drop('SeriesDescription', axis = 1, inplace = True)\n\ndf_meta['PatientSex'] = df_meta['PatientSex'].map({'F': 0, 'M': 1})\ndf_meta['ViewPosition'] = df_meta['ViewPosition'].map({'PA': 0, 'AP': 1})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_meta.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now using this cleaned data, let's first see if there is any relation of `target` variable and the categorical features like patient's age, sex and the view position of the chest x-ray.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (25, 5))\nsns.countplot(x = 'PatientAge', hue = 'Target', data = df_meta)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see from the above count plot, patients with age lying in the mid range have higher number of patients. The ratio of patients with target 0 versus target 1 is mostly same throughout the range. There does not seem to be any significant correlation between target and patient age.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(x = 'PatientSex', hue = 'Target', data = df_meta)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(x = 'ViewPosition', hue = 'Target', data = df_meta)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's also find out the correlation between the features.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_meta.corr()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see from the above correlation matrix, there is a high correlation between `ViewPosition` and `Target` features.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"At this point, we can try using the cleaned metadata for predicting the `target` column.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from lightgbm import LGBMClassifier\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.model_selection import KFold\n\nimport warnings\nwarnings.simplefilter(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def fast_lgbm_cv_scores(df, target, task, rs = 0):\n    \n    clf = LGBMClassifier(n_estimators = 10000, nthread = 4, random_state = rs)\n    metric = 'auc'\n\n    # Cross validation model\n    folds = KFold(n_splits = 2, shuffle = True, random_state = rs)\n        \n    # Create arrays and dataframes to store results\n    pred = np.zeros(df.shape[0])\n    \n    feats = df.columns.drop(target)\n    \n    feature_importance_df = pd.DataFrame(index = feats)\n    \n    for n_fold, (train_idx, valid_idx) in enumerate(folds.split(df[feats], df[target])):\n        train_x, train_y = df[feats].iloc[train_idx], df[target].iloc[train_idx]\n        valid_x, valid_y = df[feats].iloc[valid_idx], df[target].iloc[valid_idx]\n\n        clf.fit(train_x, train_y, \n                eval_set = [(valid_x, valid_y)], eval_metric = metric, \n                verbose = -1, early_stopping_rounds = 100)\n\n        pred[valid_idx] = clf.predict_proba(valid_x, num_iteration = clf.best_iteration_)[:, 1]\n        \n        feature_importance_df[n_fold] = pd.Series(clf.feature_importances_, index = feats)\n        \n        del train_x, train_y, valid_x, valid_y\n        gc.collect()\n\n    return feature_importance_df, pred, roc_auc_score(df[target], pred)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Using the 4 columns we have from the metadata - `patientId`, `x`, `y`, `width`, `height`, we can fit an LGBM Classifier and find out the area under the ROC curve. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"f_imp, _, score = fast_lgbm_cv_scores(df_meta.drop(['patientId', 'x', 'y', 'width', 'height'], axis = 1), \n                                      target = 'Target', task = 'classification')\nprint('ROC-AUC for Target = {}'.format(score))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f_imp","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see from the area under the ROC curve is quite high (~0.75) which means that there is significant information we can get from the metadata regarding the `target` column. This shows that the first approach that we proposed has some significant potential.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Model Building\n\nIn the previous section, we saw that metadata does play an important role when it comes to predicting the `target` column but including the metadata and the images in the pipeline will mean having three networks:\n1. One network for getting information out of the metadata\n2. One network for processing the chest x-ray image\n3. One network for combining the outputs of the above two networks.\n\nBefore we start with the above, let's try a more simplistic approach. In this approach, we will focus only on the chest x-ray images and train a simple CNN model to understand the performance we are getting out of it. We can then modify the layers and the activation functions to understand the effect they have on the performance. But, going by the typical computer vision solutions, it's highly unlikely that we will get excellent results from the shallow network. That's why, we will next shift to transfer learning using pre-trained models like VGG-16, VGG-19 and then to even more advanced networks like Mask-RCNN and YOLOv3. Due to the recent release of YOLOv4 which has been shown to give better results, we can also train a YOLOv4 model and compare the results. But, in the interim report, we will focus primarily on training a shallow CNN model and using VGG models for transfer learning.\n\nBefore we jump into model building, let's first generate the data which we can use later on in model training.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport csv\nimport random\nimport pydicom\nimport numpy as np\nimport pandas as pd\nfrom skimage import measure\nfrom skimage.transform import resize\nimport matplotlib.patches as patches\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\nfrom tensorflow import keras\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pneumonia_locations = {}\n# load table\nwith open(os.path.join('../input/stage_2_train_labels.csv'),\n          'r') as infile:\n    # open reader\n    reader = csv.reader(infile)\n    # skip header\n    next(reader, None)\n    # loop through rows\n    for rows in reader:\n        # retrieve information\n        filename = rows[0]\n        location = rows[1:5]\n        pneumonia = rows[5]\n        # if row contains pneumonia add label to dictionary\n        # which contains a list of pneumonia locations per filename\n        if pneumonia == '1':\n            # convert string to float to int\n            location = [int(float(i)) for i in location]\n            # save pneumonia location in dictionary\n            if filename in pneumonia_locations:\n                pneumonia_locations[filename].append(location)\n            else:\n                pneumonia_locations[filename] = [location]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load and shuffle filenames\nfolder = '../input/stage_2_train_images'\nfilenames = os.listdir(folder)\nrandom.shuffle(filenames)\n# split into train and validation filenames\nn_valid_samples = 8000\ntrain_filenames = filenames[n_valid_samples:]\nvalid_filenames = filenames[:n_valid_samples]\nprint('n train samples', len(train_filenames))\nprint('n valid samples', len(valid_filenames))\nn_train_samples = len(filenames) - n_valid_samples","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"### Data generator\n\nThe dataset is too large to fit into memory, so we need to create a generator that loads data on the fly.\n\nThe generator takes in some filenames, batch_size and other parameters.\n\nThe generator outputs a random batch of numpy images and numpy masks.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class generator(keras.utils.Sequence):    \n    def __init__(self, folder, filenames, pneumonia_locations=None, batch_size=32, image_size=256, shuffle=True, augment=False, predict=False):\n        self.folder = folder\n        self.filenames = filenames\n        self.pneumonia_locations = pneumonia_locations\n        self.batch_size = batch_size\n        self.image_size = image_size\n        self.shuffle = shuffle\n        self.augment = augment\n        self.predict = predict\n        self.on_epoch_end()\n        \n    def __load__(self, filename):\n        # load dicom file as numpy array\n        img = pydicom.dcmread(os.path.join(self.folder, filename)).pixel_array\n        # create empty mask\n        msk = np.zeros(img.shape)\n        # get filename without extension\n        filename = filename.split('.')[0]\n        # if image contains pneumonia\n        if filename in self.pneumonia_locations:\n            # loop through pneumonia\n            for location in self.pneumonia_locations[filename]:\n                # add 1's at the location of the pneumonia\n                x, y, w, h = location\n                msk[y:y+h, x:x+w] = 1\n        # resize both image and mask\n        img = resize(img, (self.image_size, self.image_size), mode='reflect')\n        msk = resize(msk, (self.image_size, self.image_size), mode='reflect') > 0.5\n        # if augment then horizontal flip half the time\n        if self.augment and random.random() > 0.5:\n            img = np.fliplr(img)\n            msk = np.fliplr(msk)\n        # add trailing channel dimension\n        img = np.expand_dims(img, -1)\n        msk = np.expand_dims(msk, -1)\n        return img, msk\n    \n    def __loadpredict__(self, filename):\n        # load dicom file as numpy array\n        img = pydicom.dcmread(os.path.join(self.folder, filename)).pixel_array\n        # resize image\n        img = resize(img, (self.image_size, self.image_size), mode='reflect')\n        # add trailing channel dimension\n        img = np.expand_dims(img, -1)\n        return img\n        \n    def __getitem__(self, index):\n        # select batch\n        filenames = self.filenames[index*self.batch_size:(index+1)*self.batch_size]\n        # predict mode: return images and filenames\n        if self.predict:\n            # load files\n            imgs = [self.__loadpredict__(filename) for filename in filenames]\n            # create numpy batch\n            imgs = np.array(imgs)\n            return imgs, filenames\n        # train mode: return images and masks\n        else:\n            # load files\n            items = [self.__load__(filename) for filename in filenames]\n            # unzip images and masks\n            imgs, msks = zip(*items)\n            # create numpy batch\n            imgs = np.array(imgs)\n            msks = np.array(msks)\n            return imgs, msks\n        \n    def on_epoch_end(self):\n        if self.shuffle:\n            random.shuffle(self.filenames)\n        \n    def __len__(self):\n        if self.predict:\n            # return everything\n            return int(np.ceil(len(self.filenames) / self.batch_size))\n        else:\n            # return full batches only\n            return int(len(self.filenames) / self.batch_size)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's also define the evaluation metrics that we are going to use.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# define iou or jaccard loss function\ndef iou_loss(y_true, y_pred):\n    #print(y_true)\n    y_true=tf.cast(y_true, tf.float32)\n    y_pred=tf.cast(y_pred, tf.float32)\n    y_true = tf.reshape(y_true, [-1])\n    y_pred = tf.reshape(y_pred, [-1])\n   \n    intersection = tf.reduce_sum(y_true * y_pred)\n    score = (intersection + 1.) / (tf.reduce_sum(y_true) + tf.reduce_sum(y_pred) - intersection + 1.)\n    return 1 - score\n\n# combine bce loss and iou loss\ndef iou_bce_loss(y_true, y_pred):\n    return 0.5 * keras.losses.binary_crossentropy(y_true, y_pred) + 0.5 * iou_loss(y_true, y_pred)\n\n# mean iou as a metric\ndef mean_iou(y_true, y_pred):\n    y_pred = tf.round(y_pred)\n    intersect = tf.reduce_sum(y_true * y_pred, axis=[1, 2, 3])\n    union = tf.reduce_sum(y_true, axis=[1, 2, 3]) + tf.reduce_sum(y_pred, axis=[1, 2, 3])\n    smooth = tf.ones(tf.shape(intersect))\n    return tf.reduce_mean((intersect + smooth) / (union - intersect + smooth))\n\ndef create_downsample(channels, inputs):\n    x = keras.layers.BatchNormalization(momentum=0.9)(inputs)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(channels, 1, padding='same', use_bias=False)(x)\n    x = keras.layers.MaxPool2D(2)(x)\n    return x\n\ndef create_resblock(channels, inputs):\n    x = keras.layers.BatchNormalization(momentum=0.9)(inputs)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(x)\n    x = keras.layers.BatchNormalization(momentum=0.9)(x)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(x)\n    return keras.layers.add([x, inputs])\n\ndef create_network(input_size, channels, n_blocks=2, depth=4):\n    # input\n    inputs = keras.Input(shape=(input_size, input_size, 1))\n    x = keras.layers.Conv2D(channels, 3, padding='same', use_bias=False)(inputs)\n    # residual blocks\n    for d in range(depth):\n        channels = channels * 2\n        x = create_downsample(channels, x)\n        for b in range(n_blocks):\n            x = create_resblock(channels, x)\n    # output\n    x = keras.layers.BatchNormalization(momentum=0.9)(x)\n    x = keras.layers.LeakyReLU(0)(x)\n    x = keras.layers.Conv2D(1, 1, activation='sigmoid')(x)\n    outputs = keras.layers.UpSampling2D(2**depth)(x)\n    model = keras.Model(inputs=inputs, outputs=outputs)\n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we will define the batch size.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 128\nIMAGE_SIZE = 128","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's now start training the model.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model = create_network(input_size=IMAGE_SIZE, channels=32, n_blocks=2, depth=4)\nmodel.compile(optimizer='adam', loss=iou_bce_loss, metrics=['accuracy', mean_iou])\n\n# cosine learning rate annealing\ndef cosine_annealing(x):\n    lr = 0.0001\n    epochs = 3\n    return lr*(np.cos(np.pi*x/epochs)+1.)/2\n\n\nlearning_rate = tf.keras.callbacks.LearningRateScheduler(cosine_annealing)\n\n# create train and validation generators\nfolder = '../input/stage_2_train_images'\ntrain_gen = generator(folder, train_filenames, pneumonia_locations, batch_size=BATCH_SIZE, \n                      image_size=IMAGE_SIZE, shuffle=True, augment=False, predict=False)\nvalid_gen = generator(folder, valid_filenames, pneumonia_locations, batch_size=BATCH_SIZE, \n                      image_size=IMAGE_SIZE, shuffle=False, predict=False)\n\nprint(model.summary())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"EPOCHS=5\nMULTI_PROCESSING = True \n\nhistory = model.fit_generator(train_gen, validation_data=valid_gen, callbacks=[learning_rate], epochs=EPOCHS, \n                              workers=4, use_multiprocessing=MULTI_PROCESSING)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's now plot the training & validation accuracy, loss and IoU values.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(12,4))\nplt.subplot(131)\nplt.plot(history.epoch, history.history[\"loss\"], label=\"Train loss\")\nplt.plot(history.epoch, history.history[\"val_loss\"], label=\"Valid loss\")\nplt.legend()\nplt.subplot(132)\nplt.plot(history.epoch, history.history[\"accuracy\"], label=\"Train accuracy\")\nplt.plot(history.epoch, history.history[\"val_accuracy\"], label=\"Valid accuracy\")\nplt.legend()\nplt.subplot(133)\nplt.plot(history.epoch, history.history[\"mean_iou\"], label=\"Train iou\")\nplt.plot(history.epoch, history.history[\"val_mean_iou\"], label=\"Valid iou\")\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, let's use the model we have trained to predict the output using the validation generator.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"i=0\nfor imgs, msks in valid_gen:    \n    # predict batch of images\n    preds = model.predict(imgs)\n    # create figure\n    f, axarr = plt.subplots(4, 8, figsize=(20,15))\n    axarr = axarr.ravel()\n    axidx = 0\n    # loop through batch\n    for img, msk, pred in zip(imgs, msks, preds):\n        i=i+1\n        #exit after 32 images\n        if i>32:\n            break\n        # plot image\n        axarr[axidx].imshow(img[:, :, 0])\n        # threshold true mask\n        comp = msk[:, :, 0] > 0.5\n        # apply connected components\n        comp = measure.label(comp)\n        # apply bounding boxes\n        predictionString = ''\n        for region in measure.regionprops(comp):\n            # retrieve x, y, height and width\n            y, x, y2, x2 = region.bbox\n            height = y2 - y\n            width = x2 - x\n            axarr[axidx].add_patch(patches.Rectangle((x,y),width,height,linewidth=2,\n                                                     edgecolor='b',facecolor='none'))\n        # threshold predicted mask\n        comp = pred[:, :, 0] > 0.5\n        # apply connected components\n        comp = measure.label(comp)\n        # apply bounding boxes\n        predictionString = ''\n        for region in measure.regionprops(comp):\n            # retrieve x, y, height and width\n            y, x, y2, x2 = region.bbox\n            height = y2 - y\n            width = x2 - x\n            axarr[axidx].add_patch(patches.Rectangle((x,y),width,height,linewidth=2,\n                                                     edgecolor='r',facecolor='none'))\n        axidx += 1\n    plt.show()\n    # only plot one batch\n    break","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## VGG 16 and VGG 19\n\nNext, let's use a pre-trained model instead of training our own CNN model from scratch.\n\nFor training the VGG-16 and VGG-19 models, we will have to change the directory structure slightly.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport csv\nimport random\nimport pydicom\nimport numpy as np\nimport pandas as pd\nfrom skimage import measure\nfrom skimage.transform import resize\nimport matplotlib.patches as patches\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\nfrom tensorflow import keras\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm import tqdm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df=pd.read_csv('../input/stage_2_train_labels.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['path']='../input/stage_2_train_images/'+df['patientId'].astype(str)+'.dcm'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"negative=df[df['Target']==0]\nprint(len(negative))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"positive=df[df['Target']==1]\nunique_positive=positive[['path','patientId']]\npath=unique_positive['path'].unique()\npatientId=unique_positive['patientId'].unique()\n\nunique_positive=pd.DataFrame({'path':path,'patientId':patientId})\nlen(unique_positive)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"try:\n    os.mkdir('/kaggle/working/data')\n    os.mkdir('/kaggle/working/data/positive')\n    os.mkdir('/kaggle/working/data/negative')\n    os.chdir('/kaggle/working')\nexcept:\n    pass","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for _,row in tqdm(unique_positive.iterrows()):\n    img=pydicom.read_file(row['path']).pixel_array\n    img=resize(img,(256,256))\n    plt.imsave('data/positive/'+row['patientId']+'.jpg',img,cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for _,row in tqdm(negative.iterrows()):\n    img=pydicom.read_file(row['path']).pixel_array\n    img=resize(img,(256,256))\n    plt.imsave('data/negative/'+row['patientId']+'.jpg',img,cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.applications.vgg19 import VGG19,preprocess_input\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"datagen=ImageDataGenerator(samplewise_center=True,samplewise_std_normalization=True,horizontal_flip=True,\n                          width_shift_range=0.05,rescale=1/255,fill_mode='nearest',height_shift_range=0.05,\n                           preprocessing_function=preprocess_input,validation_split=0.3,\n                          )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create data-generators for training and validation/testing\ntrain=datagen.flow_from_directory('data',color_mode='rgb',batch_size=32,\n                                  class_mode='binary',subset='training')\ntest=datagen.flow_from_directory('data',color_mode='rgb',batch_size=32,\n                                 class_mode='binary',subset='validation')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will remove the last layer and add our own layers to the model.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pre_trained_model = VGG19(input_shape = (256,256,3), \n                                include_top = False, \n                                weights = 'imagenet')\n\nfor layer in pre_trained_model.layers:\n    layer.trainable = False\n\n# pre_trained_model.summary()\n\nlast_layer = pre_trained_model.get_layer('block5_pool')\nprint('last layer output shape: ', last_layer.output_shape)\nlast_output = last_layer.output","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.layers import Flatten,Dense,Dropout,BatchNormalization,LeakyReLU,ReLU,GaussianDropout","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = Flatten()(last_output)\nmodel = Dense(1024)(model)\nmodel=LeakyReLU(0.1)(model)\nmodel=Dropout(0.25)(model)\nmodel=BatchNormalization()(model)\nmodel = Dense(1024)(model)\nmodel=LeakyReLU(0.1)(model)\nmodel=Dropout(0.25)(model)\nmodel=BatchNormalization()(model)\nmodel = Dense(1, activation='sigmoid')(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.models import Model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fmodel = Model( pre_trained_model.input, model) \n\nfmodel.compile(optimizer = 'adam', \n              loss = 'binary_crossentropy', \n              metrics = ['accuracy'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.callbacks import EarlyStopping,ReduceLROnPlateau","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"early=EarlyStopping(monitor='accuracy',patience=3,mode='auto')\nreduce_lr = ReduceLROnPlateau(monitor='accuracy', factor=0.5, \n                              patience=2, verbose=1,cooldown=0, mode='auto',min_delta=0.0001, min_lr=1e-5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<!--\n**Note** - The training was done on Kaggle kernel for VGG-19. Because of 30 hours restriction, only 20 steps per epoch were chosen, which resulted in a drop in accuracy. If time permits, we will re-run the training with 100 steps per epoch.\n-->","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# class_weight={0:1,1:3.3}\n# # Train model\n# fmodel.fit(train,epochs=30,callbacks=[reduce_lr],\n#            steps_per_epoch=100,validation_data=test,class_weight=class_weight)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# fmodel.save('/kaggle/working/model_vgg19.h5')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # Plot accuracy\n# plt.figure(figsize=(30,20))\n# val_acc=np.asarray(fmodel.history.history['val_accuracy'])*100\n# acc=np.asarray(fmodel.history.history['accuracy'])*100\n# acc=pd.DataFrame({'val_acc':val_acc,'acc':acc})\n# acc.plot(figsize=(20,10),yticks=range(50,100,5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # Plot loss\n# loss=fmodel.history.history['loss']\n# val_loss=fmodel.history.history['val_loss']\n# loss=pd.DataFrame({'val_loss':val_loss,'loss':loss})\n# loss.plot(figsize=(20,10))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the above graph, we can see the onset of overfitting around 30 epochs. But, from the accuracy graphs, we can see that since we stopped the training at 40 epochs, the model has not overfitted to a large extent.\n\n**Note** - Training was stopped at 40 epochs since the 30 hours quota in Kaggle was completed.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# y=[]\n\n# test.reset()\n\n# for i in tqdm(range(84)):\n#     _,tar=test.__getitem__(i)\n#     for j in tar:\n#         y.append(j)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test.reset()\n# y_pred=fmodel.predict(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# pred=[]\n# for i in y_pred:\n#     if i[0]>=0.5:\n#         pred.append(1)\n#     else:\n#         pred.append(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# from sklearn.metrics import roc_curve,auc,precision_recall_curve,classification_report","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # Classification report\n# print(classification_report(y,pred[:len(y)]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the F1-score for normal category is very high - 0.86, whereas for pneumonia, it's just 0.56. This clearly shows the effect of imbalanced dataset. Let's have a look at the area under the ROC curve to assess the model performance.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# plt.figure(figsize=(30,20))\n# fpr,tpr,_=roc_curve(y,y_pred[:len(y)])\n# area_under_curve=auc(fpr,tpr)\n# print('The area under the curve is:',area_under_curve)\n# # Plot area under curve\n# plt.plot(fpr,tpr,'b.-')\n# plt.xlabel('false positive rate')\n# plt.ylabel('true positive rate')\n# plt.plot(fpr,fpr,linestyle='--',color='black')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The area under ROC curve is 0.86 which is much higher than 0.50 (random model), this means that the model has actually learnt and is performing much better than a random model.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Next, we will replace the LeakyReLU layer with ReLU layer to see the effect it has on performance metrics.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model = Flatten()(last_output)\nmodel = Dense(1024)(model)\nmodel=ReLU(0.1)(model)\nmodel=Dropout(0.25)(model)\nmodel=BatchNormalization()(model)\nmodel = Dense(1024)(model)\nmodel=ReLU(0.1)(model)\nmodel=Dropout(0.25)(model)\nmodel=BatchNormalization()(model)\nmodel = Dense(1, activation='sigmoid')(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.models import Model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fmodel = Model( pre_trained_model.input, model) \n\nfmodel.compile(optimizer = 'adam', \n              loss = 'binary_crossentropy', \n              metrics = ['accuracy'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.callbacks import EarlyStopping,ReduceLROnPlateau","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"early=EarlyStopping(monitor='accuracy',patience=3,mode='auto')\nreduce_lr = ReduceLROnPlateau(monitor='accuracy', factor=0.5, \n                              patience=2, verbose=1,cooldown=0, mode='auto',min_delta=0.0001, min_lr=1e-5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<!--\n**Note** - The training was done on Kaggle kernel for VGG-19. Because of 30 hours restriction, only 20 steps per epoch were chosen, which resulted in a drop in accuracy. If time permits, we will re-run the training with 100 steps per epoch.\n-->","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class_weight={0:1,1:3.3}\n# Train model\nfmodel.fit(train,epochs=30,callbacks=[reduce_lr],\n           steps_per_epoch=100,validation_data=test,class_weight=class_weight)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot accuracy\nplt.figure(figsize=(30,20))\nval_acc=np.asarray(fmodel.history.history['val_accuracy'])*100\nacc=np.asarray(fmodel.history.history['accuracy'])*100\nacc=pd.DataFrame({'val_acc':val_acc,'acc':acc})\nacc.plot(figsize=(20,10),yticks=range(50,100,5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot loss\nloss=fmodel.history.history['loss']\nval_loss=fmodel.history.history['val_loss']\nloss=pd.DataFrame({'val_loss':val_loss,'loss':loss})\nloss.plot(figsize=(20,10))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y=[]\n\ntest.reset()\n\nfor i in tqdm(range(84)):\n    _,tar=test.__getitem__(i)\n    for j in tar:\n        y.append(j)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.reset()\ny_pred=fmodel.predict(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred=[]\nfor i in y_pred:\n    if i[0]>=0.5:\n        pred.append(1)\n    else:\n        pred.append(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import roc_curve,auc,precision_recall_curve,classification_report","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Classification report\nprint(classification_report(y,pred[:len(y)]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the F1-score for normal category is very high - 0.86, whereas for pneumonia, it's just 0.56. This clearly shows the effect of imbalanced dataset. Let's have a look at the area under the ROC curve to assess the model performance.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(30,20))\nfpr,tpr,_=roc_curve(y,y_pred[:len(y)])\narea_under_curve=auc(fpr,tpr)\nprint('The area under the curve is:',area_under_curve)\n# Plot area under curve\nplt.plot(fpr,tpr,'b.-')\nplt.xlabel('false positive rate')\nplt.ylabel('true positive rate')\nplt.plot(fpr,fpr,linestyle='--',color='black')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The area under ROC curve is 0.86 which is much higher than 0.50 (random model), this means that the model has actually learnt and is performing much better than a random model.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Future Work\n\nSo far, we have trained a CNN model from scratch, and used transfer learning for training with pre-trained models like VGG-16 and VGG-19. We saw that the CNN model gave a decent IoU (both training and validation) of around 0.7, whereas the VGG models gave good accuracies for the classification problem. In next steps, we will try to further utilize the power of transfer learning by training models like YOLO, Mask-RCNN, etc. which are much better suited for object detection problems, where we also need to give the bounding box coordinates in the prediction, just as is expected in the current problem statement.","execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}