{"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":"# 1. Introduction\n\n
\n

1.1 Background

\n
\n\nWelcome to this competition on predicting **Cervical Spine Fractures** from **Computed Tomography (CT) scans**. The **motivation** to use AI for this task is that a **quick** diagnosis can **reduce** the chance of **neurologic deterioration** and paralysis after trauma.\n\n
\n\n
\n\n**Dataset origin**\n\nThe dataset we are using is made up of roughly **3000 CT studies**, [from twelve locations and across six continents](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/overview/acknowledgements). Spine **radiology specialists** have provided **annotations** to indicate the presence, vertebral level and location of any cervical spine fractures.\n\nSpecial thanks to the competition **hosts** for providing such a comprehensive dataset:\n* *Radiological Society of North America (RSNA)*\n* *American Society of Neuroradiology (ASNR)*\n* *American Society of Spine Radiology (ASSR)*\n\n
\n

1.2 Evaluation metric

\n
\n\nWe need to predict the **probability of fracture** for each of the **seven cervical vertebrae** denoted by C1, C2, C3, C4, C5, C6 and C7 as well as an **overall probability** of any fractures in the cervical spine. This means there will be **8 rows per image id** in the submission file. Note that fractures in the skull base, thoracic spine, ribs, and clavicles are **ignored**.\n \nThe **competition metric** is a **weighted multi-label logarithmic loss** (averaged across all patients)\n \n$$\nL_{ij} = - w_j \\left(y_{ij} \\log(p_{ij}) + (1-y_{ij}) \\log(1-p_{ij}) \\right)\n$$\n\nwhere the **weights** [are given by](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/340392)\n \n$$\nw_{j} = \\begin{cases}\n1, \\qquad \\text{if vertebrae negative} \\\\\n2, \\qquad \\text{if vertebrae positive} \\\\\n7, \\qquad \\text{if patient negative} \\\\\n14, \\qquad \\text{if patient positive}\n\\end{cases}\n$$\n\n
\n
\n\n
\n
","metadata":{}},{"cell_type":"markdown","source":"Notice how **more weight** is put on **positive cases** and and the **most weight** on the **overall probability** of any fractures.\n\n
\n\n**Example:**\n\nBelow is a table going through an example of **how the evaluation metric is calculated** for a single case id. This would then be **averaged** across all case id's.\n\n
\n\n
\n\nNote: it is convention to use the **natural logarithm** (base e) for calculating the log loss. \n\nNB: it [has been suggested](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/344565) that the loss is then also **scaled** by the **sum of weights** on a by-patient basis. So for our example above, the weights were [2,1,1,1,2,1,1,14] which sums to 23 hence the loss would actually be 3.734/23 = 0.162. \n\n\n
\n

1.3 Code requirements

\n
\n\nThis is a **code competition**, which means that **submissions** are made through **notebooks**. Furthermore, the submission notebook is subject to these conditions:\n\n* **run-time** (CPU/GPU) **<= 9 hours**\n* **internet** access **disabled**\n* **external data is allowed**, including pre-trained models\n* the submission file must be named **submission.csv**\n \nNote that the **test set is hidden**, and will populated when you submit your notebook. \n\n
\n\n*Competition timeline:*\n \n* Start date - 28th July 2022\n* Finish date - 27th October 2022\n\n
\n

1.4 Libraries

\n
","metadata":{}},{"cell_type":"code","source":"# To read compressed dicom files:\n\n'''Online'''\n#! pip install python-gdcm\n#! pip install pylibjpeg pylibjpeg-libjpeg pydicom\n\n'''Offline (need to add dataset --> for-pydicom)'''\n!pip install -qU ../input/for-pydicom/python_gdcm-3.0.14-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl ../input/for-pydicom/pylibjpeg-1.4.0-py3-none-any.whl --find-links frozen_packages --no-index","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:19.562016Z","iopub.execute_input":"2022-08-18T22:10:19.562825Z","iopub.status.idle":"2022-08-18T22:10:31.239153Z","shell.execute_reply.started":"2022-08-18T22:10:19.562744Z","shell.execute_reply":"2022-08-18T22:10:31.237712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport matplotlib.patches as patches\nimport seaborn as sns\nsns.set(style='darkgrid', font_scale=1.6)\nimport cv2\nimport os\nfrom os import listdir\nimport re\nimport gc\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nfrom tqdm import tqdm\nfrom pprint import pprint\nfrom time import time\nimport itertools\nfrom skimage import measure \nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport nibabel as nib\nfrom glob import glob\nimport warnings\n#warnings.filterwarnings(\"ignore\", category=DeprecationWarning)\n#warnings.filterwarnings(\"ignore\", category=UserWarning)\n#warnings.filterwarnings(\"ignore\", category=FutureWarning)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:31.241508Z","iopub.execute_input":"2022-08-18T22:10:31.242333Z","iopub.status.idle":"2022-08-18T22:10:32.313522Z","shell.execute_reply.started":"2022-08-18T22:10:31.242295Z","shell.execute_reply":"2022-08-18T22:10:32.312432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Data\n\n
\n

2.1 Data frames

\n
\n","metadata":{}},{"cell_type":"code","source":"# Load metadata\ntrain_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\ntrain_bbox = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv\")\ntest_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/test.csv\")\nss = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/sample_submission.csv\")\n\n# Print dataframe shapes\nprint('train shape:', train_df.shape)\nprint('train bbox shape:', train_bbox.shape)\nprint('test shape:', test_df.shape)\nprint('ss shape:', ss.shape)\nprint('')\n\n# Show first few entries\ntrain_df.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:32.31468Z","iopub.execute_input":"2022-08-18T22:10:32.315005Z","iopub.status.idle":"2022-08-18T22:10:32.378176Z","shell.execute_reply.started":"2022-08-18T22:10:32.314952Z","shell.execute_reply":"2022-08-18T22:10:32.377129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> **train.csv** - contains metadata for train_images.\n> * StudyInstanceUID - The study ID. There is one unique study ID for each patient scan.\n> * patient_overall - The patient level outcome, i.e. if any of the vertebrae are fractured.\n> * C[1-7] - Whether the given vertebrae is fractured.","metadata":{}},{"cell_type":"code","source":"train_bbox.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:32.3808Z","iopub.execute_input":"2022-08-18T22:10:32.381142Z","iopub.status.idle":"2022-08-18T22:10:32.394498Z","shell.execute_reply.started":"2022-08-18T22:10:32.381112Z","shell.execute_reply":"2022-08-18T22:10:32.393415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> **train_bounding_boxes.csv** - contains bounding boxes of where fractures occured for a subset of the training set.\n> * StudyInstanceUID - The study ID. There is one unique study ID for each patient scan.\n> * x - x-coordinate of bounding box bottom left corner\n> * y - y-coordinate of bounding box bottom left corner\n> * width - width of bounding box\n> * height - height of bounding box\n> * slice_number - slice number of scan\n\nNote: we only have bounding boxes for a [subset](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/343105) of the training set. We'll explore the exact proportion later on.","metadata":{}},{"cell_type":"code","source":"test_df.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:32.396024Z","iopub.execute_input":"2022-08-18T22:10:32.396611Z","iopub.status.idle":"2022-08-18T22:10:32.407637Z","shell.execute_reply.started":"2022-08-18T22:10:32.396561Z","shell.execute_reply":"2022-08-18T22:10:32.406422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> **test.csv** - contains metadata for test_images.\n> * row_id - The row ID. This will match the same column in the sample submission file.\n> * StudyInstanceUID - The study ID.\n> * prediction_type - Which one of the eight target columns needs a prediction in this row.\n\nNote: The full test set will be **populated at inference time**.","metadata":{}},{"cell_type":"code","source":"ss.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:32.409481Z","iopub.execute_input":"2022-08-18T22:10:32.410014Z","iopub.status.idle":"2022-08-18T22:10:32.422252Z","shell.execute_reply.started":"2022-08-18T22:10:32.409957Z","shell.execute_reply":"2022-08-18T22:10:32.421473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> **sample_submission.csv** - A valid sample submission.\n> * row_id - The row ID. See the test.csv for what prediction needs to be filed in that row.\n> * fracture - The target column.","metadata":{}},{"cell_type":"markdown","source":"# 3. EDA\n
\n

3.1 Fracture distributions

\n
\n\n* The overall target is roughly **balanced** (52/48 split). \n* **C7** has the **highest proportion of fractures** (19%) whereas **C3** has the **lowest** (4%). \n* Several patients have **more than one** fracture.\n* If **multiple fractures** occur on a single patient, they tend to occur in vertebrae **close together**, e.g. C4 & C5 as opposed to C1 & C7.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(20,6))\nplt.subplot(1,2,1)\nax1 = sns.countplot(data=train_df, x='patient_overall')\nfor container in ax1.containers:\n ax1.bar_label(container)\nplt.title('Fractures by patient')\nplt.ylim([0,1300])\n\n# Unpivot train_df for plotting\ntrain_melt = pd.melt(train_df, id_vars = ['StudyInstanceUID', 'patient_overall'],\n value_vars = ['C1','C2','C3','C4','C5','C6','C7'],\n var_name=\"Vertebrae\",\n value_name=\"Fractured\")\n\nplt.subplot(1,2,2)\nax2 = sns.countplot(data=train_melt, x='Vertebrae', hue='Fractured')\nfor container in ax2.containers:\n ax2.bar_label(container)\nplt.title('Fractures by vertebrae')\nplt.ylim([0,2800])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:32.423549Z","iopub.execute_input":"2022-08-18T22:10:32.4241Z","iopub.status.idle":"2022-08-18T22:10:32.847515Z","shell.execute_reply.started":"2022-08-18T22:10:32.424067Z","shell.execute_reply":"2022-08-18T22:10:32.84652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12,5))\nax = sns.countplot(x = train_df[['C1','C2','C3','C4','C5','C6','C7']].sum(axis=1))\nfor container in ax.containers:\n ax.bar_label(container)\nplt.title('Number of fractures by patient')\nplt.xlabel('Number of fractures')\nplt.ylim([0,1300])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:32.849113Z","iopub.execute_input":"2022-08-18T22:10:32.849566Z","iopub.status.idle":"2022-08-18T22:10:33.050358Z","shell.execute_reply.started":"2022-08-18T22:10:32.849537Z","shell.execute_reply":"2022-08-18T22:10:33.049266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Heatmap of correlations\nplt.figure(figsize=(6,5))\nsns.heatmap(train_df[['C1','C2','C3','C4','C5','C6','C7']].corr(), cmap='bwr', vmin=-1, vmax=1)\nplt.title('Correlations')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:33.051784Z","iopub.execute_input":"2022-08-18T22:10:33.052235Z","iopub.status.idle":"2022-08-18T22:10:33.266607Z","shell.execute_reply.started":"2022-08-18T22:10:33.052205Z","shell.execute_reply":"2022-08-18T22:10:33.266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

