{"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":"# This is a fork of https://www.kaggle.com/code/andradaolteanu/rsna-fracture-detection-dicom-images-explore\n\nWith a fix for counting no of slices per case.\n\nIt was striking that it seemed that 25% of cases have less than 50 slices, since I recall no such 'short' studies, that's why I decided to check that.","metadata":{}},{"cell_type":"code","source":"# Libraries\nimport os\nimport re\nimport gc\nimport cv2\n# import wandb\nfrom PIL import Image\nimport random\nimport math\nimport shutil\nfrom glob import glob\nfrom tqdm import tqdm\nfrom pprint import pprint\nfrom time import time\nimport warnings\nimport itertools\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib as mpl\nfrom matplotlib import cm\nimport matplotlib.patches as patches\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nfrom matplotlib.offsetbox import AnnotationBbox, OffsetImage\nfrom matplotlib.colors import ListedColormap, LinearSegmentedColormap\nfrom matplotlib.patches import Rectangle\nfrom IPython.display import display_html\nplt.rcParams.update({'font.size': 16})\n\n# .dcm handling\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\n# Environment check\nwarnings.filterwarnings(\"ignore\")\n\nmy_colors = [\"#5EAFD9\", \"#449DD1\", \"#3977BB\", \n             \"#2D51A5\", \"#5C4C8F\", \"#8B4679\",\n             \"#C53D4C\", \"#E23836\", \"#FF4633\", \"#FF5746\"]\nCMAP1 = ListedColormap(my_colors)","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:51.873785Z","iopub.execute_input":"2022-08-09T09:41:51.874533Z","iopub.status.idle":"2022-08-09T09:41:52.741359Z","shell.execute_reply.started":"2022-08-09T09:41:51.874482Z","shell.execute_reply":"2022-08-09T09:41:52.740339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ⬇ Helper Functions","metadata":{}},{"cell_type":"code","source":"BASE = \"../input/rsna-2022-cervical-spine-fracture-detection\"","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.753374Z","iopub.execute_input":"2022-08-09T09:41:52.754489Z","iopub.status.idle":"2022-08-09T09:41:52.765851Z","shell.execute_reply.started":"2022-08-09T09:41:52.754436Z","shell.execute_reply":"2022-08-09T09:41:52.764639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_data():\n    '''Reads in all .csv files.'''\n    \n    train = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\n    train_bbox = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train_bounding_boxes.csv\")\n    test = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/test.csv\")\n    ss = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/sample_submission.csv\")\n    \n    return train, train_bbox, test, ss\n","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.769786Z","iopub.execute_input":"2022-08-09T09:41:52.770214Z","iopub.status.idle":"2022-08-09T09:41:52.777519Z","shell.execute_reply.started":"2022-08-09T09:41:52.770179Z","shell.execute_reply":"2022-08-09T09:41:52.776562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read in the data\ntrain, train_bbox, test, ss = read_data()","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.779256Z","iopub.execute_input":"2022-08-09T09:41:52.780066Z","iopub.status.idle":"2022-08-09T09:41:52.816673Z","shell.execute_reply.started":"2022-08-09T09:41:52.780018Z","shell.execute_reply":"2022-08-09T09:41:52.815637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. DICOM Metadata\n\n## 3.1 Retrieving Metadata for 1 file\n\n🦴 At the moment, the only information I would like to retrieve is the following:\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":{}},{"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\", \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    \n    return observation_data","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.81807Z","iopub.execute_input":"2022-08-09T09:41:52.818476Z","iopub.status.idle":"2022-08-09T09:41:52.827355Z","shell.execute_reply.started":"2022-08-09T09:41:52.818446Z","shell.execute_reply":"2022-08-09T09:41:52.826154Z"},"trusted":true},"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":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.828578Z","iopub.execute_input":"2022-08-09T09:41:52.829282Z","iopub.status.idle":"2022-08-09T09:41:52.858701Z","shell.execute_reply.started":"2022-08-09T09:41:52.829243Z","shell.execute_reply":"2022-08-09T09:41:52.857852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2 Create Metadata file","metadata":{}},{"cell_type":"code","source":"def get_metadata(max_n=10):\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(min(max_n, len(train)))):\n        dt = train.iloc[k, :]\n\n        if dt.patient_overall == 1:\n            # Get all .dcm paths for this Instance\n            dcm_paths = glob(f\"{BASE}/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=example.keys())\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":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.860319Z","iopub.execute_input":"2022-08-09T09:41:52.861072Z","iopub.status.idle":"2022-08-09T09:41:52.869365Z","shell.execute_reply.started":"2022-08-09T09:41:52.861036Z","shell.execute_reply":"2022-08-09T09:41:52.868228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> 🦴 **Note**: All data is available to download [here](https://www.kaggle.com/datasets/andradaolteanu/rsna-fracture-detection).","metadata":{}},{"cell_type":"code","source":"# Create and save the metadata\n# This cell takes ~ 1 hour to run\n# get_metadata()","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.870828Z","iopub.execute_input":"2022-08-09T09:41:52.871492Z","iopub.status.idle":"2022-08-09T09:41:52.885478Z","shell.execute_reply.started":"2022-08-09T09:41:52.871456Z","shell.execute_reply":"2022-08-09T09:41:52.884382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.3 Metadata Explore\n\nLet's play a bit with the data we have just created 😁","metadata":{}},{"cell_type":"code","source":"# Read in saved metadata\nmeta_train = pd.read_csv(\"../input/rsna-fracture-detection/meta_train.csv\")\nmeta_train[\"StudyInstanceUID\"] = meta_train[\"SOPInstanceUID\"].apply(lambda x: \".\".join(x.split(\".\")[:-2]))\n","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:52.887098Z","iopub.execute_input":"2022-08-09T09:41:52.887759Z","iopub.status.idle":"2022-08-09T09:41:54.015481Z","shell.execute_reply.started":"2022-08-09T09:41:52.887716Z","shell.execute_reply":"2022-08-09T09:41:54.013627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# !!!!!!!!!!!!!!!!!!!\n\n# Here, we need to count number of rows for each StudyInstanceUID (patient) not for each slice_idx.\n\n# !!!!!!!!!!!!!!!!!!!\n\ndata = meta_train[\"StudyInstanceUID\"].value_counts().reset_index()\ndata.columns = [\"StudyInstanceUID\", \"count\"]\ndata","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:54.017493Z","iopub.execute_input":"2022-08-09T09:41:54.01803Z","iopub.status.idle":"2022-08-09T09:41:54.077156Z","shell.execute_reply.started":"2022-08-09T09:41:54.017978Z","shell.execute_reply":"2022-08-09T09:41:54.076027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot 1\nplt.figure(figsize=(24, 12))\nsns.distplot(data[\"count\"], rug=True, 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\nsns.despine(right=True, top=True, left=True);","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:54.080291Z","iopub.execute_input":"2022-08-09T09:41:54.080683Z","iopub.status.idle":"2022-08-09T09:41:54.528736Z","shell.execute_reply.started":"2022-08-09T09:41:54.080648Z","shell.execute_reply":"2022-08-09T09:41:54.527773Z"},"trusted":true},"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 i in range(len(data))]\nperc = round(data[data[\"count\"]<=10].shape[0]/len(data), 3)*100\n\nplt.figure(figsize=(24, 12))\nsns.scatterplot(data=data, x=\"count\", y=\"y\", size=\"count\", alpha=0.65, sizes=(500, 2000),\n               hue=\"count\", palette=CMAP1)\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\n# plt.axvline(x=50, linestyle = '--', color=\"black\", lw=4)\n# plt.text(x=150, y=85, s=f\"~30% of the data is here\", color=\"black\", size=17, weight=\"bold\")\n# plt.text(x=50, y=0, s=f\"<- lower than 50 slices\", color=\"black\", size=14, weight=\"bold\")\n# plt.arrow(x=350, y=83, dx=-320, dy=0, color=\"black\", lw=4, \n#           head_width=2, head_length=8)\n\n# plt.axvline(x=900, linestyle = '--', color=\"black\", lw=4)\n# plt.text(x=600, y=85, s=f\"~25% of the data is here\", color=\"black\", size=17, weight=\"bold\")\n# plt.text(x=735, y=0, s=f\"greater than 900 slices ->\", color=\"black\", size=14, weight=\"bold\")\n# plt.arrow(x=600, y=83, dx=320, dy=0, color=\"black\", lw=4, \n#           head_width=2, head_length=8)\n\nplt.legend('',frameon=False)\nsns.despine(right=True, top=True, left=True);","metadata":{"execution":{"iopub.status.busy":"2022-08-09T09:41:54.530129Z","iopub.execute_input":"2022-08-09T09:41:54.530659Z","iopub.status.idle":"2022-08-09T09:41:55.049547Z","shell.execute_reply.started":"2022-08-09T09:41:54.530625Z","shell.execute_reply":"2022-08-09T09:41:55.048496Z"},"trusted":true},"execution_count":null,"outputs":[]}]}