{"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":"code","source":"#Importing libraries \nimport pandas as pd\nimport numpy as np\nimport pydicom \n\nimport os\nimport random\n\n#Visualisation \nimport matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings(\"ignore\");\n%matplotlib inline","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_PATH = '../input/rsna-intracranial-hemorrhage-detection/rsna-intracranial-hemorrhage-detection/'","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"os.listdir(INPUT_PATH)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# **Exploratory analysis**\n\nLet us try exploring this dataset, to get some insights and understand our data better","metadata":{}},{"cell_type":"code","source":"sub = pd.read_csv(INPUT_PATH+\"stage_2_sample_submission.csv\")\nsub.head(10)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So Image IDs of form ID_SUBTYPE, which means we would have to make predictions for each subtype under a image ID. \n`any` --> indicates there is at least one subtype present telling us that patient has IH or not\n`Label` --> indicates probability of presence","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv(INPUT_PATH+\"stage_2_train.csv\")\n# train_df.head(10)\nlabels =  train_df.Label.values","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The training dataset is provided as set image`Id` and **multiple labels**, one for each of the subtypes of hemorrhage along with an addition lable for `any`(will be true of any of the subtype labels in true). So this is a **multilable classification task**.\n\nLets split the ID into columns of images and the corresponding diagnosis(subtype)","metadata":{}},{"cell_type":"code","source":"train_df = train_df.ID.str.rsplit(\"_\",n=1,expand=True)\ntrain_df.loc[:, \"label\"] = labels\ntrain_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = train_df.rename({0 : \"image\",1 : \"subtype\"}, axis=1)","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Look into target distribution","metadata":{}},{"cell_type":"code","source":"#find the count of targets under each subtype\nsubtype_count = train_df.groupby(\"subtype\").label.value_counts().unstack()\nsubtype_count","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#calculating the % target distribution across each subtype\nsubtype_count_per =  subtype_count.loc[:,1]/train_df.groupby(\"subtype\").size() *100\n\nmulti_target_count = train_df.groupby(\"image\").label.sum()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Helper function\ndef random_colors(num_colors : int):\n    colors = []\n    for i in range(num_colors):\n        colors.append('#'+''.join([random.choice('0123456789ABCDEF') for j in range(6)]))\n    return colors ","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lets ask questions and try finding answers to it through visualisations :\n1. How many positive and negative targets do we see in the training dataset?\n2. What is the target distribution across each of the labels ?\n3. Is the dataset imbalanced ?","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(3,1, figsize=(30,50))\n\nsns.countplot(train_df.label,ax=ax[0], palette=random_colors(2))\nax[0].set_xlabel(\"Target\", fontsize=40)\nax[0].tick_params(axis='x', labelsize=25 ) \nax[0].tick_params(axis='y', labelsize=25 ) \nax[0].set_title(\"Number of positive and negative targets\",fontsize=40)\n\n\nsns.countplot(x=\"subtype\", hue=\"label\", data=train_df, ax=ax[1], palette=random_colors(6))\nax[1].set_xlabel(\"Number of targets per image\",fontsize=40)\nax[1].set_ylabel(\"Frequency\",fontsize=40)\nax[1].tick_params(axis='x', labelsize=25 ) \nax[1].tick_params(axis='y', labelsize=25 ) \nax[1].set_title(\"Target distrubution\",fontsize=40)\n\nsns.barplot(subtype_count_per.index, subtype_count_per.values, ax=ax[2], palette=random_colors(6))\nplt.xticks(rotation=45)\nax[2].set_ylabel(\"% of positive(1) occurences\",fontsize=40)\nax[2].tick_params(axis='x', labelsize=25 ) \nax[2].tick_params(axis='y', labelsize=25 ) \nax[2].set_title(\"Imbalance in target distrubution\",fontsize=40)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So what do we make of it ?\n\n- Less no of positive target values\n- Epidural type has very few positive occurences(<1%)\n- Highly Imbalanced","metadata":{}},{"cell_type":"markdown","source":"### Basic Checks before exploring images..","metadata":{}},{"cell_type":"markdown","source":"Compare counts of images provided in training dataset with the training files given to check \neverything is fine before proceeeding..","metadata":{}},{"cell_type":"code","source":"#Count of images in training dataset\ntrain_df.image.nunique()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_files = os.listdir(INPUT_PATH+\"stage_2_train\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Actual no. of image files \nlen(train_files)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ok. So no issues in training set.","metadata":{}},{"cell_type":"code","source":"test_files = os.listdir(INPUT_PATH+\"stage_2_test\")\nlen(test_files)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(train_files)/len(test_files)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So 6.2times more images in train dataset than test dataset","metadata":{}},{"cell_type":"markdown","source":"## Lets look into DICOM files\n\n### **What is a DICOM ?**\nDicom is a format that has metadata, as well as Pixeldata attached to it. Below I extract some basic info with an image. You will know about the gender and age of the patient, as well as info how the image is sampled and generated. \n\nSo lets look at some samples from our dataset","metadata":{}},{"cell_type":"code","source":"train_files[:5]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Below are some of the slices from CT scans that are stored as pixel data in DICOM files","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(40,20))\ncolumn= 5; rows = 4\nfor i in range(1, column*rows +1):\n    dcm = pydicom.dcmread(INPUT_PATH+\"stage_2_train/\"+train_files[i])\n    fig.add_subplot(rows, column, i)\n    plt.imshow(dcm.pixel_array, cmap=plt.cm.bone)\n    fig.add_subplot","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Look into the meta data that comes with the DICOM File, look for insights that can help us during processing ","metadata":{}},{"cell_type":"code","source":"print(dcm)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we see there are details about the sampling along with the patients details. \nSome of it(like Window center, Window width, Rescale Intercept) can help at better pre-processing of these DICOM files.\n\nLets look the pixel data of one sample and find out the shape of these images","metadata":{}},{"cell_type":"code","source":"image = dcm.pixel_array\nprint(type(image)) #format in which pixel data is stored\nprint(image.dtype) #datatype of the pixel values\nprint(image.shape) #shape of image(wxh)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So images are of 512x512, we'll downsample them later to deal with the large training set","metadata":{}},{"cell_type":"code","source":"plt.imshow(image, cmap=plt.cm.bone)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"DICOM images typically contain between 12–16 bits/pixel, which corresponds to approximately 4,096 to 65,536 shades of gray. But most regular computer screens are often limited to 8 bits or 256 shades of gray. \n\nMost images like the one above **display a wide range of tissue densities**(ranging from -1000HU(air) to +1000HU(bone)), but as mentioned above a **computer screen  can only display 256 shades of gray with our eye detecting only about a 6% change** in grayscale ","metadata":{}},{"cell_type":"markdown","source":"**Math around it :**\n Eye can detect only 6% change in grayscale, so `100/6 = 17 shades of gray`\n To display a DICOM(having range of approx, 2000HU) image on computer screen(can only display 256 shades pf gray) = `2000/256 = 8 --> each shade of gray would have diff of 8HU`\n \n Therefore, each variation would vary by `256/17*9 = 120HU`\n \n BUT, the difference between normal and pathologically altered tissue is usually a lot less than 120 HU \n ","metadata":{"_kg_hide-input":true}},{"cell_type":"markdown","source":"So what to do? **Windowing !**","metadata":{}},{"cell_type":"markdown","source":"# **Getting into data preprocessing**","metadata":{}},{"cell_type":"markdown","source":"The point of applying windows is **to focus down the 256 shades of grey into a narrow region of HU(Hounsfiled units) that contain the relevnat densities of tissues we are interested** in while diagonising.","metadata":{}},{"cell_type":"code","source":"## A function to correct pixel data and rescale intercercepts ob 12 bit images\ndef dcm_correction(dcm_img):\n        x = dcm_img.pixel_array + 1000\n        px_mode = 4096\n        x[x >= px_mode] = x[x >= px_mode] - px_mode #if there are extra bits in 12-bit grayscale(<=4096)\n        dcm_img.PixelData = x.tobytes()\n        dcm_img.RescaleIntercept = -1000 #setting a common value across all 12-bit US images","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_size = []\nfor i in range(len(INPUT_PATH+\"stage_2_train/\")):\n    dicom = pydicom.dcmread(INPUT_PATH+\"stage_2_train/\"+train_files[i])\n    \n    if dicom.BitsStored == 12:\n        diff_size.append(dicom)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(diff_size)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_size[1]","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_size[1].pixel_array","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dcm_correction(diff_size[1])\ndiff_size[1]","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_rescale = []\nfor i in range(len(INPUT_PATH+\"stage_2_train/\")):\n    dicom = pydicom.dcmread(INPUT_PATH+\"stage_2_train/\"+train_files[i])\n    \n    if (int(dicom.RescaleIntercept) != -1024):\n        diff_rescale.append(dicom)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_rescale[0].pixel_array","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dcm_correction(diff_rescale[0])","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_rescale[0].pixel_array","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_rescale","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"need_correct = []\nfor i in range(len(INPUT_PATH+\"stage_2_train/\")):\n    dicom = pydicom.dcmread(INPUT_PATH+\"stage_2_train/\"+train_files[i])\n    \n    if (dicom.BitsStored == 12) and (dicom.PixelRepresentation == 0):\n        need_correct.append(dicom)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(need_correct)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## **Windowing**\n\nWindowing, also known as grey-level mapping or contrast enhancement is the process in which the CT image greyscale component of an image is manipulated via the CT numbers; doing this will change the appearance of the picture **to highlight particular structures**\n\nHere's where some of the DICOM meta data comes to help -\n\n`Window width` also know as the contrast --> is the measure of the range of CT numbers that an image contains\n\n`Window center` also known as brightness -->  is the midpoint of the range of the CT numbers displayed; window level is decreased the CT image will be brighter and vice versa.\n\n\nThese two could be used to calculate the upper and lower grey levels, to produce different kinds of windows based on the kind of diagnosis.\n","metadata":{}},{"cell_type":"code","source":"#Systemic/linear windowing\ndef window_image(dcm, window_center, window_width):\n    if (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0):\n        dcm_correction(dcm)\n\n    img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept #reconstructing the image from pixels\n    img_min = window_center - window_width // 2 #lowest visible value\n    img_max = window_center + window_width // 2 #highest visible value\n    img = np.clip(img, img_min, img_max)\n\n    return img","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**So what happening above ?**\n\nIf the DICOM file is of 12-bit type(41 outliers), then we correct them before generating our windows.\n\nWe then clip the pixel intensities between the lowest and hishest visisble values, to focus only on a narrow region where the abnormality might be present. This means that every pixel value greater than the `img_max` will show up as **white** and belowe `img_min` will show up as **black**\n\nSince each window highlights particular ranges, it makes it easier for a radiologists(the DL system in our case) to see if there are any changes between normal and pathologically altered tissue. So based on the diagnosis, the model would learn to look at only certain windows of tissue desities(features). ","metadata":{}},{"cell_type":"markdown","source":"Lets now try to get a picture of how windows would help us","metadata":{}},{"cell_type":"code","source":"TRAIN_PATH = INPUT_PATH+\"stage_2_train/\"\nTEST_PATH = INPUT_PATH+\"stage_2_test/\"\n\ndef view_images(image, title=''):\n\n    dcm = pydicom.dcmread(os.path.join(TRAIN_PATH,image[3]+'.dcm'))\n\n    fig, (ax1, ax2, ax3, ax4) = plt.subplots(4,1, figsize=(10,24))\n    \n    ax1.set_title(\"Default window\")\n    ax1.imshow(dcm.pixel_array, cmap=plt.cm.bone)\n    \n    ax2.set_title(\"Brain window\")\n    brain_img = window_image(dcm, 40, 80)\n    ax2.imshow(brain_img, cmap=plt.cm.bone)\n    \n    ax3.set_title(\"Subdural window\")\n    subdural_img = window_image(dcm, 80, 200)\n    ax3.imshow(subdural_img, cmap=plt.cm.bone)\n#     ax3.annotate('', xy=(150, 380), xytext=(120, 430),\n#             arrowprops=dict(facecolor='red', shrink=0.05),\n#             )\n#     ax3.annotate('', xy=(220, 430), xytext=(190, 500),\n#             arrowprops=dict(facecolor='red', shrink=0.05),\n#             )\n    \n    ax4.set_title(\"Soft Tissue window\")\n    soft_img = window_image(dcm, 40, 380)\n    ax4.imshow(soft_img, cmap=plt.cm.bone)\n    \n    for ax in fig.axes:\n        ax.axis(\"off\")\n        \n    fig.suptitle(title)\n    plt.show()\n    ","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Too much info! what are these different windows?**\n\n[radiopedia.org](https://radiopaedia.org/articles/ct-head-an-approach?lang=gb) shows a typical workflow and well thought process that a radiologist takes when given a task to detect any abnormalties on CT scan of the brain.\n\nFor head CT, bone window and brain window are two important window settings.However, the details of soft tissues such as brain, that shows density lower than that of bones, are lost in the bone window setting. Brain window is the most frequently used setting, and the majority of evaluations of brain abnormality are done using this window setting.\n\nWhat we understand is that, while the brain matter window is able to pick most abnormalities it might cause to miss some diagnosis. So while diagnising something like a hemorrhage we need to look into other windows like the subdural that focus more on the subdural hematoma.\n\n#### **Subdurals could be tricky..** ###\n\nIf you check their definition, they usually are right next to the skull, longish in shape and follows the curvature of the skull. Hence , if you look through a brain window you might miss out on these.. hence it adviced to incoporate a subdural window.\n","metadata":{}},{"cell_type":"code","source":"train_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_images(train_df[(train_df[\"subtype\"] == 'epidural') & (train_df['label'] == 1)][:10].image.values, title='Images wth epidural')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_images(train_df[(train_df[\"subtype\"] == 'subdural') & (train_df['label'] == 1)][:10].image.values, title='Images wth subdural')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_images(train_df[(train_df[\"subtype\"] == 'subarachnoid') & (train_df['label'] == 1)][:10].image.values, title='Images wth subarachnoid')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_images(train_df[(train_df[\"subtype\"] == 'intraventricular') & (train_df['label'] == 1)][:10].image.values, title='Images wth intraventricular')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_images(train_df[(train_df[\"subtype\"] == 'intraparenchymal') & (train_df['label'] == 1)][:10].image.values, title='Images wth intraparenchymal')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Combining the windows ","metadata":{}},{"cell_type":"markdown","source":"Lets take a random sample from the dataset to show the entire preprocessing would be","metadata":{}},{"cell_type":"code","source":"test_case = os.path.join(TRAIN_PATH,'ID_12cadc6af.dcm')\n\ntest_data = pydicom.read_file(test_case)\nplt.imshow(test_data.pixel_array, cmap=plt.cm.bone)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, create the brain, subdural and soft tissue windows","metadata":{}},{"cell_type":"code","source":"brain_img = window_image(test_data, 40, 80)\nsubdural_img = window_image(dcm, 80, 200)\nsoft_img = window_image(dcm, 40, 380)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Before concatenating them, we need to make sure that all of theirvalues fall under a same range so we Standardise their values","metadata":{}},{"cell_type":"code","source":"brain_img = (brain_img - 0) / 80\n# print(brain_img)\nplt.imshow(brain_img, cmap=plt.cm.bone)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(subdural_img)\nplt.imshow(subdural_img, cmap=plt.cm.bone)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subdural_img = (subdural_img - (-20))/200\n# print(subdural_img)\nplt.imshow(subdural_img, cmap=plt.cm.bone)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(soft_img)\nplt.imshow(soft_img, cmap=plt.cm.bone)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"soft_img = (soft_img - (-150))/380\n# print(soft_img)\nplt.imshow(soft_img, cmap=plt.cm.bone)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So we brought all three windows to same scale, now we will combine them into a single 3-channel image","metadata":{}},{"cell_type":"code","source":"bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1, 2, 0)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(bsb_img, cmap=plt.cm.bone)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# REFERENCES\n[1. For DICOM related info ]( https://dicomiseasy.blogspot.com/2012/08/chapter-12-pixel-data.html)\n\n[2. To understand CT scans and workflow of radiologists](https://radiopaedia.org/articles/ct-head-an-approach?lang=gb)\n\n[3. Paper referred](https://arxiv.org/pdf/2008.00302.pdf)","metadata":{}}]}