3.2 Study Id's

\n
\n\nThe cases in the dataset have **unique id's** like '1.2.826.0.1.3680043.6200'. It turns out only the **number after the last full stop** is important. ","metadata":{}},{"cell_type":"code","source":"# Example\ntrain_df['StudyInstanceUID'][0]","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:33.26949Z","iopub.execute_input":"2022-08-18T22:10:33.269778Z","iopub.status.idle":"2022-08-18T22:10:33.27537Z","shell.execute_reply.started":"2022-08-18T22:10:33.269739Z","shell.execute_reply":"2022-08-18T22:10:33.274498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find unique numbers in study id's\nfor i in range(7):\n print(train_df['StudyInstanceUID'].map(lambda x : x.split('.')[i]).unique())","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:33.276433Z","iopub.execute_input":"2022-08-18T22:10:33.276689Z","iopub.status.idle":"2022-08-18T22:10:33.299905Z","shell.execute_reply.started":"2022-08-18T22:10:33.276666Z","shell.execute_reply":"2022-08-18T22:10:33.298933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Train images\n\n
\n

4.1 What is DICOM?

\n
\n\nA **.dcm** file follows the **Digital Imaging and Communications in Medicine** (DICOM) format. It is the standard format used for storing **medical images** and **related metadata**. It dates back to 1983, although it has been revised many times. \n\nWe can use the [pydicom library](https://pydicom.github.io/) to open and explore these files.","metadata":{}},{"cell_type":"code","source":"ex_path = \"../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.10001/101.dcm\"\ndcm_example = pydicom.dcmread(ex_path)\ndcm_example","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:33.30112Z","iopub.execute_input":"2022-08-18T22:10:33.301717Z","iopub.status.idle":"2022-08-18T22:10:33.325335Z","shell.execute_reply.started":"2022-08-18T22:10:33.301683Z","shell.execute_reply":"2022-08-18T22:10:33.324426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The **image data** is stored in an array under **'Pixel Data'**. Everything else is **metadata**.\n* The **'Rows'** and **'Columns'** values tell us the **image size**.\n* The **'Pixel Spacing'** and **'Slice Thickness'** tell us the **pixel size** and **thickness**.\n* The **'Window Center'** and **'Window Width'** give information about the **brightness** and **contrast** of the image respectively.\n* The **'Rescale Intercept'** and **'Rescale Slope'** determine the range of pixel values. ([ref](https://stackoverflow.com/questions/10193971/rescale-slope-and-rescale-intercept)).\n* **'ImagePositionPatient'** tells us the x, y, and z coordinates of the top left corner of each image in mm\n* **InstanceNumber** is the slice number.","metadata":{}},{"cell_type":"markdown","source":"
\n

4.2 Exploring the images

\n
\n\nLet's look at the what the **image data** from the dcm files look like.","metadata":{}},{"cell_type":"code","source":"# Adapted from https://www.kaggle.com/code/andradaolteanu/rsna-fracture-detection-dicom-images-explore\nbase_path = \"../input/rsna-2022-cervical-spine-fracture-detection\"\npatient_id = '1.2.826.0.1.3680043.12281'\ndcm_paths = glob(f\"{base_path}/train_images/{patient_id}/*\")\ndef atoi(text):\n return int(text) if text.isdigit() else text\ndef natural_keys(text):\n return [atoi(c) for c in re.split(r'(\\d+)', text)]\ndcm_paths.sort(key=natural_keys)\n\n# Get images\nfiles = [pydicom.dcmread(path) for path in dcm_paths]\nimages = [apply_voi_lut(file.pixel_array, file) for file in files]\n\n# Plot images\nfig, axes = plt.subplots(nrows=3, ncols=6, figsize=(24,12))\nfig.suptitle(f'ID: {patient_id}', weight=\"bold\", size=20)\n\nstart = 110\nfor i in range(start,start+18):\n img = images[i]\n file = files[i]\n slice_no = i\n\n # Plot the image\n x = (i-start) // 6\n y = (i-start) % 6\n\n axes[x, y].imshow(img, cmap=\"bone\")\n axes[x, y].set_title(f\"Slice: {slice_no}\", fontsize=14, weight='bold')\n axes[x, y].axis('off')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:33.326587Z","iopub.execute_input":"2022-08-18T22:10:33.327429Z","iopub.status.idle":"2022-08-18T22:10:39.525493Z","shell.execute_reply.started":"2022-08-18T22:10:33.327395Z","shell.execute_reply":"2022-08-18T22:10:39.524678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is great, but at the moment we don't know which **vertebrae** is being shown in each image. One way to work this out is by using the **segmentations**.","metadata":{}},{"cell_type":"markdown","source":"# 5. Segmentations\n\n
\n

5.1 What is NIfTI?

\n
\n\nA **.nii** file follows the **Neuroimaging Informatics Technology Initiative** (NIfTI) format. Compared to the DICOM, NIfTI is **simpler** and **easier** to support. \n\nTo open .nii files we can use the [nibabel library](https://nipy.org/nibabel/gettingstarted.html).","metadata":{}},{"cell_type":"code","source":"ex_path2 = f\"../input/rsna-2022-cervical-spine-fracture-detection/segmentations/{patient_id}.nii\"\nnii_example = nib.load(ex_path2)\n\n# Convert to numpy array\nseg = nii_example.get_fdata()\nseg.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:39.526486Z","iopub.execute_input":"2022-08-18T22:10:39.527297Z","iopub.status.idle":"2022-08-18T22:10:40.911234Z","shell.execute_reply.started":"2022-08-18T22:10:39.527266Z","shell.execute_reply":"2022-08-18T22:10:40.910509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each nifti file contains segmentations for **all slices** in a scan. However, we need to be careful about the **orientation** of the segmentations.\n\n> Please be aware that the NIFTI files consist of segmentation in the sagittal plane, while the DICOM files are in the axial plane.\n\nThe correct way to deal with this is explained in this [discussion post](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/340612). We need to use the following line:","metadata":{}},{"cell_type":"code","source":"# Align orientation with images\nseg = seg[:, ::-1, ::-1].transpose(2, 1, 0)\nseg.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:40.912357Z","iopub.execute_input":"2022-08-18T22:10:40.913166Z","iopub.status.idle":"2022-08-18T22:10:40.918857Z","shell.execute_reply.started":"2022-08-18T22:10:40.913136Z","shell.execute_reply":"2022-08-18T22:10:40.918024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n
\n

5.2 Exploring the masks

\n
","metadata":{"execution":{"iopub.status.busy":"2022-08-12T13:03:43.13158Z","iopub.execute_input":"2022-08-12T13:03:43.131927Z","iopub.status.idle":"2022-08-12T13:03:43.137208Z","shell.execute_reply.started":"2022-08-12T13:03:43.131896Z","shell.execute_reply":"2022-08-12T13:03:43.136019Z"}}},{"cell_type":"code","source":"# Plot images\nfig, axes = plt.subplots(nrows=3, ncols=6, figsize=(24,12))\nfig.suptitle(f'ID: {patient_id}', weight=\"bold\", size=20)\n\nfor i in range(start,start+18):\n mask = seg[i]\n slice_no = i\n\n # Plot the image\n x = (i-start) // 6\n y = (i-start) % 6\n\n axes[x, y].imshow(mask, cmap='inferno')\n axes[x, y].set_title(f\"Slice: {slice_no}\", fontsize=14, weight='bold')\n axes[x, y].axis('off')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:40.919871Z","iopub.execute_input":"2022-08-18T22:10:40.920143Z","iopub.status.idle":"2022-08-18T22:10:42.486349Z","shell.execute_reply.started":"2022-08-18T22:10:40.920118Z","shell.execute_reply":"2022-08-18T22:10:42.485352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compare these with the previous images. These masks give us the **location** of the vertebrae, which is very helpful because we know the fractures can only occur in these regions.\n\nThey also tells us which **vertebrae** are in the images. By looking at the unique values in each slice, we find a 0 for the background and another number like 2 for vertebrae C2. ","metadata":{}},{"cell_type":"code","source":"np.unique(seg[116])","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:42.487446Z","iopub.execute_input":"2022-08-18T22:10:42.488329Z","iopub.status.idle":"2022-08-18T22:10:42.499295Z","shell.execute_reply.started":"2022-08-18T22:10:42.488292Z","shell.execute_reply":"2022-08-18T22:10:42.498375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Unfortunately, we don't have segmentations for the whole train set. Checkout my other [notebook](https://www.kaggle.com/code/samuelcortinhas/extracting-vertebrae-c1-c7), where I exlore how to **extract the vertebrae** for the **whole train set**.","metadata":{}},{"cell_type":"code","source":"# Number of cases with masks\nseg_paths = glob(f\"{base_path}/segmentations/*\")\nprint(f'Number of cases with segmentations: {len(seg_paths)}, ({np.round(100*len(seg_paths)/len(train_df),1)}%)')","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:42.500387Z","iopub.execute_input":"2022-08-18T22:10:42.500719Z","iopub.status.idle":"2022-08-18T22:10:42.529384Z","shell.execute_reply.started":"2022-08-18T22:10:42.500694Z","shell.execute_reply":"2022-08-18T22:10:42.528146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"An idea to keep in mind is that maybe we could train a **segmentation model** like UNet to predict the segmentation masks for the rest of the train and all of the test images. ","metadata":{}},{"cell_type":"markdown","source":"# 6. Extract metadata\n\n
\n

6.1 Save metadata

\n
\n\nThis section uses code from Andrada Olteanu.","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/code/andradaolteanu/rsna-fracture-detection-dicom-images-explore\ndef get_observation_data(path):\n '''\n Get information from the .dcm files\n '''\n\n dataset = pydicom.read_file(path)\n \n # Dictionary to store the information from the image\n observation_data = {\n \"Rows\" : dataset.get(\"Rows\"),\n \"Columns\" : dataset.get(\"Columns\"),\n \"SOPInstanceUID\" : dataset.get(\"SOPInstanceUID\"),\n \"ContentDate\" : dataset.get(\"ContentDate\"),\n \"SliceThickness\" : dataset.get(\"SliceThickness\"),\n \"InstanceNumber\" : dataset.get(\"InstanceNumber\"),\n \"ImagePositionPatient\" : dataset.get(\"ImagePositionPatient\"),\n \"ImageOrientationPatient\" : dataset.get(\"ImageOrientationPatient\"),\n }\n\n # String columns\n str_columns = [\"SOPInstanceUID\", \"ContentDate\", \n \"SliceThickness\", \"InstanceNumber\"]\n for k in str_columns:\n observation_data[k] = str(dataset.get(k)) if k in dataset else None\n\n return observation_data\n\ndef get_metadata():\n '''\n Retrieves the desired metadata from the .dcm files and saves it into dataframe.\n '''\n \n exceptions = 0\n dicts = []\n\n for k in tqdm(range(len(train_df))):\n if (k % 100)==0:\n print(f'Iteration: {k}')\n \n dt = train_df.iloc[k, :]\n\n # Get all .dcm paths for this Instance\n dcm_paths = glob(f\"{base_path}/train_images/{dt.StudyInstanceUID}/*\")\n\n for path in dcm_paths:\n try:\n # Get datasets\n dataset = get_observation_data(path)\n dicts.append(dataset)\n except Exception as e:\n exceptions += 1\n continue\n\n # Convert into df\n meta_train_data = pd.DataFrame(data=dicts, columns=md_example.keys())\n \n # Export information\n meta_train_data.to_csv(\"meta_train.csv\", index=False)\n \n print(f\"Metadata created. Number of total fails: {exceptions}.\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:42.530807Z","iopub.execute_input":"2022-08-18T22:10:42.531142Z","iopub.status.idle":"2022-08-18T22:10:42.541618Z","shell.execute_reply.started":"2022-08-18T22:10:42.531111Z","shell.execute_reply":"2022-08-18T22:10:42.540806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Example\nmd_example = get_observation_data(ex_path)\npprint(md_example)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:42.542796Z","iopub.execute_input":"2022-08-18T22:10:42.543227Z","iopub.status.idle":"2022-08-18T22:10:42.563809Z","shell.execute_reply.started":"2022-08-18T22:10:42.543197Z","shell.execute_reply":"2022-08-18T22:10:42.563025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create and save the metadata (~ 2 hours)\n#get_metadata()","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:42.564953Z","iopub.execute_input":"2022-08-18T22:10:42.565555Z","iopub.status.idle":"2022-08-18T22:10:42.578105Z","shell.execute_reply.started":"2022-08-18T22:10:42.565527Z","shell.execute_reply":"2022-08-18T22:10:42.577279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can find the dataset containing the metadata [here](https://www.kaggle.com/datasets/samuelcortinhas/rsna-2022-spine-fracture-detection-metadata). (Andrada also has a version containing only the positive patients)","metadata":{}},{"cell_type":"markdown","source":"
\n

6.2 Explore metadata

\n
","metadata":{}},{"cell_type":"code","source":"# Read in saved metadata\nmeta_train = pd.read_csv(\"../input/rsna-2022-spine-fracture-detection-metadata/meta_train.csv\")\nmeta_train[\"StudyInstanceUID\"] = meta_train[\"SOPInstanceUID\"].apply(lambda x: \".\".join(x.split(\".\")[:-2]))\nprint('meta_train shape:', meta_train.shape)\nmeta_train.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:42.579136Z","iopub.execute_input":"2022-08-18T22:10:42.579524Z","iopub.status.idle":"2022-08-18T22:10:45.279922Z","shell.execute_reply.started":"2022-08-18T22:10:42.579501Z","shell.execute_reply":"2022-08-18T22:10:45.278904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Image size**","metadata":{}},{"cell_type":"code","source":"# Image size\nmeta_train[\"ImageSize\"] = meta_train[\"Rows\"].astype(str) + \" x \" + meta_train[\"Columns\"].astype(str)\n\n# Plot image sizes\nplt.figure(figsize=(18, 5))\nax = sns.countplot(data=meta_train, x=\"ImageSize\")\nfor container in ax.containers:\n ax.bar_label(container)\nplt.ylim([0,800000])\nplt.title('Image sizes in train images', fontsize=25, y=1.02)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:45.281055Z","iopub.execute_input":"2022-08-18T22:10:45.281336Z","iopub.status.idle":"2022-08-18T22:10:46.532306Z","shell.execute_reply.started":"2022-08-18T22:10:45.28131Z","shell.execute_reply":"2022-08-18T22:10:46.531557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Almost all images have size **512x512**. We should therefore **resize** the other images to size 512x512. ","metadata":{}},{"cell_type":"markdown","source":"**Content Date**","metadata":{}},{"cell_type":"code","source":"# Unique values\nmeta_train['ContentDate'].unique()","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:46.533256Z","iopub.execute_input":"2022-08-18T22:10:46.53397Z","iopub.status.idle":"2022-08-18T22:10:46.541784Z","shell.execute_reply.started":"2022-08-18T22:10:46.533945Z","shell.execute_reply":"2022-08-18T22:10:46.54089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can drop this feature as it is constant.","metadata":{}},{"cell_type":"code","source":"meta_train.drop('ContentDate', axis=1, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:46.542913Z","iopub.execute_input":"2022-08-18T22:10:46.543162Z","iopub.status.idle":"2022-08-18T22:10:46.694019Z","shell.execute_reply.started":"2022-08-18T22:10:46.543139Z","shell.execute_reply":"2022-08-18T22:10:46.692984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Slice Thickness**","metadata":{}},{"cell_type":"code","source":"meta_train[\"SliceThickness\"].unique()","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:46.695795Z","iopub.execute_input":"2022-08-18T22:10:46.696192Z","iopub.status.idle":"2022-08-18T22:10:46.711334Z","shell.execute_reply.started":"2022-08-18T22:10:46.696161Z","shell.execute_reply":"2022-08-18T22:10:46.710545Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot slice thickness\nplt.figure(figsize=(18, 5))\nax = sns.countplot(data=meta_train, x=\"SliceThickness\")\nfor container in ax.containers:\n ax.bar_label(container)\nax.set_xticklabels(['0.488','0.500','0.600','0.60...2','0.625','0.664','0.670','0.75','0.800','0.900','1.000'])\nplt.ylim([0,420000])\nplt.title('Slice thickness distribution', fontsize=25, y=1.02)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:46.712394Z","iopub.execute_input":"2022-08-18T22:10:46.713601Z","iopub.status.idle":"2022-08-18T22:10:46.975885Z","shell.execute_reply.started":"2022-08-18T22:10:46.71355Z","shell.execute_reply":"2022-08-18T22:10:46.974941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Number of slices per scan**","metadata":{}},{"cell_type":"code","source":"# Slice counts\nslice_counts = meta_train[\"StudyInstanceUID\"].value_counts().reset_index()\nslice_counts.columns = [\"StudyInstanceUID\", \"count\"]\n\n# Distribution of slices counts\nplt.figure(figsize=(18, 5))\nsns.histplot(data=slice_counts, x=\"count\", kde=True, bins=40)\nplt.title(\"Number of slices by scan\", size=25, y=1.02)\nplt.xlabel(\"Number of Slices\", size = 18)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:46.981509Z","iopub.execute_input":"2022-08-18T22:10:46.981798Z","iopub.status.idle":"2022-08-18T22:10:47.271372Z","shell.execute_reply.started":"2022-08-18T22:10:46.981772Z","shell.execute_reply":"2022-08-18T22:10:47.270436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Image Position Patient**","metadata":{}},{"cell_type":"code","source":"# Extract x, y, z coordinates of position vector\nmeta_train['ImagePositionPatient_x'] = meta_train['ImagePositionPatient'].apply(lambda x: float(x.replace(',','').replace(']','').replace('[','').split()[0]))\nmeta_train['ImagePositionPatient_y'] = meta_train['ImagePositionPatient'].apply(lambda x: float(x.replace(',','').replace(']','').replace('[','').split()[1]))\nmeta_train['ImagePositionPatient_z'] = meta_train['ImagePositionPatient'].apply(lambda x: float(x.replace(',','').replace(']','').replace('[','').split()[2]))","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:47.272459Z","iopub.execute_input":"2022-08-18T22:10:47.272718Z","iopub.status.idle":"2022-08-18T22:10:49.584149Z","shell.execute_reply.started":"2022-08-18T22:10:47.272693Z","shell.execute_reply":"2022-08-18T22:10:49.582874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot position coordinates\nplt.figure(figsize=(20,5))\nplt.subplot(1,3,1)\nsns.histplot(meta_train['ImagePositionPatient_x'])\nplt.ylim([0,20000])\nplt.title('x-coordinate', fontsize=25, y=1.02)\n\nplt.subplot(1,3,2)\nsns.histplot(meta_train['ImagePositionPatient_y'], color='C1')\nplt.ylabel('')\nplt.ylim([0,20000])\nplt.title('y-coordinate', fontsize=25, y=1.02)\n\nplt.subplot(1,3,3)\nsns.histplot(meta_train['ImagePositionPatient_z'], color='C2')\nplt.ylabel('')\nplt.ylim([0,20000])\nplt.title('z-coordinate', fontsize=25, y=1.02)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:49.585289Z","iopub.execute_input":"2022-08-18T22:10:49.585537Z","iopub.status.idle":"2022-08-18T22:10:53.334404Z","shell.execute_reply.started":"2022-08-18T22:10:49.585513Z","shell.execute_reply":"2022-08-18T22:10:53.333487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can use the 'ImagePositionPatient_z' feature to infer the vertebrae present in each slice. It seems non-trivial though - I need to think about this.","metadata":{}},{"cell_type":"markdown","source":"**Image Orientation Patient**","metadata":{}},{"cell_type":"markdown","source":"This feature tells us coordinates for the position of the patient in the scanner. It could tell us if the images are slightly rotated for example. ","metadata":{}},{"cell_type":"code","source":"meta_train['ImageOrientationPatient'].nunique()","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:53.335333Z","iopub.execute_input":"2022-08-18T22:10:53.335635Z","iopub.status.idle":"2022-08-18T22:10:53.412808Z","shell.execute_reply.started":"2022-08-18T22:10:53.335609Z","shell.execute_reply":"2022-08-18T22:10:53.41196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_train['ImageOrientationPatient'].unique()[:5]","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:53.41405Z","iopub.execute_input":"2022-08-18T22:10:53.414337Z","iopub.status.idle":"2022-08-18T22:10:53.48938Z","shell.execute_reply.started":"2022-08-18T22:10:53.414311Z","shell.execute_reply":"2022-08-18T22:10:53.488248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Clean meta data**","metadata":{}},{"cell_type":"code","source":"# Clean metadata\nmeta_train_clean = meta_train.drop(['SOPInstanceUID','ImagePositionPatient','ImageOrientationPatient','ImageSize'], axis=1)\nmeta_train_clean.rename(columns={\"Rows\": \"ImageHeight\", \"Columns\": \"ImageWidth\",\"InstanceNumber\": \"Slice\"}, inplace=True)\nmeta_train_clean = meta_train_clean[['StudyInstanceUID','Slice','ImageHeight','ImageWidth','SliceThickness','ImagePositionPatient_x','ImagePositionPatient_y','ImagePositionPatient_z']]\nmeta_train_clean.sort_values(by=['StudyInstanceUID','Slice'], inplace=True)\nmeta_train_clean.reset_index(drop=True, inplace=True)\n\n# Export information\n#meta_train_clean.to_csv(\"meta_train_clean.csv\", index=False)\n\n# Preview first few columns\nmeta_train_clean.head(3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:53.49103Z","iopub.execute_input":"2022-08-18T22:10:53.491964Z","iopub.status.idle":"2022-08-18T22:10:53.783877Z","shell.execute_reply.started":"2022-08-18T22:10:53.491922Z","shell.execute_reply":"2022-08-18T22:10:53.782742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The **cleaned version** of this metadata is also in my dataset [here](https://www.kaggle.com/datasets/samuelcortinhas/rsna-2022-spine-fracture-detection-metadata).","metadata":{}},{"cell_type":"markdown","source":"# 7. Bounding boxes\n\n
\n

7.1 Box distributions

\n
\n\nWe are only given bounding boxes for a **subset** of the data. In particular, only **12%** of patients in the train set have any bounding box measurements.\n\nThis information is useful in telling us exactly where the fractures have occured. We could consider training an **object localisation** algorithm to provide bounding boxes for the whole train set. ","metadata":{}},{"cell_type":"code","source":"print(f'Patients with bounding box measurements: {train_bbox[\"StudyInstanceUID\"].nunique()} ({np.round(100*train_bbox[\"StudyInstanceUID\"].nunique()/len(train_df),1)} %)')","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:53.78597Z","iopub.execute_input":"2022-08-18T22:10:53.786339Z","iopub.status.idle":"2022-08-18T22:10:53.794099Z","shell.execute_reply.started":"2022-08-18T22:10:53.786307Z","shell.execute_reply":"2022-08-18T22:10:53.793186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The bounding boxes are only provided for patients with fractures (but not all of them).","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/code/leventelippenszky/rsna-eda-dicom-segmentations-bboxes-3d-plot\ntrain_df_bbox = train_df[train_df[\"StudyInstanceUID\"].isin(train_bbox[\"StudyInstanceUID\"])]\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 5))\nsns.countplot(x=\"patient_overall\", data=train_df_bbox, ax=ax1)\nax1.set_title(\"Fracture overall (patients with bboxes)\")\n\ntrain_df_bbox_melt = pd.melt(train_df_bbox, id_vars=[\"StudyInstanceUID\", \"patient_overall\"], var_name=\"cervical_vertebrae\", value_name=\"fracture\")\nsns.countplot(x=\"cervical_vertebrae\", hue=\"fracture\", data=train_df_bbox_melt, ax=ax2)\nax2.set_title(\"Fracture of cervical vertebrae (patients with bboxes)\")\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:53.795481Z","iopub.execute_input":"2022-08-18T22:10:53.795814Z","iopub.status.idle":"2022-08-18T22:10:54.098428Z","shell.execute_reply.started":"2022-08-18T22:10:53.795783Z","shell.execute_reply":"2022-08-18T22:10:54.097313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(16,5))\nsns.histplot(train_bbox[\"StudyInstanceUID\"].value_counts().values, kde=True, bins=40)\nplt.title('Number of slices with bounding boxes per patient')\nplt.xlabel('Number of bboxes')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:54.100873Z","iopub.execute_input":"2022-08-18T22:10:54.101924Z","iopub.status.idle":"2022-08-18T22:10:54.376868Z","shell.execute_reply.started":"2022-08-18T22:10:54.101881Z","shell.execute_reply":"2022-08-18T22:10:54.376203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(16,6))\nplt.subplot(1,2,1)\nsns.scatterplot(data=train_bbox, x='x', y='y')\nplt.title('Anchor coordinates')\n\nplt.subplot(1,2,2)\nsns.scatterplot(data=train_bbox, x='width', y='height')\nplt.title('Width and heights')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-18T22:10:54.377766Z","iopub.execute_input":"2022-08-18T22:10:54.378534Z","iopub.status.idle":"2022-08-18T22:10:54.763106Z","shell.execute_reply.started":"2022-08-18T22:10:54.378504Z","shell.execute_reply":"2022-08-18T22:10:54.761666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

7.2 Examples of fractures

\n
","metadata":{}},{"cell_type":"code","source":"def plot_fracture(slice_num,bbox_id,ax_id1,ax_id2):\n file = pydicom.dcmread(f\"{base_path}/train_images/{bbox_id}/{slice_num}.dcm\")\n img = apply_voi_lut(file.pixel_array, file)\n info = train_bbox[(train_bbox['StudyInstanceUID']==bbox_id)&(train_bbox['slice_number']==slice_num)]\n rect = patches.Rectangle((float(info.x), float(info.y)), float(info.width), float(info.height), linewidth=3, edgecolor='r', facecolor='none')\n\n axes[ax_id1,ax_id2].imshow(img, cmap=\"bone\")\n axes[ax_id1,ax_id2].add_patch(rect)\n axes[ax_id1,ax_id2].set_title(f\"ID:{bbox_id}, Slice: {slice_num}\", fontsize=20, weight='bold',y=1.02)\n axes[ax_id1,ax_id2].axis('off')\n\nfig, axes = plt.subplots(nrows=2, ncols=2, figsize=(24,24))\nplot_fracture(119,'1.2.826.0.1.3680043.25651',0,0)\nplot_fracture(156,'1.2.826.0.1.3680043.23817',0,1)\nplot_fracture(325,'1.2.826.0.1.3680043.12031',1,0)\nplot_fracture(151,'1.2.826.0.1.3680043.11899',1,1)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:54.76431Z","iopub.execute_input":"2022-08-18T22:10:54.764579Z","iopub.status.idle":"2022-08-18T22:10:55.863844Z","shell.execute_reply.started":"2022-08-18T22:10:54.764554Z","shell.execute_reply":"2022-08-18T22:10:55.862936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 8. Baseline solution\n\n
\n

8.1 Submission

\n
\n\nWe use the **means** but **scale them up** a bit to take into account the competition metric puts **more weight on positive cases**. (See version 14 of this notebook for a derivation of the scaling function.)","metadata":{}},{"cell_type":"code","source":"# Put more weight on positive predictions\ndef scale_up(q):\n return 2*q/(1+q)\n\npreds = train_df.mean(numeric_only=True).map(scale_up).to_dict()\nss['fractured'] = test_df['prediction_type'].map(preds)\nss.to_csv('submission.csv', index=False)\nss.head(3)","metadata":{"execution":{"iopub.status.busy":"2022-08-18T22:10:55.864916Z","iopub.execute_input":"2022-08-18T22:10:55.865495Z","iopub.status.idle":"2022-08-18T22:10:55.881222Z","shell.execute_reply.started":"2022-08-18T22:10:55.865465Z","shell.execute_reply":"2022-08-18T22:10:55.880323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# -1. References\n\n* [RSNA Fracture Detection: DICOM & Images Explore](https://www.kaggle.com/code/andradaolteanu/rsna-fracture-detection-dicom-images-explore) by [Andrada Olteanu](https://www.kaggle.com/andradaolteanu).\n* [Fork of great EDA with fix for slice count dist](https://www.kaggle.com/code/kretes/fork-of-great-eda-with-fix-for-slice-count-dist) by [Tomasz Bartczak](https://www.kaggle.com/kretes).\n* [What are .DCM and .NII files and how to read them](https://www.kaggle.com/code/datark1/what-are-dcm-and-nii-files-and-how-to-read-them) by [Robert Kwiatkowski](https://www.kaggle.com/datark1).\n* [RSNA CSF - Cervical Spine Fracture EDA](https://www.kaggle.com/code/allunia/rsna-csf-cervical-spine-fracture-eda) by [Laura Fink](https://www.kaggle.com/allunia).\n* [Explaining Data and Submission in detail](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/340612) by [Harshit Sheoran](https://www.kaggle.com/harshitsheoran).\n* [Cervical Spine Fracture Detection Quick EDA](https://www.kaggle.com/code/ipythonx/cervical-spine-fracture-detection-quick-eda) by [M.Innat](https://www.kaggle.com/ipythonx).\n* [RSNA EDA: DICOM, segmentations, bboxes, 3D plot](https://www.kaggle.com/code/leventelippenszky/rsna-eda-dicom-segmentations-bboxes-3d-plot/notebook) by [_lev_lipinski](https://www.kaggle.com/leventelippenszky).","metadata":{}}]}