{"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":"# Table of contents\n\n1. [Prepare to start](#prepare)\n2. [How are the targets distributed?](#targets)\n3. [Working with dicom files](#work_with_dicom)\n    * [Exploring an example dicom file](#example_dicom)\n    * [Transforming the data to Hounsfield Units (HU)](#hu)\n    * [Patient 1.2.826.0.1.3680043.10001](#example_patient)\n    * [The distribution of voxel sizes](#voxel_sizes)\n4. [Ideas for validation](#validation)","metadata":{}},{"cell_type":"code","source":"show_3d_reconstruction = False","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prepare to start <a class=\"anchor\" id=\"prepare\"></a>\n\nLet's load the packages first...","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport seaborn as sns\nsns.set()\n\nimport pydicom\n\nimport cv2\n\nfrom skimage import measure \nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\n\nfrom os import listdir\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-07T18:10:16.452193Z","iopub.execute_input":"2022-08-07T18:10:16.452629Z","iopub.status.idle":"2022-08-07T18:10:16.462837Z","shell.execute_reply.started":"2022-08-07T18:10:16.452594Z","shell.execute_reply":"2022-08-07T18:10:16.461846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And the data! ;-)","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\ntrain.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:09.47249Z","iopub.execute_input":"2022-08-07T18:06:09.472959Z","iopub.status.idle":"2022-08-07T18:06:09.504893Z","shell.execute_reply.started":"2022-08-07T18:06:09.472913Z","shell.execute_reply":"2022-08-07T18:06:09.503903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:09.506825Z","iopub.execute_input":"2022-08-07T18:06:09.507686Z","iopub.status.idle":"2022-08-07T18:06:09.528936Z","shell.execute_reply.started":"2022-08-07T18:06:09.507637Z","shell.execute_reply":"2022-08-07T18:06:09.527733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.StudyInstanceUID.nunique()","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:09.531979Z","iopub.execute_input":"2022-08-07T18:06:09.533476Z","iopub.status.idle":"2022-08-07T18:06:09.550004Z","shell.execute_reply.started":"2022-08-07T18:06:09.533427Z","shell.execute_reply":"2022-08-07T18:06:09.548732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/test.csv\")\ntest.shape","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:09.551865Z","iopub.execute_input":"2022-08-07T18:06:09.55234Z","iopub.status.idle":"2022-08-07T18:06:09.578501Z","shell.execute_reply.started":"2022-08-07T18:06:09.552295Z","shell.execute_reply":"2022-08-07T18:06:09.576464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:09.58029Z","iopub.execute_input":"2022-08-07T18:06:09.58067Z","iopub.status.idle":"2022-08-07T18:06:09.594875Z","shell.execute_reply.started":"2022-08-07T18:06:09.580636Z","shell.execute_reply":"2022-08-07T18:06:09.593458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# How are the targets distributed? <a class=\"anchor\" id=\"targets\"></a>","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(2,4,figsize=(25,12))\n\nsns.countplot(train.patient_overall, ax=ax[0,0], palette=\"Reds_r\")\nax[0,0].set_title(\"Patient overall target count\");\nsns.countplot(train.C1, ax=ax[0,1], palette=\"Blues_r\")\nax[0,1].set_title(\"C1 target count\");\nsns.countplot(train.C2, ax=ax[0,2], palette=\"Blues_r\")\nax[0,2].set_title(\"C2 target count\");\nsns.countplot(train.C3, ax=ax[0,3], palette=\"Blues_r\")\nax[0,3].set_title(\"C3 target count\");\n\nsns.countplot(train.C4, ax=ax[1,0], palette=\"Blues_r\")\nax[1,0].set_title(\"C4 target count\");\nsns.countplot(train.C5, ax=ax[1,1], palette=\"Blues_r\")\nax[1,1].set_title(\"C5 target count\");\nsns.countplot(train.C6, ax=ax[1,2], palette=\"Blues_r\")\nax[1,2].set_title(\"C6 target count\");\nsns.countplot(train.C7, ax=ax[1,3], palette=\"Blues_r\")\nax[1,3].set_title(\"C7 target count\");\n\nfor n in range(4):\n    ax[0,n].set_xlabel(\"\")\n    ax[1,n].set_xlabel(\"\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-07T18:06:09.597354Z","iopub.execute_input":"2022-08-07T18:06:09.598269Z","iopub.status.idle":"2022-08-07T18:06:10.592196Z","shell.execute_reply.started":"2022-08-07T18:06:09.598206Z","shell.execute_reply":"2022-08-07T18:06:10.59123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Insights\n\n* The patient overall target is roughly balanced compared to the underlysing subtarget structure. \n* The overall target tells us whether any of the vertebrae are fructured whereas the others show us in detail which of the vertebral columns are affected.\n\nThis opens the possibility that a patient has more than one fracture. Let's count the number of fractures over all vertebral columns C1 to C7 to find out how multiple fractures are distributed:","metadata":{}},{"cell_type":"code","source":"fracture_counts = train[train.patient_overall==1].drop([\"StudyInstanceUID\", \"patient_overall\"], axis=1).sum(axis=1)\nplt.figure(figsize=(10,5))\nsns.countplot(fracture_counts, palette=\"Blues_r\")\nplt.title(\"How are the number of fractures per patient distributed?\")\nplt.xlabel(\"Number of factures\")\nplt.ylabel(\"Patient counts\");","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:10.593373Z","iopub.execute_input":"2022-08-07T18:06:10.594075Z","iopub.status.idle":"2022-08-07T18:06:10.824039Z","shell.execute_reply.started":"2022-08-07T18:06:10.59404Z","shell.execute_reply":"2022-08-07T18:06:10.823063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A lot of patients have only one fracture. But for those with multiple fractures it can be interesting to explore whether specific combinations are more likely than others. Perhaps I will come back to that part later.","metadata":{}},{"cell_type":"markdown","source":"# Working with dicom files <a class=\"anchor\" id=\"work_with_dicom\"></a>\n\n## Exploring an example dicom file <a class=\"anchor\" id=\"explore_dicom\"></a>\n","metadata":{}},{"cell_type":"code","source":"example = \"../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.10001/104.dcm\"\nexample_file = pydicom.dcmread(example)\nexample_file","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:10.825335Z","iopub.execute_input":"2022-08-07T18:06:10.826371Z","iopub.status.idle":"2022-08-07T18:06:10.855561Z","shell.execute_reply.started":"2022-08-07T18:06:10.826332Z","shell.execute_reply":"2022-08-07T18:06:10.854561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Our image data is stored in the array \"Pixel Data\". \n* The Rows and the Columns attribute yield us the image size.\n* The pixel spacing and the slice thickness tell us how much physical distance is covered by one pixel in mm. We can use them to compute the voxel size which holds the volumne of one pixel.  \n* The thrid value of the Image Position (Patient) helps us to order the files the right way from bottom to top.","metadata":{}},{"cell_type":"markdown","source":"## Transforming the data to Hounsfield Units (HU) <a class=\"anchor\" id=\"hu\"></a>\n\n* Before we can work with our pixel array data, **we need to transform its values to Hounsfield Units**. \n* A CT-scan holds information about the radiodensity of an object that was exposed to x-rays.  \n* The spectral composition of the x-rays depends on the measurement settings like the tube voltage and aquisition parameters. **To make different measurements comparable we need to normalizing to values of water (HU:0) and air (HU:-1000).** \n* This can be done by using the Rescale Intercept and the Rescale Slope stored in our dicom file.","metadata":{}},{"cell_type":"markdown","source":"Here is how you can do it! ;-)","metadata":{}},{"cell_type":"code","source":"image = example_file.pixel_array.flatten()\nrescaled_image = image * example_file.RescaleSlope + example_file.RescaleIntercept","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:10.858293Z","iopub.execute_input":"2022-08-07T18:06:10.858976Z","iopub.status.idle":"2022-08-07T18:06:10.869825Z","shell.execute_reply.started":"2022-08-07T18:06:10.858937Z","shell.execute_reply":"2022-08-07T18:06:10.868511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a look now at the raw pixel array distribution and after HU transformation:","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,5))\nax[0].imshow(example_file.pixel_array, cmap=\"bone\")\nax[0].axis(\"off\")\nsns.distplot(image.flatten(), ax=ax[1]);\nsns.distplot(rescaled_image.flatten(), ax=ax[2])\nax[1].set_title(\"Raw pixel array distributions\")\nax[2].set_title(\"HU unit distributions\");","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-07T18:06:10.871601Z","iopub.execute_input":"2022-08-07T18:06:10.872344Z","iopub.status.idle":"2022-08-07T18:06:13.440684Z","shell.execute_reply.started":"2022-08-07T18:06:10.872296Z","shell.execute_reply":"2022-08-07T18:06:13.439282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Insights\n\n* Take a look at the raw pixel array distribution. You can see that we have a lot of values at 0 (which belongs to air) and at 1000 (which belongs to water). Furthermore we have a *strange* peak at -2000. These values belong to the \"outside scanner\"-region which is colored black in the raw ct-scan image. It's also possible that these values have been set to even lower values like -3000.  \n* After transforming to Hounsfield Units the distribution has a peak at -3000 (outside scanner), at roughly -1000 (air) and a peak at 0 (water). ","metadata":{}},{"cell_type":"markdown","source":"Now we need to make a decision. Do we like to set the outside scanner region to values of air or not?! Instead of setting its values to air we could also set it to -2000 by default and add an aperture augmentation to increase variability during training. Here is some code you can use to transform to Hounsfield Units and setting the outside scanner region to air:","metadata":{}},{"cell_type":"code","source":"def set_outside_scanner_to_air(raw_pixelarrays): \n    # Let's threshold between air (0) and this default (-2000) using -1000\n    raw_pixelarrays[raw_pixelarrays <= -1000] = 0\n    return raw_pixelarrays\n\ndef transform_to_hu(slices):\n    images = np.stack([file.pixel_array for file in slices])\n    images = images.astype(np.int16)\n\n    images = set_outside_scanner_to_air(images)\n    \n    # convert to HU\n    for n in range(len(slices)):\n        \n        intercept = slices[n].RescaleIntercept\n        slope = slices[n].RescaleSlope\n        \n        if slope != 1:\n            images[n] = slope * images[n].astype(np.float64)\n            images[n] = images[n].astype(np.int16)\n            \n        images[n] += np.int16(intercept)\n    \n    return np.array(images, dtype=np.int16)\n\ndef load_scans(study_id_path):\n    slices = [pydicom.dcmread(study_id_path + \"/\" + file) for file in listdir(study_id_path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    \n    hu_images = transform_to_hu(slices)\n    return hu_images","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:54:52.631128Z","iopub.execute_input":"2022-08-07T18:54:52.63154Z","iopub.status.idle":"2022-08-07T18:54:52.641889Z","shell.execute_reply.started":"2022-08-07T18:54:52.631506Z","shell.execute_reply":"2022-08-07T18:54:52.640614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* The dicom file contains much more information that is worth it to explore but to get started with a model, it might be sufficient to use the images that had been transformed to Hounsfield Units. \n","metadata":{}},{"cell_type":"markdown","source":"## Patient 1.2.826.0.1.3680043.10001 <a class=\"anchor\" id=\"example_patient\"></a>\n\nLet's stay with our example patient *1.2.826.0.1.3680043.10001*. ","metadata":{}},{"cell_type":"code","source":"train[train.StudyInstanceUID==\"1.2.826.0.1.3680043.10001\"]","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:09:11.406328Z","iopub.execute_input":"2022-08-07T18:09:11.406802Z","iopub.status.idle":"2022-08-07T18:09:11.420988Z","shell.execute_reply.started":"2022-08-07T18:09:11.406763Z","shell.execute_reply":"2022-08-07T18:09:11.420128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient_path = \"../input/rsna-2022-cervical-spine-fracture-detection/train_images/1.2.826.0.1.3680043.10001/\"","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:54:00.575027Z","iopub.execute_input":"2022-08-07T18:54:00.575973Z","iopub.status.idle":"2022-08-07T18:54:00.580868Z","shell.execute_reply.started":"2022-08-07T18:54:00.575931Z","shell.execute_reply":"2022-08-07T18:54:00.579912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient_files = listdir(patient_path)","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:54:00.735344Z","iopub.execute_input":"2022-08-07T18:54:00.735759Z","iopub.status.idle":"2022-08-07T18:54:00.74326Z","shell.execute_reply.started":"2022-08-07T18:54:00.735726Z","shell.execute_reply":"2022-08-07T18:54:00.742069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ok, before we start, we should make sure that all slices of the patient show the same image size and the same slice thickness, pixelspacing and therefore the same voxel size.","metadata":{}},{"cell_type":"code","source":"def check_entry(l):\n    if len(np.unique(l))==1:\n        l = np.unique(l)[0]\n    else:\n        l = None\n    return l\n\ndef get_patient_meta(patient_path):\n    rows = []\n    columns = []\n    pixelsp_x = []\n    pixelsp_y = []\n    thicknesses = []\n    patient_files = listdir(patient_path)\n    for file in patient_files:\n        path = \"{}{}\".format(patient_path, file)\n        dcm_file = pydicom.dcmread(path)\n        rows.append(dcm_file.Rows)\n        columns.append(dcm_file.Columns)\n        pixelsp_x.append(dcm_file.PixelSpacing[0])\n        pixelsp_y.append(dcm_file.PixelSpacing[1])\n        thicknesses.append(dcm_file.SliceThickness)\n    \n    rows = check_entry(rows)\n    columns = check_entry(columns)\n    pixelsp_x = check_entry(pixelsp_x)\n    pixelsp_y = check_entry(pixelsp_y)\n    thicknesses = check_entry(thicknesses)\n    \n    return {\"image_size\": [rows, columns],\n            \"voxel\": [pixelsp_x, pixelsp_y, thicknesses]}","metadata":{"execution":{"iopub.status.busy":"2022-08-07T19:45:17.078203Z","iopub.execute_input":"2022-08-07T19:45:17.079154Z","iopub.status.idle":"2022-08-07T19:45:17.091898Z","shell.execute_reply.started":"2022-08-07T19:45:17.079094Z","shell.execute_reply":"2022-08-07T19:45:17.089894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"get_patient_meta(patient_path)","metadata":{"execution":{"iopub.status.busy":"2022-08-07T19:45:18.534182Z","iopub.execute_input":"2022-08-07T19:45:18.534638Z","iopub.status.idle":"2022-08-07T19:45:19.156814Z","shell.execute_reply.started":"2022-08-07T19:45:18.5346Z","shell.execute_reply":"2022-08-07T19:45:19.155537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Great! For this patient we can be sure that all slices have the same properties! :-)\n\nPlease keep in mind that this does not need to be true for all patients. Usally these properties differ from patient to patient as it depends on the scanner type, the doctors personal preferences and the specifications of the clinic.  ","metadata":{}},{"cell_type":"markdown","source":"## The distribution of voxel sizes <a class=\"anchor\" id=\"voxel_sizes\"></a>\n\n**To be continued** :-)","metadata":{}},{"cell_type":"code","source":"for study_id in train.StudyInstanceUID:\n    path = \"../input/rsna-2022-cervical-spine-fracture-detection/train_images/{}/\".format(study_id)\n    all_patient_files = listdir(path)\n    train.loc[train.StudyInstanceUID==study_id, \"num_slices\"] = len(all_patient_files)\n    train.loc[train.StudyInstanceUID==study_id, \"num_slices\"] = len(all_patient_files)","metadata":{"execution":{"iopub.status.busy":"2022-08-07T19:34:28.274602Z","iopub.execute_input":"2022-08-07T19:34:28.27497Z","iopub.status.idle":"2022-08-07T19:34:32.903317Z","shell.execute_reply.started":"2022-08-07T19:34:28.274936Z","shell.execute_reply":"2022-08-07T19:34:32.902131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-07T19:34:28.254417Z","iopub.execute_input":"2022-08-07T19:34:28.255189Z","iopub.status.idle":"2022-08-07T19:34:28.272526Z","shell.execute_reply.started":"2022-08-07T19:34:28.255121Z","shell.execute_reply":"2022-08-07T19:34:28.271061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.distplot(train.num_slices)","metadata":{"execution":{"iopub.status.busy":"2022-08-07T19:38:17.327181Z","iopub.execute_input":"2022-08-07T19:38:17.327617Z","iopub.status.idle":"2022-08-07T19:38:17.748053Z","shell.execute_reply.started":"2022-08-07T19:38:17.327583Z","shell.execute_reply":"2022-08-07T19:38:17.746357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ideas for validation","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}