{"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":"\n\n

- Exploratory Phase -

\n\n> **Goal**: Detect and localize cervical spine fractures within CT scans.\n\n*Competition timeline:*\n \n- Start date - 28th July 2022\n- Finish date - 27th October 2022\n\nAll the details about the competition can be found [here](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection).\n\n### References\n\n- [🦴 RSNA Fracture Detection: DICOM & Images Explore](https://www.kaggle.com/code/andradaolteanu/rsna-fracture-detection-dicom-images-explore)\n- [🦴 RSNA Fracture Detection - in-depth EDA](https://www.kaggle.com/code/samuelcortinhas/rsna-fracture-detection-in-depth-eda)","metadata":{"papermill":{"duration":0.018901,"end_time":"2022-10-19T01:39:07.157041","exception":false,"start_time":"2022-10-19T01:39:07.13814","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"markdown","source":"# 1. Introduction\n\n
\n

1.1 Background

\n
\n\n### Cervical Fractures (Broken Neck)\n\nThere are **7 bones that make up the cervical vertebrae** (neck). They support the head and connect it to the shoulders and body. A fracture, or break, in one of the cervical vertebrae is commonly called a broken neck ([source](https://orthoinfo.aaos.org/en/diseases--conditions/cervical-fracture-broken-neck/#:~:text=A%20fracture%2C%20or%20break%2C%20in,Athletes%20are%20also%20at%20risk.)).\n\n
","metadata":{"papermill":{"duration":0.016761,"end_time":"2022-10-19T01:39:07.192662","exception":false,"start_time":"2022-10-19T01:39:07.175901","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"markdown","source":"
\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, & \\text{if vertebrae negative} \\\\\n2, & \\text{if vertebrae positive} \\\\\n7, & \\text{if patient negative} \\\\\n14, & \\text{if patient positive}\n\\end{cases}\n$$\n\n
\n
\n\n
\n
","metadata":{"papermill":{"duration":0.016355,"end_time":"2022-10-19T01:39:07.226296","exception":false,"start_time":"2022-10-19T01:39:07.209941","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"markdown","source":"Notice how **more weight** is put on **positive cases** and the **most weight** on the **overall probability** of any fractures.\n\n---\n\n**Example:**\n\nThe table below illustrates **how the evaluation metric is derived** for a particular case id. This is then **averaged** over all case ids.\n\n| Category | Label | Prediction | Loss |\n|:----------:|:-------:|:------------:|:--------------------------------:|\n| C1 | 1 | 0.8 | $$-2\\log(0.8) \\approx 0.446$$ |\n| C2 | 0 | 0.1 | $$-\\log(1-0.1) \\approx 0.105$$ |\n| C3 | 0 | 0.2 | $$-\\log(1-0.2) \\approx 0.223$$ |\n| C4 | 0 | 0 | $$-\\log(1-0) \\approx 0$$ |\n| C5 | 1 | 0.6 | $$-2\\log(0.6) \\approx 1.022$$ |\n| C6 | 0 | 0.3 | $$-\\log(1-0.3) \\approx 0.357$$ |\n| C7 | 0 | 0.1 | $$-\\log(1-0.1) \\approx 0.105$$ |\n| Overall | 1 | 0.9 | $$-14\\log(0.9) \\approx 1.475$$ |\n| | | | $$\\text{Total} = 3.734$$ |\n\nNote: it is conventiontional 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 \\approx 0.162$. ","metadata":{"papermill":{"duration":0.016807,"end_time":"2022-10-19T01:39:07.260056","exception":false,"start_time":"2022-10-19T01:39:07.243249","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"markdown","source":"
\n

1.3 Libraries

\n
","metadata":{"papermill":{"duration":0.016907,"end_time":"2022-10-19T01:39:07.293671","exception":false,"start_time":"2022-10-19T01:39:07.276764","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[],"_kg_hide-input":true,"_kg_hide-output":true}},{"cell_type":"code","source":"!pip install -qU 'python-gdcm' pydicom pylibjpeg 'opencv-python-headless'","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:07.329008Z","iopub.status.busy":"2022-10-19T01:39:07.328395Z","iopub.status.idle":"2022-10-19T01:39:24.830101Z","shell.execute_reply":"2022-10-19T01:39:24.828927Z"},"papermill":{"duration":17.522708,"end_time":"2022-10-19T01:39:24.832927","exception":false,"start_time":"2022-10-19T01:39:07.310219","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[],"vscode":{"languageId":"shellscript"}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import re\nimport cv2\nimport random\nfrom glob import glob\nfrom pprint import pprint\nimport warnings\nimport itertools\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.patches as patches\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import ListedColormap\nfrom IPython.display import display_html\n\nplt.rcParams.update({'font.size': 16})\n\n# .dcm handling\nimport pydicom\nimport nibabel as nib\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\n# Environment check\nwarnings.filterwarnings('ignore')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:24.869447Z","iopub.status.busy":"2022-10-19T01:39:24.868938Z","iopub.status.idle":"2022-10-19T01:39:26.297238Z","shell.execute_reply":"2022-10-19T01:39:26.295929Z"},"papermill":{"duration":1.450037,"end_time":"2022-10-19T01:39:26.300105","exception":false,"start_time":"2022-10-19T01:39:24.850068","status":"completed"},"pycharm":{"is_executing":true,"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Custom colors\nclass clr:\n S = '\\033[1m' + '\\033[94m'\n E = '\\033[0m'\n\n\nmy_colors = [\n '#5EAFD9',\n '#449DD1',\n '#3977BB',\n '#2D51A5',\n '#5C4C8F',\n '#8B4679',\n '#C53D4C',\n '#E23836',\n '#FF4633',\n '#FF5746',\n]\nCMAP1 = ListedColormap(my_colors)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"collapsed":false,"execution":{"iopub.execute_input":"2022-10-19T01:39:26.336801Z","iopub.status.busy":"2022-10-19T01:39:26.336343Z","iopub.status.idle":"2022-10-19T01:39:26.342627Z","shell.execute_reply":"2022-10-19T01:39:26.341435Z"},"jupyter":{"outputs_hidden":false},"papermill":{"duration":0.027361,"end_time":"2022-10-19T01:39:26.345162","exception":false,"start_time":"2022-10-19T01:39:26.317801","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

1.4 Helper Functions

\n
","metadata":{"papermill":{"duration":0.016565,"end_time":"2022-10-19T01:39:26.378828","exception":false,"start_time":"2022-10-19T01:39:26.362263","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[],"_kg_hide-input":true,"_kg_hide-output":true}},{"cell_type":"code","source":"def show_values_on_bars(axs, h_v='v', space=0.4):\n \"\"\"Plots the value at the end of the seaborn bar plot.\n axs: the ax of the plot\n h_v: weather or not the bar plot is vertical/ horizontal\"\"\"\n\n def _show_on_single_plot(ax):\n if h_v == 'v':\n for p in ax.patches:\n _x = p.get_x() + p.get_width() / 2\n _y = p.get_y() + p.get_height()\n value = int(p.get_height())\n ax.text(_x, _y, format(value, ','), ha='center')\n elif h_v == 'h':\n for p in ax.patches:\n _x = p.get_x() + p.get_width() + float(space)\n _y = p.get_y() + p.get_height()\n value = int(p.get_width())\n ax.text(_x, _y, format(value, ','), ha='left')\n\n if isinstance(axs, np.ndarray):\n for idx, ax in np.ndenumerate(axs):\n _show_on_single_plot(ax)\n else:\n _show_on_single_plot(axs)\n\n\ndef atoi(text):\n return int(text) if text.isdigit() else text\n\n\ndef natural_keys(text):\n \"\"\"\n alist.sort(key=natural_keys) sorts in human order\n https://nedbatchelder.com/blog/200712/human_sorting.html\n (See Toothy's implementation in the comments)\n \"\"\"\n return [atoi(c) for c in re.split(r'(\\d+)', text)]","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:26.414289Z","iopub.status.busy":"2022-10-19T01:39:26.41389Z","iopub.status.idle":"2022-10-19T01:39:26.425913Z","shell.execute_reply":"2022-10-19T01:39:26.424793Z"},"papermill":{"duration":0.032719,"end_time":"2022-10-19T01:39:26.428383","exception":false,"start_time":"2022-10-19T01:39:26.395664","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Data\n\n
\n

2.1 Data frames

\n
\n\n- **Train Dataframe**:\n - `patient_overall` - whether or not there is at least 1 vertebrae that is fractured\n - `C1-C7` - which of the vertebraes (if any) are fractured\n - 2019 *unique* IDs in total\n- **Train BBox Dataframe**:\n - contains bbox coordinates with the exact location of the fractures\n- **Train Images folder**:\n - 2019 folders containing the CT scans of the patient in `.dcm` format","metadata":{"papermill":{"duration":0.0174,"end_time":"2022-10-19T01:39:26.462681","exception":false,"start_time":"2022-10-19T01:39:26.445281","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"BASE = '../input/rsna-2022-cervical-spine-fracture-detection'","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:26.498542Z","iopub.status.busy":"2022-10-19T01:39:26.498138Z","iopub.status.idle":"2022-10-19T01:39:26.503264Z","shell.execute_reply":"2022-10-19T01:39:26.502014Z"},"papermill":{"duration":0.025478,"end_time":"2022-10-19T01:39:26.505532","exception":false,"start_time":"2022-10-19T01:39:26.480054","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_data():\n \"\"\"Reads in all .csv files.\"\"\"\n\n train = pd.read_csv(\n '../input/rsna-2022-cervical-spine-fracture-detection/train.csv'\n )\n train_bbox = pd.read_csv(\n '../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv'\n )\n test = pd.read_csv('../input/rsna-2022-cervical-spine-fracture-detection/test.csv')\n ss = pd.read_csv(\n '../input/rsna-2022-cervical-spine-fracture-detection/sample_submission.csv'\n )\n\n return train, train_bbox, test, ss\n\n\ndef get_csv_info(csv, name='Default'):\n \"\"\"Prints main information for the specified .csv file.\"\"\"\n\n print(f'{clr.S}=== {name} ==={clr.E}')\n print(f'{clr.S}Shape:{clr.E}', csv.shape)\n print(\n f'{clr.S}Missing Values:{clr.E}',\n csv.isna().sum().sum(),\n 'total missing datapoints.',\n )\n print(f'{clr.S}Columns:{clr.E}', list(csv.columns), '\\n')\n\n display_html(csv.head())\n print('\\n')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:26.540701Z","iopub.status.busy":"2022-10-19T01:39:26.54031Z","iopub.status.idle":"2022-10-19T01:39:26.549656Z","shell.execute_reply":"2022-10-19T01:39:26.54875Z"},"papermill":{"duration":0.029635,"end_time":"2022-10-19T01:39:26.551823","exception":false,"start_time":"2022-10-19T01:39:26.522188","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read in the data\ntrain, train_bbox, test, ss = read_data()\n\n# Print useful information on it\nfor csv, name in zip(\n [train, train_bbox, test, ss], ['Train', 'Train Bbox', 'Test', 'Sample Submission']\n):\n get_csv_info(csv, name)","metadata":{"_kg_hide-input":true,"_kg_hide-output":false,"execution":{"iopub.execute_input":"2022-10-19T01:39:26.587484Z","iopub.status.busy":"2022-10-19T01:39:26.587089Z","iopub.status.idle":"2022-10-19T01:39:26.683244Z","shell.execute_reply":"2022-10-19T01:39:26.682075Z"},"papermill":{"duration":0.118691,"end_time":"2022-10-19T01:39:26.687489","exception":false,"start_time":"2022-10-19T01:39:26.568798","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Exploratory Data Analysis","metadata":{"papermill":{"duration":0.017963,"end_time":"2022-10-19T01:39:26.725147","exception":false,"start_time":"2022-10-19T01:39:26.707184","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"markdown","source":"
\n

3.1 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.\n\nBelow we can see that the **first 6 strings** are **identical** for all cases. This means that only the **last string** is the **unique identifier** for each case.","metadata":{"papermill":{"duration":0.0179,"end_time":"2022-10-19T01:39:26.761175","exception":false,"start_time":"2022-10-19T01:39:26.743275","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Find unique numbers in study id's\nfor i in range(7):\n print(train['StudyInstanceUID'].map(lambda x: x.split('.')[i]).unique())","metadata":{"execution":{"iopub.execute_input":"2022-10-19T01:39:26.852806Z","iopub.status.busy":"2022-10-19T01:39:26.851789Z","iopub.status.idle":"2022-10-19T01:39:26.874896Z","shell.execute_reply":"2022-10-19T01:39:26.873252Z"},"papermill":{"duration":0.045542,"end_time":"2022-10-19T01:39:26.877652","exception":false,"start_time":"2022-10-19T01:39:26.83211","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[],"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

3.2 Fractures Analysis

\n
\n\n**Main Takeaways**:\n\n- *classes are balanced* - there is a ~50%/50% ratio between patients with and without fractures\n- *C7* - the bone that is most likely to be fractured\n- *C3* - the bone that is least likely to be fractured","metadata":{"papermill":{"duration":0.017615,"end_time":"2022-10-19T01:39:26.913364","exception":false,"start_time":"2022-10-19T01:39:26.895749","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"dt = pd.melt(\n train,\n id_vars=['StudyInstanceUID', 'patient_overall'],\n var_name='Vertebra',\n value_name='Flag',\n)\n\n# Plot\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(24, 12))\nfig.suptitle('Fractures - Main Analysis', weight='bold', size=25)\n\nsns.countplot(\n data=train, x='patient_overall', ax=ax1, palette=[my_colors[0], my_colors[4]]\n)\nshow_values_on_bars(ax1, h_v='v', space=0.4)\nax1.set_title('Patient Overall Fracture Flag [frequency]', weight='bold', size=19)\nax1.set_xlabel('Fracture Flag', size=18, weight='bold')\nax1.set_ylabel('')\nax1.set_yticks([])\n\nhatches = itertools.cycle(['', '//'])\nfor bar in ax1.patches:\n hatch = next(hatches)\n bar.set_hatch(hatch)\n\nsns.countplot(\n data=dt, x='Vertebra', hue='Flag', ax=ax2, palette=[my_colors[0], my_colors[4]]\n)\nshow_values_on_bars(ax2, h_v='v', space=0.4)\nax2.set_title('Fracture Flag on Vertebra Type [frequency]', weight='bold', size=19)\nax2.set_xlabel('Vertebra', size=18, weight='bold')\nax2.set_ylabel('')\nax2.set_yticks([])\n\nfor i, bar in enumerate(ax2.patches):\n hatch = ''\n if i in [7, 8, 9, 10, 11, 12, 13]:\n hatch = '//'\n bar.set_hatch(hatch)\n\nsns.despine(right=True, top=True, left=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:26.95343Z","iopub.status.busy":"2022-10-19T01:39:26.953022Z","iopub.status.idle":"2022-10-19T01:39:27.427546Z","shell.execute_reply":"2022-10-19T01:39:27.426296Z"},"papermill":{"duration":0.496931,"end_time":"2022-10-19T01:39:27.430366","exception":false,"start_time":"2022-10-19T01:39:26.933435","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

3.3 Number of Injuries per Patient

\n
\n\n**Main Takeaways**:\n\n- *single-bone injuries* - of the 961 cases with at least one fracture, more than half have a single-bone injury.\n- *multiple bone injuries* - relatively rare occurrences involve damage to four or more bones at the same time.","metadata":{"papermill":{"duration":0.018777,"end_time":"2022-10-19T01:39:27.46889","exception":false,"start_time":"2022-10-19T01:39:27.450113","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"train['total_fractures'] = train.iloc[:, 2:].sum(axis=1)\n\n# Plot\nplt.figure(figsize=(24, 12))\naxs = sns.countplot(data=train, x='total_fractures', palette=my_colors)\nshow_values_on_bars(axs, h_v='v', space=0.4)\nplt.title('Number of Bones Fractured per Patient', weight='bold', size=19)\nplt.xlabel('Total Bones Fractured', size=18, weight='bold')\nplt.ylabel('')\nplt.yticks([])\n\n# Hatch\nfor i, bar in enumerate(axs.patches):\n hatch = ''\n if i == 0:\n hatch = '-'\n bar.set_hatch(hatch)\n\n# Arrow\nstyle = 'Simple, tail_width=2, head_width=14, head_length=16'\nkw = dict(arrowstyle=style, color=my_colors[2])\narrow = patches.FancyArrowPatch(\n (1.8, 900), (1, 660), connectionstyle='arc3,rad=.10', **kw\n)\nplt.gca().add_patch(arrow)\nplt.text(\n x=2,\n y=900,\n s=f'{round(623 / 961 * 100, 2)}% cases have',\n color=my_colors[2],\n size=17,\n weight='bold',\n)\nplt.text(\n x=2, y=860, s='only 1 bone fractured', color=my_colors[2], size=17, weight='bold'\n)\n\nplt.text(x=2, y=820, s='at the same time.', color=my_colors[2], size=17, weight='bold')\n\n# Line\nplt.axvline(x=3.5, linestyle='--', color=my_colors[5], lw=2)\nplt.text(\n x=3.6,\n y=250,\n s=f'{round(35 / 961 * 100, 2)}% cases have 4+ bones fractured at the same time.',\n color=my_colors[5],\n size=17,\n weight='bold',\n)\n\nsns.despine(right=True, top=True, left=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:27.509274Z","iopub.status.busy":"2022-10-19T01:39:27.508835Z","iopub.status.idle":"2022-10-19T01:39:27.842103Z","shell.execute_reply":"2022-10-19T01:39:27.840949Z"},"papermill":{"duration":0.356777,"end_time":"2022-10-19T01:39:27.844738","exception":false,"start_time":"2022-10-19T01:39:27.487961","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Image Data [.dcm]\n\n
\n

4.1 Understanding Images

\n
\n\nLet's start with a few samples to get a sense of how these images look. Following that, we may extract further details from the `.dcm` files and investigate general features (like image size, number of slices per patient etc.).","metadata":{"papermill":{"duration":0.020107,"end_time":"2022-10-19T01:39:27.885598","exception":false,"start_time":"2022-10-19T01:39:27.865491","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"def show_dcm_images(patient_id, rows=4, cols=8, random=0):\n \"\"\"Show .dcm images based on id.\"\"\"\n\n N = rows * cols\n\n fig, axes = plt.subplots(nrows=rows, ncols=cols, figsize=(24, 15))\n fig.suptitle(f'ID: {patient_id}', weight='bold', size=20)\n\n # Get .dcm paths\n dcm_paths = glob(f'{BASE}/train_images/{patient_id}/*')\n print(f'{clr.S}Number of TOTAL Slices:{clr.E}', len(dcm_paths))\n dcm_paths.sort(key=natural_keys)\n dcm_paths = dcm_paths[random: (random + N)]\n # Get corresponding datasets and images\n datasets = [pydicom.dcmread(path) for path in dcm_paths]\n images = [apply_voi_lut(dataset.pixel_array, dataset) for dataset in datasets]\n\n # Loop through the information\n for data, img, i in zip(datasets, images, range(N)):\n slice_no = data.SOPInstanceUID.split('.')[-1]\n\n # Plot the image\n x = i // cols\n y = i % cols\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,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:27.928465Z","iopub.status.busy":"2022-10-19T01:39:27.928059Z","iopub.status.idle":"2022-10-19T01:39:27.93929Z","shell.execute_reply":"2022-10-19T01:39:27.937982Z"},"papermill":{"duration":0.03561,"end_time":"2022-10-19T01:39:27.941557","exception":false,"start_time":"2022-10-19T01:39:27.905947","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Example 1","metadata":{"papermill":{"duration":0.020334,"end_time":"2022-10-19T01:39:27.982433","exception":false,"start_time":"2022-10-19T01:39:27.962099","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"show_dcm_images('1.2.826.0.1.3680043.10001', rows=4, cols=8)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:28.025321Z","iopub.status.busy":"2022-10-19T01:39:28.024622Z","iopub.status.idle":"2022-10-19T01:39:31.009304Z","shell.execute_reply":"2022-10-19T01:39:31.007877Z"},"papermill":{"duration":3.022042,"end_time":"2022-10-19T01:39:31.02488","exception":false,"start_time":"2022-10-19T01:39:28.002838","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Example 2","metadata":{"papermill":{"duration":0.044456,"end_time":"2022-10-19T01:39:31.114607","exception":false,"start_time":"2022-10-19T01:39:31.070151","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"show_dcm_images('1.2.826.0.1.3680043.14994', rows=4, cols=8, random=23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:31.206793Z","iopub.status.busy":"2022-10-19T01:39:31.20636Z","iopub.status.idle":"2022-10-19T01:39:34.151382Z","shell.execute_reply":"2022-10-19T01:39:34.150164Z"},"papermill":{"duration":3.008111,"end_time":"2022-10-19T01:39:34.166633","exception":false,"start_time":"2022-10-19T01:39:31.158522","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

4.2 DICOM Metadata

\n
\n\n### Retrieving Metadata for a single file\n\nFor the time being, the only information I'd like to receive is the following:\n\n- Rows -> the height of the CT scan/image\n- Columns -> the width of the CT scan/image\n- SOPInstanceUID -> Unique identifier containing the `StudyInstanceUID` + slice number\n- ContentDate -> the date the image pixel data creation started\n- SliceThickness -> gives the thickness of the imaged slice (*TODO: maybe pair with `Spacing Between Slices` - gives the distance between two adjacent slices*)\n- InstanceNumber -> slice number\n- ImagePositionPatient -> the x, y, and z coordinates of the upper left hand corner (center of the first voxel transmitted) of the image, in *mm*\n- ImageOrientationPatient -> the direction cosines of the first row and the first column with respect to the patient\n\n> **Note**: all attribute explanations [are here](https://dicom.innolitics.com/ciods/rt-dose/image-plane/00200037).","metadata":{"papermill":{"duration":0.062758,"end_time":"2022-10-19T01:39:34.295845","exception":false,"start_time":"2022-10-19T01:39:34.233087","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"def 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', '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","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:34.426063Z","iopub.status.busy":"2022-10-19T01:39:34.424959Z","iopub.status.idle":"2022-10-19T01:39:34.433111Z","shell.execute_reply":"2022-10-19T01:39:34.432236Z"},"papermill":{"duration":0.076259,"end_time":"2022-10-19T01:39:34.435258","exception":false,"start_time":"2022-10-19T01:39:34.358999","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# An example\npath = '../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.10001/109.dcm'\nexample = get_observation_data(path)\npprint(example)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:34.564467Z","iopub.status.busy":"2022-10-19T01:39:34.563851Z","iopub.status.idle":"2022-10-19T01:39:34.581881Z","shell.execute_reply":"2022-10-19T01:39:34.580606Z"},"papermill":{"duration":0.085932,"end_time":"2022-10-19T01:39:34.584187","exception":false,"start_time":"2022-10-19T01:39:34.498255","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

4.3 Exploring Metadata

\n
\n\nLet's start exploring the metadata.","metadata":{"papermill":{"duration":0.06495,"end_time":"2022-10-19T01:39:34.716359","exception":false,"start_time":"2022-10-19T01:39:34.651409","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"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(\n lambda x: '.'.join(x.split('.')[:-2])\n)\n\n# Information\nprint(f'{clr.S}Total number of .dcm files: (CT scans){clr.E}', len(meta_train))\nprint(\n f'{clr.S}[sanity check] Total Number of Study Instances{clr.E}',\n meta_train['StudyInstanceUID'].nunique()\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:34.848665Z","iopub.status.busy":"2022-10-19T01:39:34.848267Z","iopub.status.idle":"2022-10-19T01:39:38.211877Z","shell.execute_reply":"2022-10-19T01:39:38.210693Z"},"papermill":{"duration":3.432887,"end_time":"2022-10-19T01:39:38.214555","exception":false,"start_time":"2022-10-19T01:39:34.781668","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Distribution of .dcm files on Study Instance\n\n**Please note**:\n\n- The distribution is skewed to the right\n- 78% of all observations ranged from 200 to 400 slices\n- Only 6% of study cases contain less than 200 slices\n- There are 50 Study Instances with more than 600 slices","metadata":{"papermill":{"duration":0.064072,"end_time":"2022-10-19T01:39:38.344631","exception":false,"start_time":"2022-10-19T01:39:38.280559","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Get the data\ndata = meta_train['StudyInstanceUID'].value_counts().reset_index()\ndata.columns = ['StudyInstanceUID', 'count']","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:38.474622Z","iopub.status.busy":"2022-10-19T01:39:38.473636Z","iopub.status.idle":"2022-10-19T01:39:38.551849Z","shell.execute_reply":"2022-10-19T01:39:38.550955Z"},"papermill":{"duration":0.145397,"end_time":"2022-10-19T01:39:38.554131","exception":false,"start_time":"2022-10-19T01:39:38.408734","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot 1\nplt.figure(figsize=(24, 12))\nsns.distplot(\n data['count'],\n rug=True,\n bins=10,\n rug_kws={'color': my_colors[1]},\n kde_kws={'color': my_colors[8], 'lw': 5, 'alpha': 0.7},\n hist_kws={'histtype': 'step', 'linewidth': 3, 'alpha': 1, 'color': my_colors[1]},\n)\n\nplt.title('Distribution of .dcm files on Study Instance', weight='bold', size=25)\nplt.xlabel('Number of Slices', size=18, weight='bold')\nplt.ylabel('Frequency')\n\n# Arrow\nstyle = 'Simple, tail_width=1, head_width=12, head_length=14'\nkw = dict(arrowstyle=style, color='black')\narrow = patches.FancyArrowPatch(\n (600, 0.0033), (300, 0.0020), connectionstyle='arc3,rad=-.10', **kw\n)\nplt.gca().add_patch(arrow)\n\nplt.text(\n x=400,\n y=0.0035,\n s='78% of observations have 200 to 400 slices.',\n color='black',\n size=17,\n)\n\nsns.despine(right=True, top=True, left=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:38.68743Z","iopub.status.busy":"2022-10-19T01:39:38.687048Z","iopub.status.idle":"2022-10-19T01:39:39.086307Z","shell.execute_reply":"2022-10-19T01:39:39.085028Z"},"papermill":{"duration":0.470328,"end_time":"2022-10-19T01:39:39.089327","exception":false,"start_time":"2022-10-19T01:39:38.618999","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"random.seed(25)\n\n# Add a fictive y\n# this 'y' doesn't mean anything, it's just for\n# showcasing purposes\ndata['y'] = [random.randint(0, 100) for _ in range(len(data))]\nperc = round(data[data['count'] <= 10].shape[0] / len(data), 3) * 100\n\nplt.figure(figsize=(24, 12))\nsns.scatterplot(\n data=data,\n x='count',\n y='y',\n size='count',\n alpha=0.65,\n sizes=(500, 2000),\n hue='count',\n palette=CMAP1,\n)\n\nplt.title('Distribution of .dcm files on Study Instance', weight='bold', size=25)\nplt.xlabel('Number of Slices', size=18, weight='bold')\nplt.ylabel('')\nplt.yticks([])\n\nplt.axvline(x=200, linestyle='--', color='black', lw=4)\nplt.axvline(x=400, linestyle='--', color='black', lw=4)\nplt.text(x=225, y=50, s='~80% of data is here', color='black', size=20, weight='bold')\n\nplt.axvline(x=600, linestyle='--', color='black', lw=4)\nplt.text(\n x=610,\n y=7,\n s='only 50 instances with slices >=600',\n color='black',\n size=14,\n weight='bold',\n)\n\nplt.arrow(x=600, y=5, dx=200, dy=0, color='black', lw=4, head_width=2, head_length=8)\n\nplt.legend('', frameon=False)\nsns.despine(right=True, top=True, left=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:39.226463Z","iopub.status.busy":"2022-10-19T01:39:39.225554Z","iopub.status.idle":"2022-10-19T01:39:39.812933Z","shell.execute_reply":"2022-10-19T01:39:39.811772Z"},"papermill":{"duration":0.668577,"end_time":"2022-10-19T01:39:39.826459","exception":false,"start_time":"2022-10-19T01:39:39.157882","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Image sizes\n\n**Note**:\n\n- 99.6% of images have a fixed size of 512 by 512.\n- The rest of the slices (the other 0.3%) should be resized to 512x512 as well.","metadata":{"papermill":{"duration":0.083734,"end_time":"2022-10-19T01:39:39.99847","exception":false,"start_time":"2022-10-19T01:39:39.914736","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Create column with full image size (height x width)\nmeta_train['ImageSize'] = (\n meta_train['Rows'].astype(str) + ' x ' + meta_train['Columns'].astype(str)\n)\n\n# Plot\nplt.figure(figsize=(24, 10))\naxs = sns.countplot(data=meta_train, x='ImageSize', palette=my_colors[5:])\nshow_values_on_bars(axs, h_v='v', space=0.4)\nplt.title('Frequency of Image Sizes in .dcm data', weight='bold', size=19)\nplt.xlabel('Image Size (height x width)', size=18, weight='bold')\nplt.ylabel('')\nplt.yticks([])\n\n# Hatch\nfor i, bar in enumerate(axs.patches):\n hatch = ''\n if i == 0:\n hatch = '/'\n bar.set_hatch(hatch)\n\n# Arrow\nstyle = 'Simple, tail_width=2, head_width=14, head_length=16'\nkw = dict(arrowstyle=style, color=my_colors[4])\narrow = patches.FancyArrowPatch(\n (0.8, 300000), (0.2, 200000), connectionstyle='arc3,rad=.10', **kw\n)\nplt.gca().add_patch(arrow)\nplt.text(\n x=0.8,\n y=300000,\n s='99.6% of images have size 512x512',\n color=my_colors[4],\n size=17,\n weight='bold',\n)\n\nplt.text(\n x=1,\n y=70000,\n s='The other 0.32% should be resized as 512x512',\n color=my_colors[6],\n size=17,\n weight='bold',\n)\n\nplt.axhline(xmin=0.4, xmax=0.95, y=60000, linestyle='--', color=my_colors[6], lw=2)\n\nsns.despine(right=True, top=True, left=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:40.17145Z","iopub.status.busy":"2022-10-19T01:39:40.171023Z","iopub.status.idle":"2022-10-19T01:39:41.865845Z","shell.execute_reply":"2022-10-19T01:39:41.864659Z"},"papermill":{"duration":1.782543,"end_time":"2022-10-19T01:39:41.868437","exception":false,"start_time":"2022-10-19T01:39:40.085894","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5. Segmentations\n\n
\n

5.1 What is NIfTI?

\n
\n\nThe segmentation data is provided in the `.nii` files, which follow the **Neuroimaging Informatics Technology Initiative** (NIfTI) format.\n\nTo open .nii files we use the [nibabel library](https://nipy.org/nibabel/gettingstarted.html).\n\nWe immediately notice that the data consists of segmentation in the sagittal plane (in contrast to the axial plane of the CT scans).","metadata":{"papermill":{"duration":0.085512,"end_time":"2022-10-19T01:39:42.039146","exception":false,"start_time":"2022-10-19T01:39:41.953634","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"patient_id = '1.2.826.0.1.3680043.12281'\nex_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":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:42.215151Z","iopub.status.busy":"2022-10-19T01:39:42.214657Z","iopub.status.idle":"2022-10-19T01:39:43.632798Z","shell.execute_reply":"2022-10-19T01:39:43.631606Z"},"papermill":{"duration":1.509104,"end_time":"2022-10-19T01:39:43.63547","exception":false,"start_time":"2022-10-19T01:39:42.126366","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"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 them is explained in the [discussion post](https://www.kaggle.com/competitions/rsna-2022-cervical-spine-fracture-detection/discussion/340612). We need to align the orientation of the segmentations with the one of the DICOM images.","metadata":{"papermill":{"duration":0.10079,"end_time":"2022-10-19T01:39:43.824943","exception":false,"start_time":"2022-10-19T01:39:43.724153","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Align orientation with images\nseg = seg[:, ::-1, ::-1].transpose(2, 1, 0)\nseg.shape","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:43.998156Z","iopub.status.busy":"2022-10-19T01:39:43.997139Z","iopub.status.idle":"2022-10-19T01:39:44.006397Z","shell.execute_reply":"2022-10-19T01:39:44.004781Z"},"papermill":{"duration":0.098023,"end_time":"2022-10-19T01:39:44.008971","exception":false,"start_time":"2022-10-19T01:39:43.910948","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n
\n

5.2 Exploring the masks

\n
","metadata":{"papermill":{"duration":0.084714,"end_time":"2022-10-19T01:39:44.179571","exception":false,"start_time":"2022-10-19T01:39:44.094857","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"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\nstart = 110\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.execute_input":"2022-10-19T01:39:44.352411Z","iopub.status.busy":"2022-10-19T01:39:44.351948Z","iopub.status.idle":"2022-10-19T01:39:45.773622Z","shell.execute_reply":"2022-10-19T01:39:45.772555Z"},"papermill":{"duration":1.511068,"end_time":"2022-10-19T01:39:45.776013","exception":false,"start_time":"2022-10-19T01:39:44.264945","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"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 really useful because we know that the fractures can only occur in these areas.\n\nThey also tell us which **vertebrae** are represented in the image. By looking at the unique values in each slice, we find a 0 for the background and another number for vertebrae, e.g. 2 for vertebrae C2. ","metadata":{"papermill":{"duration":0.085837,"end_time":"2022-10-19T01:39:45.95018","exception":false,"start_time":"2022-10-19T01:39:45.864343","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"np.unique(seg[116])","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:46.127954Z","iopub.status.busy":"2022-10-19T01:39:46.12654Z","iopub.status.idle":"2022-10-19T01:39:46.141078Z","shell.execute_reply":"2022-10-19T01:39:46.139949Z"},"papermill":{"duration":0.106335,"end_time":"2022-10-19T01:39:46.143663","exception":false,"start_time":"2022-10-19T01:39:46.037328","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Unfortunately, we don't have segmentations for the whole training set.","metadata":{"papermill":{"duration":0.088192,"end_time":"2022-10-19T01:39:46.31885","exception":false,"start_time":"2022-10-19T01:39:46.230658","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Number of cases with masks\nseg_paths = glob(f'{BASE}/segmentations/*')\nprint(\n f'Number of cases with segmentations: {len(seg_paths)}, ({np.round(100 * len(seg_paths) / len(train), 1)}%)'\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:46.495107Z","iopub.status.busy":"2022-10-19T01:39:46.494098Z","iopub.status.idle":"2022-10-19T01:39:46.513418Z","shell.execute_reply":"2022-10-19T01:39:46.512053Z"},"papermill":{"duration":0.109852,"end_time":"2022-10-19T01:39:46.51608","exception":false,"start_time":"2022-10-19T01:39:46.406228","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We could consider training a segmentation model like *UNet* to predict the segmentation masks for the remainder of the train and all of the test pictures.","metadata":{"papermill":{"duration":0.087557,"end_time":"2022-10-19T01:39:46.69125","exception":false,"start_time":"2022-10-19T01:39:46.603693","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"markdown","source":"# 6. Bounding Boxes\n\nLet's now analyse **only the slices/images containing bounding box annotations** (these slices contain 1 or more ruptures).\n\nWe are only given bounding boxes for a **subset** of the data. In particular, images of only **11.6%** of patients in the training set have any bounding box annotations.\n\nThis information is essential in determining the exact location of the fractures. To give bounding boxes for the whole training set, we may train an **object detection** algorithm.","metadata":{"papermill":{"duration":0.087902,"end_time":"2022-10-19T01:39:46.866939","exception":false,"start_time":"2022-10-19T01:39:46.779037","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"print(\n f'{clr.S}Number of Study Instances that contain bounding boxes:{clr.E} {train_bbox[\"StudyInstanceUID\"].nunique()}'\n f' ({np.round(100 * train_bbox[\"StudyInstanceUID\"].nunique() / len(train), 1)} %)'\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:47.046911Z","iopub.status.busy":"2022-10-19T01:39:47.045804Z","iopub.status.idle":"2022-10-19T01:39:47.054611Z","shell.execute_reply":"2022-10-19T01:39:47.053315Z"},"papermill":{"duration":0.101276,"end_time":"2022-10-19T01:39:47.057092","exception":false,"start_time":"2022-10-19T01:39:46.955816","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

6.1 Number of Bounding Boxes per Study Instance

\n
\n\n**Note**:\n\n- ~50% of all slices have about 25 or less bboxes per study.\n - ! ~16% of slices have less than 10 bboxes per study instance.\n- there are 6 studies that are outliers - with more than 100 bboxes per study.\n\n*! Keep in mind: number of bboxes == number of fractures*.","metadata":{"papermill":{"duration":0.086935,"end_time":"2022-10-19T01:39:47.405117","exception":false,"start_time":"2022-10-19T01:39:47.318182","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Plot 1\nplt.figure(figsize=(24, 12))\nsns.distplot(\n train_bbox['StudyInstanceUID'].value_counts().values,\n rug=True,\n bins=20,\n rug_kws={'color': my_colors[1]},\n kde_kws={'color': my_colors[8], 'lw': 5, 'alpha': 0.7},\n hist_kws={'histtype': 'step', 'linewidth': 3, 'alpha': 1, 'color': my_colors[1]},\n)\n\nplt.title('Number of Bounding Boxes per Study Instance', weight='bold', size=25)\nplt.xlabel('Number of BBoxes / Study Instance', size=18, weight='bold')\nplt.ylabel('Frequency')\n\n# Arrow\nstyle = 'Simple, tail_width=1, head_width=12, head_length=14'\nkw = dict(arrowstyle=style, color='black')\narrow = patches.FancyArrowPatch(\n (70, 0.022), (20, 0.018), connectionstyle='arc3,rad=.10', **kw\n)\nplt.gca().add_patch(arrow)\nplt.text(\n x=70,\n y=0.022,\n s='~50% of slices have less than 25 bboxes/study.',\n color='black',\n size=17,\n)\n\nstyle = 'Simple, tail_width=1, head_width=12, head_length=14'\nkw = dict(arrowstyle=style, color='black')\narrow = patches.FancyArrowPatch(\n (160, 0.008), (150, 0.001), connectionstyle='arc3,rad=.10', **kw\n)\nplt.gca().add_patch(arrow)\nplt.text(\n x=130, y=0.0085, s='Only 6 slices have 100+ bboxes/study.', color='black', size=17\n)\n\nsns.despine(right=True, top=True, left=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:47.582559Z","iopub.status.busy":"2022-10-19T01:39:47.581476Z","iopub.status.idle":"2022-10-19T01:39:47.925001Z","shell.execute_reply":"2022-10-19T01:39:47.923759Z"},"papermill":{"duration":0.435998,"end_time":"2022-10-19T01:39:47.928056","exception":false,"start_time":"2022-10-19T01:39:47.492058","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

6.2 Bounding boxes on images

\n
\n\nLet's now plot some bounding boxes on the CT scans.","metadata":{"papermill":{"duration":0.090681,"end_time":"2022-10-19T01:39:48.109118","exception":false,"start_time":"2022-10-19T01:39:48.018437","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Compute x2 and y2\ntrain_bbox['x2'] = train_bbox['x'] + train_bbox['width']\ntrain_bbox['y2'] = train_bbox['y'] + train_bbox['height']\n\n# Rename columns\ntrain_bbox.rename(columns={'x': 'x1', 'y': 'y1'}, inplace=True)\n\n# Change to int\ntrain_bbox['x1'] = train_bbox['x1'].apply(lambda x: int(x))\ntrain_bbox['x2'] = train_bbox['x2'].apply(lambda x: int(x))\ntrain_bbox['y1'] = train_bbox['y1'].apply(lambda x: int(x))\ntrain_bbox['y2'] = train_bbox['y2'].apply(lambda x: int(x))\n\ntrain_bbox.head(3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:48.294026Z","iopub.status.busy":"2022-10-19T01:39:48.293245Z","iopub.status.idle":"2022-10-19T01:39:48.333155Z","shell.execute_reply":"2022-10-19T01:39:48.331985Z"},"papermill":{"duration":0.134805,"end_time":"2022-10-19T01:39:48.335489","exception":false,"start_time":"2022-10-19T01:39:48.200684","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def show_dcm_bboxes(study_instance_UID, rows, cols):\n # Set defaults\n BASE_PATH = '../input/rsna-2022-cervical-spine-fracture-detection/train_images'\n N = rows * cols\n\n fig, axes = plt.subplots(nrows=rows, ncols=cols, figsize=(24, 25))\n\n # Get dataframe & .dcm paths\n bbox_df = (\n train_bbox[train_bbox['StudyInstanceUID'] == study_instance_UID]\n .reset_index(drop=True)\n .head(N)\n )\n png_paths = [\n f'{BASE_PATH}/{study_instance_UID}/{slice_no}.dcm'\n for slice_no in bbox_df['slice_number']\n ]\n\n for path, k in zip(png_paths, range(N)):\n dataset = pydicom.dcmread(path)\n img = apply_voi_lut(dataset.pixel_array, dataset)\n slice_no = png_paths[k].split('/')[-1].split('.')[0]\n\n # BBoxes\n bbox_list = (\n bbox_df.loc[k, 'x1'],\n bbox_df.loc[k, 'y1'],\n bbox_df.loc[k, 'x2'],\n bbox_df.loc[k, 'y2'],\n )\n x1, y1, x2, y2 = [int(x) for x in bbox_list]\n\n # Plot the image\n x_plot = k // cols\n y_plot = k % cols\n\n cv2.rectangle(img, (x1, y1), (x2, y2), (226, 56, 54), 5)\n axes[x_plot, y_plot].imshow(img, cmap='bone')\n axes[x_plot, y_plot].set_title(f'Slice: {slice_no}', fontsize=14, weight='bold')\n axes[x_plot, y_plot].axis('off')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:48.518401Z","iopub.status.busy":"2022-10-19T01:39:48.517984Z","iopub.status.idle":"2022-10-19T01:39:48.530131Z","shell.execute_reply":"2022-10-19T01:39:48.529147Z"},"papermill":{"duration":0.106226,"end_time":"2022-10-19T01:39:48.532431","exception":false,"start_time":"2022-10-19T01:39:48.426205","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Example 1","metadata":{"papermill":{"duration":0.102994,"end_time":"2022-10-19T01:39:48.726124","exception":false,"start_time":"2022-10-19T01:39:48.62313","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"show_dcm_bboxes(study_instance_UID='1.2.826.0.1.3680043.5783', rows=3, cols=3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:48.926985Z","iopub.status.busy":"2022-10-19T01:39:48.926506Z","iopub.status.idle":"2022-10-19T01:39:52.76457Z","shell.execute_reply":"2022-10-19T01:39:52.762216Z"},"papermill":{"duration":3.966486,"end_time":"2022-10-19T01:39:52.793644","exception":false,"start_time":"2022-10-19T01:39:48.827158","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Example 2","metadata":{"papermill":{"duration":0.127894,"end_time":"2022-10-19T01:39:53.052912","exception":false,"start_time":"2022-10-19T01:39:52.925018","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"show_dcm_bboxes(study_instance_UID='1.2.826.0.1.3680043.19778', rows=3, cols=3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:53.314819Z","iopub.status.busy":"2022-10-19T01:39:53.314133Z","iopub.status.idle":"2022-10-19T01:39:57.641938Z","shell.execute_reply":"2022-10-19T01:39:57.640701Z"},"papermill":{"duration":4.488056,"end_time":"2022-10-19T01:39:57.670214","exception":false,"start_time":"2022-10-19T01:39:53.182158","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Example 3","metadata":{"papermill":{"duration":0.169356,"end_time":"2022-10-19T01:39:58.016591","exception":false,"start_time":"2022-10-19T01:39:57.847235","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"show_dcm_bboxes(study_instance_UID='1.2.826.0.1.3680043.31077', rows=2, cols=2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:39:58.357199Z","iopub.status.busy":"2022-10-19T01:39:58.356509Z","iopub.status.idle":"2022-10-19T01:40:02.690647Z","shell.execute_reply":"2022-10-19T01:40:02.689153Z"},"papermill":{"duration":4.537316,"end_time":"2022-10-19T01:40:02.722289","exception":false,"start_time":"2022-10-19T01:39:58.184973","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 7. 3D CT Scans\n\nIf is also possible to visualize the CT scans in 3D using the *NIfTI* (*.nii*) segmentation files available in the dataset.\n\n> 🙏 This part is taken from [Seungwon Song](https://www.kaggle.com/songseungwon)'s notebook: [Break Time☕️ Just Enjoy drawing 3D cervical spine](https://www.kaggle.com/code/songseungwon/break-time-just-enjoy-drawing-3d-cervical-spine/notebook).","metadata":{"papermill":{"duration":0.255035,"end_time":"2022-10-19T01:40:05.853579","exception":false,"start_time":"2022-10-19T01:40:05.598544","status":"completed"},"pycharm":{"name":"#%% md\n"},"tags":[]}},{"cell_type":"code","source":"# Get .nii paths\nnii_paths = glob('../input/rsna-2022-cervical-spine-fracture-detection/segmentations/*')\nprint(f'{clr.S}Total paths in [segmentations] folder:{clr.E}', len(nii_paths))","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:40:06.361906Z","iopub.status.busy":"2022-10-19T01:40:06.361467Z","iopub.status.idle":"2022-10-19T01:40:06.368158Z","shell.execute_reply":"2022-10-19T01:40:06.367211Z"},"papermill":{"duration":0.263815,"end_time":"2022-10-19T01:40:06.370785","exception":false,"start_time":"2022-10-19T01:40:06.10697","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class DrawMaskSample:\n def __init__(self, nii_paths, i, ax):\n self.xyz_li = None\n self.i = i\n self.ax = ax\n self.nii_sample = nib.load(nii_paths[self.i]).get_fdata()\n self.get_xyz()\n self.draw_sample_3d()\n\n def get_xyz(self):\n self.xyz_li = []\n cnt = 0\n max_cnt = self.nii_sample.shape[-1]\n for iter_z, (iter_img) in enumerate(self.nii_sample.transpose(2, 0, 1)):\n for iter_x, iter_arr_y in enumerate(iter_img):\n iter_arr_y = np.where(iter_arr_y)[0]\n if len(iter_arr_y) >= 1:\n iter_arr_y = list(\n {iter_arr_y.max(), iter_arr_y.min()})\n xyz = [\n (iter_x, iter_y, iter_z)\n for iter_y in iter_arr_y\n if np.any(iter_y)\n ]\n self.xyz_li.append(xyz)\n cnt += 1\n\n def draw_sample_3d(self):\n xyz_matrix = np.array(list(itertools.chain.from_iterable(self.xyz_li)))\n X = xyz_matrix[:, 0]\n Y = xyz_matrix[:, 1]\n Z = xyz_matrix[:, 2]\n\n self.ax.scatter(X, Y, Z, s=1, alpha=0.05, color='#e3dac9')\n xlim, ylim, zlim = self.nii_sample.shape\n self.ax.set_xlim(0, xlim)\n self.ax.set_ylim(0, ylim)\n self.ax.set_zlim(0, zlim)\n self.ax.set_title(f'sample - ({self.i + 1})')\n self.ax.set_yticklabels([])\n self.ax.set_xticklabels([])\n self.ax.set_zticklabels([])\n self.ax.set_facecolor('#081921')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.execute_input":"2022-10-19T01:40:06.874473Z","iopub.status.busy":"2022-10-19T01:40:06.874048Z","iopub.status.idle":"2022-10-19T01:40:06.887197Z","shell.execute_reply":"2022-10-19T01:40:06.886367Z"},"papermill":{"duration":0.268557,"end_time":"2022-10-19T01:40:06.889453","exception":false,"start_time":"2022-10-19T01:40:06.620896","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(24, 24))\nfor i in range(9):\n ax = fig.add_subplot(int(f'33{i + 1}'), projection='3d')\n DrawMaskSample(nii_paths, i, ax)\n\nplt.tight_layout()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.execute_input":"2022-10-19T01:40:07.395793Z","iopub.status.busy":"2022-10-19T01:40:07.395353Z","iopub.status.idle":"2022-10-19T01:40:47.68244Z","shell.execute_reply":"2022-10-19T01:40:47.681086Z"},"papermill":{"duration":40.828768,"end_time":"2022-10-19T01:40:47.968807","exception":false,"start_time":"2022-10-19T01:40:07.140039","status":"completed"},"pycharm":{"name":"#%%\n"},"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 8. Summary\n\n- The images are CT scans, **not X-rays**.\n- The competition is about predicting the probability of a fracture for each of the seven cervical vertebrae and an overall probability of a fracture in the cervical spine. Therefore, **it is not an object detection problem**.\n - However, a part of the training data (235 out of 2019 study instances, i.e. 11.6% of patients) contains bounding boxes' annotations for the fractures. This information could be used to train an object detection model. The test set is not publicly available, so we don't know if it contains bounding boxes or not. In the worst case, the training set would have to be used to perform training, validation and testing.\n- The image size is mostly 512x512 pixels, but there is a fraction of images (around 0.3%) with different dimensions (512x519 and 768x768). These images should be resized to 512x512.","metadata":{}}]}