{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<img src=\"https://storage.googleapis.com/kaggle-competitions/kaggle/26680/logos/header.png?t=2021-04-23-22-04-05)\" />","metadata":{}},{"cell_type":"markdown","source":"# Project Description","metadata":{}},{"cell_type":"markdown","source":"## Understanding the challenge\n\nThis challenge, as well as the dataset itself, is composed of two levels. The first is the image level which contains the chest radiographs, and above it we have the study level, which contains the general conclusion from all the patient radiographs.  \nOn the study level, each study is classified by specialists as `Negative for Pneumonia`, or as `Typical Appearance`, `Indeterminate  Appearance`, or `Atypical Appearance` to Covid-19.\nThe grading system is based on [this paper](https://journals.lww.com/thoracicimaging/Fulltext/2020/11000/Review_of_Chest_Radiograph_Findings_of_COVID_19.4.aspx) which proposes a new reporting language for chest radiographs (CXR) findings related to COVID-19, as described in the following table (Table 1 in the paper):\n\n\n> |Radiographic Classification | CXR Findings | Suggested Reporting Language|\n|:--------------------------|:-------------|:----------------------------|\n|Typical appearance|Multifocal bilateral, peripheral opacities Opacities with rounded morphology Lower lung–predominant distribution|“Findings typical of COVID-19 pneumonia are present. However, these can overlap with other infections, drug reactions, and other causes of acute lung injury”|\n| Indeterminate appearance | Absence of typical findings AND Unilateral, central or upper lung predominant distribution | “Findings indeterminate for COVID-19 pneumonia and which can occur with a variety of infections and noninfectious conditions” |\n| Atypical appearance | Pneumothorax or pleural effusion Pulmonary edema Lobar consolidation Solitary lung nodule or mass Diffuse tiny nodules Cavity\t|“Findings atypical or uncommonly reported for COVID-19 pneumonia. Consider alternative diagnoses” |\n|Negative for pneumonia | No lung opacities | “No findings of pneumonia. However, chest radiographic findings can be absent early in the course of COVID-19 pneumonia”|\n      \nAlthough these findings refer to the CXR themselves, on this challenge we were provided with these labels only at the study level, while each study can have many images.\nOn image level, each image has a list of bounding boxes of findings. The bounding boxes can contain findings from different types, as described by the competition hosts:\n> Bounding boxes were placed on lung opacities, whether typical or indeterminate. Bounding boxes were also placed on some atypical findings including solitary lobar consolidation, nodules/masses, and cavities. Bounding boxes were not placed on pleural effusions, or pneumothoraces. No bounding boxes were placed for the negative for pneumonia category.\n\nThe dataset doesn't distinguish between the findings type. The findings were given the label opacity, and the prediction in the submission for the findings class should be always  `opacity`.  \nThe details of the study grading method according to the findings in the images described in the table above. Even though the exact meaning of the terminology is definitely beyond my understanding, one thing we can learn from this table is that the classifying is based on the nature of findings, as well as on their region in the lungs. This is crucial for a better understanding of what our model is supposed to learn.","metadata":{}},{"cell_type":"markdown","source":"## Evalutaion\n\nThe evaluation at the study level is quite simple: we will check our prediction accuracy. But at the image level, our predictions will be bounding boxes. Probably the bounding boxes will not match exactly to the labeled ones, and we do not care if there are minor differences. So how will we decide whether our predictions are consistent with the labels or not?\nCome to think of it, the most important thing here is how much our predictions area intersect with the ground truth labels. For ideal prediction, the predicted area will match the ground truth exactly. Meaning, the intersection and the area of each prediction and the ground truth are equals. In a more realistic case, our prediction isa bit smaller or larger than the ground truth, or span out in one direction and too short in another. In all these cases, the more the intersection area is large with respect to both the predicted area and the ground truth label area, the more we can regard the prediction as correct. This is the rationale behind the PASCAL VOC2010 IoU (Intersection over Union) evaluation method, which is used in this competition: a bounding box prediction is considered correct if the rate of the intersection over union of  the prediction and the ground truth is greater than $0.5$. i.e, we demand\n$$IoU = \\frac{A_y \\cap A_\\hat{y}}{A_y \\cup A_\\hat{y}} > 0.5$$  \nWhere $A_y$ is the ground truth bounding box area and $A_\\hat{y}$ is the area of the predicted box. ","metadata":{}},{"cell_type":"markdown","source":"# Notebook initialize","metadata":{}},{"cell_type":"markdown","source":"## Basic Imports","metadata":{}},{"cell_type":"code","source":"!pip install numpy --upgrade\n!pip install python-gdcm","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-06-06T22:05:03.34935Z","iopub.execute_input":"2021-06-06T22:05:03.349801Z","iopub.status.idle":"2021-06-06T22:05:03.357168Z","shell.execute_reply.started":"2021-06-06T22:05:03.349756Z","shell.execute_reply":"2021-06-06T22:05:03.355891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pathlib import Path\nimport sys\nfrom ast import literal_eval\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib as mlp\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\n\nimport seaborn as sns\n\nprint(sys.version)\n\nmlp.rcParams['figure.figsize'] = (15, 7)","metadata":{"id":"ndzKZLIlqLVU","executionInfo":{"status":"ok","timestamp":1621713123565,"user_tz":-180,"elapsed":983,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"execution":{"iopub.status.busy":"2021-06-06T22:05:14.616481Z","iopub.execute_input":"2021-06-06T22:05:14.616981Z","iopub.status.idle":"2021-06-06T22:05:15.36612Z","shell.execute_reply.started":"2021-06-06T22:05:14.616943Z","shell.execute_reply":"2021-06-06T22:05:15.36506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Enviroment Settings","metadata":{}},{"cell_type":"code","source":"is_colab = 'google.colab' in sys.modules\nif is_colab:\n    from google.colab import drive\n    drive.mount('/content/drive')\n    path = Path('/content/drive/MyDrive/covid19-detection/data') \nelse:\n    path = Path('/kaggle/input/siim-covid19-detection') ","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:05:32.251529Z","iopub.execute_input":"2021-06-06T22:05:32.252312Z","iopub.status.idle":"2021-06-06T22:05:32.259168Z","shell.execute_reply.started":"2021-06-06T22:05:32.252255Z","shell.execute_reply":"2021-06-06T22:05:32.258071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EDA","metadata":{}},{"cell_type":"markdown","source":"The dataset is composed of three parts. The CRX files in DICOM format, and two metadata tables: one for the image level and another for the study level. Let's explore first the image-level metadata of the training set.","metadata":{"id":"ZqZmqPidqXjM"}},{"cell_type":"code","source":"image_df = pd.read_csv(path/'train_image_level.csv', index_col='id')","metadata":{"id":"Qr5i_4VLqR__","executionInfo":{"status":"ok","timestamp":1621713136142,"user_tz":-180,"elapsed":890,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"execution":{"iopub.status.busy":"2021-06-06T22:05:35.913905Z","iopub.execute_input":"2021-06-06T22:05:35.914551Z","iopub.status.idle":"2021-06-06T22:05:35.979762Z","shell.execute_reply.started":"2021-06-06T22:05:35.914496Z","shell.execute_reply":"2021-06-06T22:05:35.97844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_df.head()","metadata":{"id":"OL-6oimSqiIm","executionInfo":{"status":"ok","timestamp":1621713143097,"user_tz":-180,"elapsed":303,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"f602bceb-0acf-480f-fb37-14e6244f5413","execution":{"iopub.status.busy":"2021-06-06T22:05:38.811784Z","iopub.execute_input":"2021-06-06T22:05:38.812169Z","iopub.status.idle":"2021-06-06T22:05:38.837427Z","shell.execute_reply.started":"2021-06-06T22:05:38.812139Z","shell.execute_reply":"2021-06-06T22:05:38.836402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Before doing anything else, we'd like to change this terrible column name `StudyInstanceUID` to a more reasonable one.","metadata":{}},{"cell_type":"code","source":"image_df = image_df.rename(columns={'StudyInstanceUID':'study_id'})","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:20:29.59206Z","iopub.execute_input":"2021-06-06T22:20:29.59245Z","iopub.status.idle":"2021-06-06T22:20:29.60007Z","shell.execute_reply.started":"2021-06-06T22:20:29.592419Z","shell.execute_reply":"2021-06-06T22:20:29.59902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now it's much better.\n\nFor each image we are provided with image id, study id, the findings bounding boxes, and labels for each bounding box. Let's examine the label column first.\nThe content of the label column corresponds to the submission's desired format. It contains a description of an unlimited number of finding, separated by whitespace. Each of the descriptions contains 6 fields, also separated by whitespace, as follows:\n  \n`finding_label confidence xmin ymin xmax ymax  `\nThe content of this row is as this pattern, repeated as the number of this image's findings. So if we have $k$ finding for specific image, the label row will be:\n`finding_label_1 confidence_1 xmin_1 ymin_1 xmax_1 ymax_1 finding_label_2 ...  finding_label_k confidence_k xmin_k ymin_k xmax_k ymax_k `\n\nLet's extract these values.","metadata":{"id":"E6226mLzq1Ut"}},{"cell_type":"code","source":"def extract_label(row):\n    values = row.label.split()\n    if len(values) % 6 != 0:\n        # corrupted row\n        print(f'row #{row.index}: wrong number of paramerers in label field')\n    return [dict(zip('id finding_id label confident xmin ymin xmax ymax'.split(), [row.name, i] + values[6*i:6*(i+1)])) for i in range(len(values) // 6)]","metadata":{"id":"kJBYLM2SrYDm","executionInfo":{"status":"ok","timestamp":1621713148402,"user_tz":-180,"elapsed":343,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"execution":{"iopub.status.busy":"2021-06-06T22:07:16.891446Z","iopub.execute_input":"2021-06-06T22:07:16.89192Z","iopub.status.idle":"2021-06-06T22:07:16.898806Z","shell.execute_reply.started":"2021-06-06T22:07:16.891883Z","shell.execute_reply":"2021-06-06T22:07:16.897522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"findings = pd.DataFrame.from_dict(image_df.apply(extract_label, axis=1).sum()).set_index('id')\nfindings.head()","metadata":{"id":"flFjvZWABw-F","executionInfo":{"status":"ok","timestamp":1621713150597,"user_tz":-180,"elapsed":423,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"a22859d1-cb5f-4587-e8e9-69422921b5a5","execution":{"iopub.status.busy":"2021-06-06T22:07:31.32519Z","iopub.execute_input":"2021-06-06T22:07:31.325752Z","iopub.status.idle":"2021-06-06T22:07:31.668767Z","shell.execute_reply.started":"2021-06-06T22:07:31.325701Z","shell.execute_reply":"2021-06-06T22:07:31.667764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can see the domains of these values","metadata":{}},{"cell_type":"code","source":"print(f'The unique findings label values are {findings.label.unique()}')","metadata":{"id":"gf-fRXPyCKur","executionInfo":{"status":"ok","timestamp":1621713153038,"user_tz":-180,"elapsed":422,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"c885d40e-cbb2-4806-cdf1-2f9ac7e876e9","execution":{"iopub.status.busy":"2021-06-06T22:11:38.668115Z","iopub.execute_input":"2021-06-06T22:11:38.668546Z","iopub.status.idle":"2021-06-06T22:11:38.677373Z","shell.execute_reply.started":"2021-06-06T22:11:38.66851Z","shell.execute_reply":"2021-06-06T22:11:38.676072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'The unique findings confdence values are {findings.confident.unique()}')","metadata":{"id":"_XsSyVXNCvw5","executionInfo":{"status":"ok","timestamp":1621713155197,"user_tz":-180,"elapsed":5,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"6b9d2339-3182-4d12-c5d0-7b679cf4efef","execution":{"iopub.status.busy":"2021-06-06T22:11:42.065782Z","iopub.execute_input":"2021-06-06T22:11:42.066193Z","iopub.status.idle":"2021-06-06T22:11:42.074207Z","shell.execute_reply.started":"2021-06-06T22:11:42.066158Z","shell.execute_reply":"2021-06-06T22:11:42.072935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The labels are only `none` and `opacity`, and the confidence in the training set is always 1 (since this is a labeled dataset). All bounding box data is provided in the `boxes` field, this field is `Nan` when there are no findings (as we can see in the second row in the head of the Dataframe printed above). So in fact, all the data we need exists in the `boxes` field.  \nThus, we can extract the findings data directly from the boxes fields and examine some of their properties.","metadata":{"id":"Qf4UaR6-C01i"}},{"cell_type":"code","source":"findings = image_df.apply(lambda x: [{'id':x.name, 'finding_id': i, **box} for i, box in enumerate(literal_eval(x.boxes))] if type(x.boxes) == str else [{'id': x.name}], axis=1)\nfindings = pd.DataFrame.from_dict(findings.sum()).set_index('id')","metadata":{"id":"PoYbLAM7ROG1","executionInfo":{"status":"ok","timestamp":1621713156863,"user_tz":-180,"elapsed":999,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"e95cf2b8-c928-4e36-f584-c6df1ec42f14","execution":{"iopub.status.busy":"2021-06-06T22:15:40.368043Z","iopub.execute_input":"2021-06-06T22:15:40.368416Z","iopub.status.idle":"2021-06-06T22:15:40.891026Z","shell.execute_reply.started":"2021-06-06T22:15:40.368385Z","shell.execute_reply":"2021-06-06T22:15:40.88983Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'Total number of findings: {findings.shape[0]}')\nfindings.head()","metadata":{"id":"_GsRb9K0RN__","executionInfo":{"status":"ok","timestamp":1621713156864,"user_tz":-180,"elapsed":15,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"5ce06293-c3bb-4bfe-8a94-b3ba47b92b3f","execution":{"iopub.status.busy":"2021-06-06T22:15:44.311642Z","iopub.execute_input":"2021-06-06T22:15:44.312031Z","iopub.status.idle":"2021-06-06T22:15:44.32878Z","shell.execute_reply.started":"2021-06-06T22:15:44.311999Z","shell.execute_reply":"2021-06-06T22:15:44.327651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Before further exploration of findings properties, let's explore the study-level metadata:","metadata":{"id":"KR3oQnrgUEeU"}},{"cell_type":"code","source":"study_df = pd.read_csv(path/'train_study_level.csv', index_col='id')","metadata":{"id":"nfxdSSm5UlK3","executionInfo":{"status":"ok","timestamp":1621713157335,"user_tz":-180,"elapsed":479,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"execution":{"iopub.status.busy":"2021-06-06T22:16:08.07106Z","iopub.execute_input":"2021-06-06T22:16:08.071822Z","iopub.status.idle":"2021-06-06T22:16:08.09488Z","shell.execute_reply.started":"2021-06-06T22:16:08.071782Z","shell.execute_reply":"2021-06-06T22:16:08.093869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"study_df.head()","metadata":{"id":"bDwtYcZPUl1b","executionInfo":{"status":"ok","timestamp":1621713157341,"user_tz":-180,"elapsed":14,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"a1d4b7b5-31cd-450a-d4fa-026f0351b4d1","execution":{"iopub.status.busy":"2021-06-06T22:16:09.413406Z","iopub.execute_input":"2021-06-06T22:16:09.413797Z","iopub.status.idle":"2021-06-06T22:16:09.425279Z","shell.execute_reply.started":"2021-06-06T22:16:09.413764Z","shell.execute_reply":"2021-06-06T22:16:09.424291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this dataframe each study is classified into one of 4 classes: `Negative for Pneumonia`, `Typical Appearance`, `Indeterminate Appearance`, and `Atypical Appearance`. It is important to know how these classes are distributed over the dataset.  \nIn the evaluation section in the competition details in Kaggle it's said that\n> Studies in the test set may contain more than one label. They are as follows: `negative`, `typical`, `indeterminate`, `atypical`\n\nAccordingly,  this is a multilabel classification task.    \nIn contrast, in [a post](https://www.kaggle.com/c/siim-covid19-detection/discussion/240250) in the competition discussion section, the hosts indicated that\n> Per the grading schema, chest radiographs are classified into one of four categories, which are **mutually exclusive**\n\nSince the two descriptions contradict, it is worth inspecting the training set to see which labels can be assigned to an image together.","metadata":{"id":"0URL3Dc62NsR"}},{"cell_type":"code","source":"(study_df.apply(sum, axis=1) == 1).all()","metadata":{"id":"nDFxpTQi6-Xy","executionInfo":{"status":"ok","timestamp":1621713204161,"user_tz":-180,"elapsed":472,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"56604cc3-6a20-4a16-8e92-40e2438d0294","execution":{"iopub.status.busy":"2021-06-06T22:18:59.529273Z","iopub.execute_input":"2021-06-06T22:18:59.529702Z","iopub.status.idle":"2021-06-06T22:18:59.624058Z","shell.execute_reply.started":"2021-06-06T22:18:59.529667Z","shell.execute_reply":"2021-06-06T22:18:59.622952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So in the training set, each of the studies has a single label attached to it, and this classification is in fact one-hot encoded classification for each study to one of those 4 classes.  \n\nSince the results of our check on the training set supports the second post, and since inherently, by their meanings, the labels seem to be mutually exclusive , we will leave it as a single-label classification task. For convenience , we will store the labels in one row rather than in one-hot encoding format and join it with the images dataframe.","metadata":{"id":"1t0uD9V0Ulth"}},{"cell_type":"code","source":"image_df.study_id += '_study' \nstudy_labels = study_df.idxmax(axis=1).rename('study_label')\n# study_labels.index = study_labels.index.str.extract('(^[^_]*)').apply(lambda x:x[0], axis=1)\nimage_df = image_df.merge(study_labels, left_on='study_id', right_index=True)","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:20:33.980091Z","iopub.execute_input":"2021-06-06T22:20:33.980511Z","iopub.status.idle":"2021-06-06T22:20:34.022744Z","shell.execute_reply.started":"2021-06-06T22:20:33.980474Z","shell.execute_reply":"2021-06-06T22:20:34.021646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next we will check how these classes are distributed over the dataset.","metadata":{"id":"ik856etENKHy"}},{"cell_type":"code","source":"label_count = study_df.sum()\nplt.figure(figsize=(8, 8))\nplt.pie(label_count, labels=label_count.index, wedgeprops={'edgecolor': 'black'}, autopct='%1.f%%', textprops={'fontsize': 16}, explode=[.01]*4, shadow=True)\nplt.title('Label Distribution', fontdict={'fontsize':22});","metadata":{"id":"j-sJPhSkLCHD","executionInfo":{"status":"ok","timestamp":1621713212198,"user_tz":-180,"elapsed":549,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"a61de634-c487-4027-8c6d-c5271f0aac0a","execution":{"iopub.status.busy":"2021-06-06T22:20:49.173279Z","iopub.execute_input":"2021-06-06T22:20:49.173716Z","iopub.status.idle":"2021-06-06T22:20:49.376502Z","shell.execute_reply.started":"2021-06-06T22:20:49.173676Z","shell.execute_reply":"2021-06-06T22:20:49.375725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have 47% certain Covid cases (typical appearance), 36% non-covid (28% negative for pneumonia and 8% atypical to covid), and 17% obscure cases. From a covid vs non-covid point of view, the dataset is quite balanced. But from the classification point of view,  almost 50% of the cases are from one class and only 8% of the cases are from another class.  \nTo better understand the distribution, let’s see that classification distribution in absolute numbers:","metadata":{"id":"PGSZw5wdOW7j"}},{"cell_type":"code","source":"g = sns.barplot(x=label_count.index, y=label_count)\nfor label, count in zip(range(label_count.index.shape[0]), label_count):\n    g.text(label, count, count, color='black', ha=\"center\", fontdict=dict(fontsize=14, weight='bold'))","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:21:20.471399Z","iopub.execute_input":"2021-06-06T22:21:20.471997Z","iopub.status.idle":"2021-06-06T22:21:20.658417Z","shell.execute_reply.started":"2021-06-06T22:21:20.471959Z","shell.execute_reply":"2021-06-06T22:21:20.65749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, let's look at some properties of the findings. The number of findings varies for each image. Each of them is an opacity (or another type of the above-mentioned findings) in the CXR, and for each of them we provided the bounding box of the opacity area. Let’s look at the number of findings per image and the main statistical properties of their areas: sum, mean, max, etc.","metadata":{}},{"cell_type":"code","source":"findings['area'] = findings.width * findings.height\nfindings_props = findings.groupby('id').area.agg(['count','sum', 'min', 'max', 'mean', 'std'])","metadata":{"id":"Y4Glmi06Q5My","executionInfo":{"status":"ok","timestamp":1621713215480,"user_tz":-180,"elapsed":389,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"execution":{"iopub.status.busy":"2021-06-06T22:21:43.64396Z","iopub.execute_input":"2021-06-06T22:21:43.644408Z","iopub.status.idle":"2021-06-06T22:21:43.668486Z","shell.execute_reply.started":"2021-06-06T22:21:43.644366Z","shell.execute_reply":"2021-06-06T22:21:43.667559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# findings_props.index = findings_props.index.str.extract('(^[^_]*)').apply(lambda x:x[0], axis=1)\nimage_props = image_df.join(findings_props)\nimage_props.head()","metadata":{"id":"5JqABYuxRsqz","executionInfo":{"status":"ok","timestamp":1621713217073,"user_tz":-180,"elapsed":5,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"a3d3409f-8c7e-4523-8562-d6acb92f98c8","execution":{"iopub.status.busy":"2021-06-06T22:21:44.896859Z","iopub.execute_input":"2021-06-06T22:21:44.897497Z","iopub.status.idle":"2021-06-06T22:21:44.927799Z","shell.execute_reply.started":"2021-06-06T22:21:44.897459Z","shell.execute_reply":"2021-06-06T22:21:44.926691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The findings count for each class label is:","metadata":{"execution":{"iopub.status.busy":"2021-05-31T23:00:56.910211Z","iopub.execute_input":"2021-05-31T23:00:56.910549Z","iopub.status.idle":"2021-05-31T23:00:56.918248Z","shell.execute_reply.started":"2021-05-31T23:00:56.91052Z","shell.execute_reply":"2021-05-31T23:00:56.916473Z"}}},{"cell_type":"code","source":"plt.figure(figsize=(10,6))\nplt.gca().yaxis.grid(linestyle='--')\nsns.violinplot(data=image_props, x='study_label', y='count')\nplt.title('Findings Count', fontdict={'fontsize':16})\nplt.show()","metadata":{"id":"WpBFBm6ebEx1","executionInfo":{"status":"ok","timestamp":1621713219835,"user_tz":-180,"elapsed":425,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"f58cf779-149b-4c83-a1c9-defdc4c0c2ac","execution":{"iopub.status.busy":"2021-06-06T22:22:03.769713Z","iopub.execute_input":"2021-06-06T22:22:03.770139Z","iopub.status.idle":"2021-06-06T22:22:04.013315Z","shell.execute_reply.started":"2021-06-06T22:22:03.770099Z","shell.execute_reply":"2021-06-06T22:22:04.012444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is clear now that all the negative results have no findings at all, as stated in the grading method table. On the other hand, for each of the other four types it seems that there are instances with no opacity findings, contrary to these grades descriptions in the above table. But is this really the case? We saw earlier that we have more images than studied. That is, some studies have more than one image. So it seems that in some cases the prognosis is based on findings that are determined only in one of the scans. Let's verify this conclusion.","metadata":{"id":"am-ljE10qyp6"}},{"cell_type":"code","source":"study_findings_props = image_props.groupby(['study_id'])\nstudy_findings_count = study_findings_props['count'].agg('sum').to_frame().join(study_labels)\nplt.figure(figsize=(10,6))\nplt.gca().yaxis.grid(linestyle='--')\nsns.violinplot(data=study_findings_count, x='study_label', y='count')\nplt.title('Findings Count', fontdict={'fontsize':16})\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:22:32.840955Z","iopub.execute_input":"2021-06-06T22:22:32.84143Z","iopub.status.idle":"2021-06-06T22:22:33.263051Z","shell.execute_reply.started":"2021-06-06T22:22:32.84139Z","shell.execute_reply":"2021-06-06T22:22:33.261951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It becomes clear that the above conclusion is correct. Almost all of the clear covid-19 cases have 2 findings, and a couple of them with 3 findings. The indeterminate cases also have at least one finding each, and only the non-covid cases sometimes have no findings, even when positive to pneumonia. But according to our table, `No Findings` means `Positive to Pneumonia`, so we'll put these instances aside for now.","metadata":{}},{"cell_type":"code","source":"studies_to_remove = study_df[study_df.index.isin(study_findings_props['count'].sum()[study_findings_props['count'].sum() == 0].index) & \n         (study_df['Negative for Pneumonia'] == 0)]\nstudy_df['removed'] = False\nstudy_df.loc[studies_to_remove.index, 'removed'] = True\nprint(f'Total number of removed rows: {study_df.loc[studies_to_remove.index].shape[0]}')\nprint(f'\\n\\nRemoved rows by label:\\n')\nprint(study_df.loc[studies_to_remove.index].iloc[:, :-1].sum().to_string())","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:23:32.274643Z","iopub.execute_input":"2021-06-06T22:23:32.275072Z","iopub.status.idle":"2021-06-06T22:23:32.299019Z","shell.execute_reply.started":"2021-06-06T22:23:32.275038Z","shell.execute_reply":"2021-06-06T22:23:32.297085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_props = image_props.drop(image_props.loc[image_props.study_id.isin(studies_to_remove.index)].index)\n\nstudy_findings_props = image_props.groupby(['study_id'])\nstudy_findings_count = study_findings_props['count'].agg('sum').to_frame().join(study_labels)\nplt.figure(figsize=(10,6))\n# plt.gca().yaxis.grid(linestyle='--')\nsns.violinplot(data=study_findings_count, x='study_label', y='count')\nplt.title('Findings Count', fontdict={'fontsize':16})\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:23:34.70818Z","iopub.execute_input":"2021-06-06T22:23:34.708832Z","iopub.status.idle":"2021-06-06T22:23:34.967595Z","shell.execute_reply.started":"2021-06-06T22:23:34.708773Z","shell.execute_reply":"2021-06-06T22:23:34.96673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now all the positive cases have findings. Let's inspect other findings properties:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(20,10))\n\nfor i, prop in enumerate(['sum', 'min', 'max', 'mean', 'std'], start=1):\n    plt.subplot(2, 3, i)\n#     plt.gca().yaxis.grid(linestyle='--')\n    sns.violinplot(data=image_props, x='study_label', y=prop, order=study_labels.unique())\n    plt.xticks(rotation=10)\n    title = f'Findings {prop.title()}'\n    if prop.title() != 'Sum': title += ' Area'\n    plt.title(title, fontdict={'fontsize':16})\nplt.tight_layout()\nplt.show()","metadata":{"id":"0t8PYKKwupkt","executionInfo":{"status":"ok","timestamp":1621713225209,"user_tz":-180,"elapsed":1534,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"67c05f6c-87cc-4b66-e24a-c1a05ec1dde3","execution":{"iopub.status.busy":"2021-06-06T22:24:50.729578Z","iopub.execute_input":"2021-06-06T22:24:50.73017Z","iopub.status.idle":"2021-06-06T22:24:51.904467Z","shell.execute_reply.started":"2021-06-06T22:24:50.730118Z","shell.execute_reply":"2021-06-06T22:24:51.903622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One can see that clear covid cases strongly tend to have  larger findings area. The indeterminate cases also tend to have larger findings areas than the atypical ones, but this difference is much less significant. Let's inspect these features again, but now at the study level.","metadata":{"id":"v-9rwN79wBPc"}},{"cell_type":"code","source":"study_props = (image_props[['study_id']]\n               .join(findings)\n               .groupby('study_id').area.agg(['count','sum', 'min', 'max', 'mean', 'std'])\n               .merge(study_labels, left_on='study_id', right_index=True))\n\nplt.figure(figsize=(20,10))\n\nfor i, prop in enumerate(['sum', 'min', 'max', 'mean', 'std'], start=1):\n    plt.subplot(2, 3, i)\n#     plt.gca().yaxis.grid(linestyle='--')\n    sns.violinplot(data=study_props, x='study_label', y=prop, order=study_labels.unique(), pallete=['blue', 'green', 'green', 'red'])\n    plt.xticks(rotation=10)\n    title = f'Findings {prop.title()}'\n    if prop.title() != 'Sum': title += ' Area'\n    plt.title(title, fontdict={'fontsize':16})\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:25:05.172185Z","iopub.execute_input":"2021-06-06T22:25:05.172872Z","iopub.status.idle":"2021-06-06T22:25:06.283903Z","shell.execute_reply.started":"2021-06-06T22:25:05.172818Z","shell.execute_reply":"2021-06-06T22:25:06.282642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It seems that there is no significant difference between the image level and study level. But, this leads us to two important questions: How many images are related to one study on average? In case that a study had more than one image, how many images the prognosis is based on?","metadata":{}},{"cell_type":"code","source":"image_df.groupby('study_id').label.count().unique()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:25:54.471435Z","iopub.execute_input":"2021-06-06T22:25:54.471906Z","iopub.status.idle":"2021-06-06T22:25:54.491114Z","shell.execute_reply.started":"2021-06-06T22:25:54.471864Z","shell.execute_reply":"2021-06-06T22:25:54.489852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(image_df.groupby('study_id').label.agg(images_count='count').value_counts().to_string())\n# g = sns.countplot(data=image_df.groupby('study_id').label.agg(images_count='count'), x='images_count')","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:25:55.989213Z","iopub.execute_input":"2021-06-06T22:25:55.989657Z","iopub.status.idle":"2021-06-06T22:25:56.013493Z","shell.execute_reply.started":"2021-06-06T22:25:55.989618Z","shell.execute_reply":"2021-06-06T22:25:56.012392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In most cases there’s one image per study. But in the cases with multiple images, what is the difference between the images? Is the prognosis based on all of the images?  \nIt is straightforward to get an answer to the second password from the data. We simply count the number of images labeled with findings.","metadata":{}},{"cell_type":"code","source":"multiple_images = image_df.groupby('study_id').filter((lambda x: x.label.count() > 1))\nmultiple_images_with_findings = multiple_images[multiple_images.boxes.notna()]\nprint('Images with finding in Study Counts')\nprint(multiple_images_with_findings.groupby('study_id').boxes.agg(count='count').value_counts().to_string())","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:26:26.391787Z","iopub.execute_input":"2021-06-06T22:26:26.392199Z","iopub.status.idle":"2021-06-06T22:26:27.465093Z","shell.execute_reply.started":"2021-06-06T22:26:26.392162Z","shell.execute_reply":"2021-06-06T22:26:27.463978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So it is established that , there is never more than one image labeled with findings per study. Now it will interesting to see the difference between the images in one study. To do that, we have to pay attention to the third and most important part of our dataset - the DICOM files.","metadata":{}},{"cell_type":"markdown","source":"#### DICOM files\nThe data is provided in [DICOM fromat](https://en.wikipedia.org/wiki/DICOM), which is the standard in medical imaging information and related data. This format packs each medical image with related data, such as Patient Id, Name, Sex, etc. In our case, the data de-identified for privacy reasons, but we still may have important data in the metadata provided in the DICOM file.\nLet's pick a file and see what it looks like.","metadata":{"id":"d9ZxtKX41Ixu"}},{"cell_type":"code","source":"from pydicom import dcmread\n\ndef extract_id(full_id): return full_id[:full_id.index('_')]\ndef get_file_path(study_id, image_id, dataset='train'):\n    study_id, image_id = extract_id(study_id), extract_id(image_id)\n    return [*(path/dataset/study_id).glob(f'**/{image_id}.dcm')][0] \n\n#     return [*(path/dataset/row.study_id.str.extract(\"(^[^_]*)\").values[0,0]).glob(f'**/{row.index.str.extract(\"(^[^_]*)\").values[0,0]}.dcm')][0] \nsample = image_df.sample(random_state=14)\nfpath = get_file_path(sample.study_id.values[0], sample.index.values[0])\nds = dcmread(fpath)\n\nprint(ds)\nplt.imshow(ds.pixel_array, cmap=plt.cm.gray)\nplt.colorbar()\nplt.show()","metadata":{"id":"K90QUFlb1nSi","executionInfo":{"status":"ok","timestamp":1621713238518,"user_tz":-180,"elapsed":426,"user":{"displayName":"איתמר דביר","photoUrl":"","userId":"06992624979939551967"}},"outputId":"56e7be2d-8530-446d-f9c6-90e17225e8ef","execution":{"iopub.status.busy":"2021-06-06T22:28:42.383792Z","iopub.execute_input":"2021-06-06T22:28:42.384215Z","iopub.status.idle":"2021-06-06T22:28:44.624885Z","shell.execute_reply.started":"2021-06-06T22:28:42.384173Z","shell.execute_reply":"2021-06-06T22:28:44.62353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a taste from our data:","metadata":{}},{"cell_type":"code","source":"# import gdcm\ndef show_dicoms(df, ncols=5, size=5, annotate=False):\n    n = df.shape[0]\n    nrows = int(np.ceil(n/ncols))\n\n    fig = plt.figure(figsize=(size*ncols, size*nrows))\n    for i, row in enumerate(df.itertuples()):\n        fpath = get_file_path(row.study_id, row.Index)\n        ds = dcmread(fpath)\n        plt.subplot(nrows, ncols,i+1)\n        plt.imshow(ds.pixel_array, cmap=plt.cm.gray)\n        if annotate:\n            if isinstance(row.boxes, str):\n                for box in literal_eval(row.boxes):\n                    rect = patches.Rectangle((box['x'], box['y']), \n                                             box['width'], box['height'],\n                                            color='r', fill=False)\n                    plt.gca().add_patch(rect)\n        plt.xticks([]), plt.yticks([])\n    plt.subplots_adjust()\n    return fig\n#     plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:28:51.872758Z","iopub.execute_input":"2021-06-06T22:28:51.873205Z","iopub.status.idle":"2021-06-06T22:28:51.882996Z","shell.execute_reply.started":"2021-06-06T22:28:51.873165Z","shell.execute_reply":"2021-06-06T22:28:51.882175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_dicoms(image_df.sample(25, random_state=25))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:32:14.416813Z","iopub.execute_input":"2021-06-06T22:32:14.417245Z","iopub.status.idle":"2021-06-06T22:32:38.962973Z","shell.execute_reply.started":"2021-06-06T22:32:14.417208Z","shell.execute_reply":"2021-06-06T22:32:38.962034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Many of the images are cropped, rotated, and have different lummination level. Lungs are contained in all of the images, but the location of the lungs in the image is not constant, The images margin size are varying, and the images may contain other body parts - neck, stomach, hands, etc. To get a better understanding of the matter in hand, it will be helpful to see CXR from the different labels with the annotated bounding boxes drawn on the image.","metadata":{}},{"cell_type":"code","source":"gb = image_df.groupby('study_label')\nfor name, group in gb:\n    if not name.startswith('Neg'): group = group.loc[group.boxes.notna()]\n    fig = show_dicoms(group.sample(9, random_state=4), 3, 8, annotate=True)\n    fig.suptitle(name, fontsize=18)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:33:26.961389Z","iopub.execute_input":"2021-06-06T22:33:26.962116Z","iopub.status.idle":"2021-06-06T22:34:16.21777Z","shell.execute_reply.started":"2021-06-06T22:33:26.962072Z","shell.execute_reply":"2021-06-06T22:34:16.216837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The results are not clear for the non-expert eye. Although sometimes there is a kind of opacity in the boxes, in other cases there is no clear difference between the area inside and the area outside the box. This will affect the algorithm development and verification processes since I will not be able to rely on my own knowledge and intuition.","metadata":{}},{"cell_type":"markdown","source":"Now let's take a look at the metadata provided by the DICOMs. The attributes that may interest us are the body part examined, sex, image size, pixel spacing (represent the physical size of the image), modality (the scanning method), and image type. Let's extract these features to a pandas dataframe for later use.","metadata":{}},{"cell_type":"code","source":"from tqdm.notebook import tqdm\ndata = {}\ndef append_dcm_properties(row, props):\n#     print(row)\n    try:\n        fpath = get_file_path(row.study_id, row.Index)\n#         fpath = [*(path/'train'/row.study_id[:row.study_id.index(\"_\")]).glob(f'**/{row.Index[:row.Index.index(\"_\")]}.dcm')][0] \n        ds = dcmread(fpath)\n        data[row.Index] = {prop: getattr(ds, prop.replace(' ', '')) for prop in props}\n    except Exception as e:\n        print(f'**/{row.name[:row.name.index(\"_\")]}.dcm')\n        raise\n        \nprops = ['Image Type', 'Modality','Body Part Examined', 'Photometric Interpretation',\n                                       'Patient Sex', 'Imager Pixel Spacing', 'Rows', 'Columns']\n# image_df.apply(lambda x: append_dcm_properties(x, props) , axis=1)\nfor row in tqdm(image_df.itertuples(), total=image_df.shape[0]):\n    append_dcm_properties(row, props)","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:35:46.964847Z","iopub.execute_input":"2021-06-06T22:35:46.965246Z","iopub.status.idle":"2021-06-06T22:53:33.416632Z","shell.execute_reply.started":"2021-06-06T22:35:46.965213Z","shell.execute_reply":"2021-06-06T22:53:33.415408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dm = pd.DataFrame.from_dict(data, orient='index')\ndm.to_csv('DICOM_metadata.csv')\ndm.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:54:46.825609Z","iopub.execute_input":"2021-06-06T22:54:46.826063Z","iopub.status.idle":"2021-06-06T22:54:46.875071Z","shell.execute_reply.started":"2021-06-06T22:54:46.826021Z","shell.execute_reply":"2021-06-06T22:54:46.873914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's inspect the value ranges of our new data:","metadata":{}},{"cell_type":"code","source":"for column in dm[:]:\n    print(f'{column}:' )\n    print(dm[column].unique())\n    print('-'*100)","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:55:32.69995Z","iopub.execute_input":"2021-06-06T22:55:32.700324Z","iopub.status.idle":"2021-06-06T22:55:32.720843Z","shell.execute_reply.started":"2021-06-06T22:55:32.700292Z","shell.execute_reply":"2021-06-06T22:55:32.719338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The first thing to inspect is the sex field.\nHow is our data splitted between the sexes? How the sex is related with the COVID19 prognoses?","metadata":{}},{"cell_type":"code","source":"sns.countplot(data=dm, x='Patient Sex');","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:55:41.549875Z","iopub.execute_input":"2021-06-06T22:55:41.55024Z","iopub.status.idle":"2021-06-06T22:55:41.698563Z","shell.execute_reply.started":"2021-06-06T22:55:41.55021Z","shell.execute_reply":"2021-06-06T22:55:41.69753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,7))\nsns.countplot(data=dm.join(image_df), x='study_label', hue='Patient Sex')","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:55:43.850808Z","iopub.execute_input":"2021-06-06T22:55:43.851214Z","iopub.status.idle":"2021-06-06T22:55:44.077077Z","shell.execute_reply.started":"2021-06-06T22:55:43.851179Z","shell.execute_reply":"2021-06-06T22:55:44.076362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see here the well-known fact that statistically women suffer less from Covid-19. Although our data is quite balanced with respect to sex, women suffer much less from pneumonia of any kind. In the typical covid cases (clear/severe covid cases) there are about two-thirds cases of women than men.","metadata":{}},{"cell_type":"markdown","source":"In `Body Part Examined` column, we have the unique values\n>`'CHEST' 'PORT CHEST' 'TORAX' nan 'T?RAX' 'Pecho' 'THORAX' 'ABDOMEN'\n 'SKULL' '2- TORAX' 'TÒRAX' 'PECHO'`\n \n'CHEST', 'THORAX' (which exists in many versions), and 'PECHO' are all the same, whether you prefer English or Spanish. Let's try to see what else we can extract from the metadata. (Why do we have SKULLs here???)","metadata":{}},{"cell_type":"code","source":"body_parts = dm.loc[~dm['Body Part Examined'].isin(['CHEST', 'TORAX','THORAX', 'T?RAX',\n                                                     'TÒRAX','PECHO', 'Pecho'])]\ngb = body_parts.join(image_df).groupby('Body Part Examined')\nfor name, group in gb:\n    fig = show_dicoms(group.sample(6, random_state=4), 3, 8)\n    fig.suptitle(name, fontsize=18)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-06T22:56:02.262016Z","iopub.execute_input":"2021-06-06T22:56:02.262717Z","iopub.status.idle":"2021-06-06T22:56:36.699864Z","shell.execute_reply.started":"2021-06-06T22:56:02.262663Z","shell.execute_reply":"2021-06-06T22:56:36.698991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The `ABDOMEN` images seem to contain the body's lower part too. Besides that, there does not seem to be a significant difference between the images group (in particular, we have no `SKULL`s here).\n","metadata":{}}]}