{"cells":[{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n\nfrom IPython.display import HTML","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# What is intracranial hemorrhage?\n\nHmm, let's watch a video! :-)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"HTML('<iframe width=\"800\" height=\"500\" src=\"https://www.youtube.com/embed/Kb_wzb7-rvE\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Thank you Armando! I have no medical background and this was really great to understand the topic better! :-)"},{"metadata":{},"cell_type":"markdown","source":"## Our goal\n\nWe are asked to predict the occurence and the subtype of intracranial hemorrhage. For this purpose we have to make 6 decisions per image: 5 subtypes and if there is an occurence (any)."},{"metadata":{},"cell_type":"markdown","source":"## Table of contents\n\n1. [Prepare to start](#prepare)\n2. [Exploratory analysis](#explore)\n    * [Sample submission](#sample_submission)\n    * [Evaluation metric](#evaluate)\n    * [Target distribution](#targets)\n    * [Number of samples](#num_samples)\n3. [Preprocessing dicom files](#dicom)\n    * [What is given by a dicom file?](#dicomfile)\n    * [Why can't we see something without windowing?](#aboutwindows)\n    * [What is covered by raw pixel values?](#pixelarray)\n    * [What does pixelspacing mean?](#pixelspacing)\n    * [The doctors windows](#docwindows)\n    * [My custom window](#customwindow)\n    * [The image shape](#imageshape)"},{"metadata":{},"cell_type":"markdown","source":"# Prepare to start <a class=\"anchor\" id=\"prepare\"></a>"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set()\n\nimport pydicom\n\nfrom os import listdir\n\nfrom skimage.transform import resize\nfrom imgaug import augmenters as iaa\n\nfrom sklearn.model_selection import train_test_split\n\nfrom keras.applications import ResNet50, VGG16\nfrom keras.applications.resnet50 import preprocess_input as preprocess_resnet_50\nfrom keras.applications.vgg16 import preprocess_input as preprocess_vgg_16\n\nfrom keras.utils import Sequence\n\nlistdir(\"../input/rsna-intracranial-hemorrhage-detection/\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"INPUT_PATH = \"../input/rsna-intracranial-hemorrhage-detection/\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Exploratory analysis <a class=\"anchor\" id=\"explore\"></a>\n\n## Sample Submission <a class=\"anchor\" id=\"sample_submission\"></a>"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"submission = pd.read_csv(INPUT_PATH + \"stage_1_sample_submission.csv\")\nsubmission.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can clearly see that we have to make several predictions for one image id:\n\n* epidural \n* subdural \n* subarachnoid \n* intraparenchymal \n* intraventricular \n* any - this one indicates that at least one subtype is present, hence it tells us if the patient has IH or not."},{"metadata":{},"cell_type":"markdown","source":"## Evaluation metric <a class=\"anchor\" id=\"evaluate\"></a>\n\nThe **weighted multilabel logarithmic loss** is used to score our model performance. In the description it's given that the log loss is first taken for each subtype $s$ given an image id $n$:\n\n$$l_{n,s} = t_{n,s} \\cdot \\ln(y_{n,s}) + (1-t_{n,s}) \\cdot \\ln(1-y_{n,s}) $$  \n\nThen they seem to be added, whereas the any-subtype obtains a higher weight than the others:\n\n$$ l_{n} = \\sum_{s} w_{s} \\cdot l_{n,s} $$\n\nAnd finally this loss is averaged over all samples:\n\n$$ Loss = -\\frac{1}{N} \\cdot \\sum_{n} l_{n} $$\n\nHopefully I got this right! ;-) \n\nInterestingly the competition host has not provided the weights (or did they?). Hence this is going to be an LB probing hyperparameter to keep in mind."},{"metadata":{},"cell_type":"markdown","source":"## Target distribution <a class=\"anchor\" id=\"targets\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"traindf = pd.read_csv(INPUT_PATH + \"stage_1_train.csv\")\ntraindf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"label = traindf.Label.values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"traindf = traindf.ID.str.rsplit(\"_\", n=1, expand=True)\ntraindf.loc[:, \"label\"] = label\ntraindf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"traindf = traindf.rename({0: \"id\", 1: \"subtype\"}, axis=1)\nsubtype_counts = traindf.groupby(\"subtype\").label.value_counts().unstack()\nsubtype_counts = subtype_counts.loc[:, 1] / traindf.groupby(\"subtype\").size() * 100","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"traindf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We need the same for our test data later:"},{"metadata":{"trusted":true},"cell_type":"code","source":"testdf = submission.ID.str.rsplit(\"_\", n=1, expand=True)\ntestdf = testdf.rename({0: \"id\", 1: \"subtype\"}, axis=1)\ntestdf.loc[:, \"label\"] = 0\ntestdf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"multi_target_count = traindf.groupby(\"id\").label.sum()\n\nfig, ax = plt.subplots(1,3,figsize=(20,5))\n\nsns.countplot(traindf.label, ax=ax[0], palette=\"Reds\")\nax[0].set_xlabel(\"Binary label\")\nax[0].set_title(\"How often do we observe a positive label?\");\n\nsns.countplot(multi_target_count, ax=ax[1])\nax[1].set_xlabel(\"Numer of targets per image\")\nax[1].set_ylabel(\"Frequency\")\nax[1].set_title(\"Multi-Hot occurences\")\n\nsns.barplot(x=subtype_counts.index, y=subtype_counts.values, ax=ax[2], palette=\"Set2\")\nplt.xticks(rotation=45); \nax[2].set_title(\"How much binary imbalance do we have?\")\nax[2].set_ylabel(\"% of positive occurences (1)\");\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* The first image already shows that we will have much more zero occurences than positive target values.\n* Going into details of each subtype we can see that we have to deal with high class imbalance. \n* Epidural is the worst case. For this type we only have a few (< 1%) of positive occurrences. It will be difficult to train a model that is robust enough and does not tend to overfit."},{"metadata":{},"cell_type":"markdown","source":"## Number of samples <a class=\"anchor\" id=\"num_samples\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"traindf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"traindf.id.nunique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check whether this matches the number of train images we have:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_dir = INPUT_PATH + \"stage_1_train_images/\"\ntrain_files = listdir(train_dir)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_size = len(train_files)\ntrain_size","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, as expected."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_dir = INPUT_PATH + \"stage_1_test_images/\"\ntest_files = listdir(test_dir)\ntest_size = len(test_files)\ntest_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_size/test_size","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, in this stage 1 we have almost 8.6 times more training images than the test data."},{"metadata":{},"cell_type":"markdown","source":"# Preprocessing dicom files <a class=\"anchor\" id=\"dicom\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_files[0:10]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, the filename is given by the *ID_alphanum* column. As we have used str.rsplit of pandas we are lucky and can easily load images given the id column: "},{"metadata":{"trusted":true},"cell_type":"code","source":"subtypes = traindf.subtype.unique()\nsubtypes","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## What is given by a dicom file? <a class=\"anchor\" id=\"dicomfile\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"dataset = pydicom.dcmread(train_dir + \"ID_c5c23af94.dcm\")\nprint(dataset)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Why can't we see something without windowing? <a class=\"anchor\" id=\"aboutwindows\"></a>\n\nI haven't worked often with dicom images so far, so I still get confused why we have someting like window center, width and rescale parameters. So I decided to watch a video about understanding CT images. You can speed up a bit if you like by starting at ~7min:\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/KZld-5W99cI\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, I learnt... \n\n* that hounsfield units are a measurement to describe radiodensity. \n* different tissues have different HUs.\n* Our eye can only detect ~6% change in greyscale (16 shades of grey).\n* Given 2000 HU of one image (-1000 to 1000), this means that 1 greyscale covers 8 HUs.\n* Consequently there can happen a change of 120 HUs unit our eye is able to detect an intensity change in the image.\n* The example of a hemorrhage in the brain shows relevant HUs in the range of 8-70. We won't be able to see important changes in the intensity to detect the hemorrhage. \n* This is the reason why we have to focus 256 shades of grey into a small range/window of HU units. (WINDOW)\n* The level means where this window is centered.\n"},{"metadata":{},"cell_type":"markdown","source":"## What is covered by raw pixel values? <a class=\"anchor\" id=\"pixelarray\"></a>\n\nNo! If we browse through the dicom files, we can see that this is not true. Let's see how different dicom datasets differ in the distribution of pixel array values:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,1,figsize=(20,10))\nfor file in train_files[0:10]:\n    dataset = pydicom.dcmread(train_dir + file)\n    image = dataset.pixel_array.flatten()\n    rescaled_image = image * dataset.RescaleSlope + dataset.RescaleIntercept\n    sns.distplot(image.flatten(), ax=ax[0]);\n    sns.distplot(rescaled_image.flatten(), ax=ax[1])\nax[0].set_title(\"Raw pixel array distributions for 10 examples\")\nax[1].set_title(\"HU unit distributions for 10 examples\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Crazy, that we observe cases with -1000 and -2000. What does that mean? We need to find it out!\n* Nonetheless the mode is located at 0 for all example raw pixel value arrays. This is likely to correspond to air. Consequently our raw pixel values are not given as HU units. \n* But we can transform the image to HU units by scaling with the slope and intercept.\n* Then we can see that mode is located at -1000 (HU of air).\n\nTo understand this part, the following tutorial was extremely helpful for me: https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial. Thank you a lot @Guido Zuidhof. Let's try to understand the -2000 and -3000 cases of rescaled images:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(4,10,figsize=(20,12))\n\nfor n in range(10):\n    dataset = pydicom.dcmread(train_dir + train_files[n])\n    image = dataset.pixel_array\n    rescaled_image = image * dataset.RescaleSlope + dataset.RescaleIntercept\n    mask2000 = np.where((rescaled_image <= -1500) & (rescaled_image > -2500), 1, 0)\n    mask3000 = np.where(rescaled_image <= -2500, 1, 0)\n    ax[0,n].imshow(rescaled_image)\n    rescaled_image[rescaled_image < -1024] = -1024\n    ax[1,n].imshow(mask2000)\n    ax[2,n].imshow(mask3000)\n    ax[3,n].imshow(rescaled_image)\n    ax[0,n].grid(False)\n    ax[1,n].grid(False)\n    ax[2,n].grid(False)\n    ax[3,n].grid(False)\nax[0,0].set_title(\"Rescaled image\")\nax[1,0].set_title(\"Mask -2000\")\nax[2,0].set_title(\"Mask -3000\");\nax[3,0].set_title(\"Background to air\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Yeah! Great! We can see that both cases (-2000 and -3000) correspond to the outside region ot cylindrical CT scanners. Consequently we can set these values to -1000 (air in HU) without worries."},{"metadata":{},"cell_type":"markdown","source":"## What does pixel spacing mean? <a class=\"anchor\" id=\"pixelspacing\"></a>\n\nWhen browsing through the dicom files, we can see that a value called pixel spacing changes as well! I don't know what that means, but perhaps we can understand it by looking at some extremes."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"pixelspacing_w = []\npixelspacing_h = []\nspacing_filenames = []\nfor file in train_files[0:1000]:\n    dataset = pydicom.dcmread(train_dir + file)\n    spacing = dataset.PixelSpacing\n    pixelspacing_w.append(spacing[0])\n    pixelspacing_h.append(spacing[1])\n    spacing_filenames.append(file)\n\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(pixelspacing_w, ax=ax[0], color=\"Limegreen\", kde=False)\nax[0].set_title(\"Pixel spacing width \\n distribution\")\nax[0].set_ylabel(\"Frequency given 1000 images\")\nsns.distplot(pixelspacing_h, ax=ax[1], color=\"Mediumseagreen\", kde=False)\nax[1].set_title(\"Pixel spacing height \\n distribution\");\nax[1].set_ylabel(\"Frequency given 1000 images\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"min_file = spacing_filenames[np.argmin(pixelspacing_w)]\nmax_file = spacing_filenames[np.argmax(pixelspacing_w)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def rescale_pixelarray(dataset):\n    image = dataset.pixel_array\n    rescaled_image = image * dataset.RescaleSlope + dataset.RescaleIntercept\n    rescaled_image[rescaled_image < -1024] = -1024\n    return rescaled_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,10))\n\ndataset_min = pydicom.dcmread(train_dir + min_file)\nimage_min = rescale_pixelarray(dataset_min)\n\ndataset_max = pydicom.dcmread(train_dir + max_file)\nimage_max = rescale_pixelarray(dataset_max)\n\nax[0].imshow(image_min, cmap=\"Spectral\")\nax[0].set_title(\"Pixel spacing w: \" + str(np.min(pixelspacing_w)))\nax[1].imshow(image_max, cmap=\"Spectral\");\nax[1].set_title(\"Pixel spacing w: \" + str(np.max(pixelspacing_w)))\nax[0].grid(False)\nax[1].grid(False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Hmm... still don't know what this means. Reading in the [docs](https://dicom.innolitics.com/ciods/ct-image/image-plane/00280030):\n\n*All pixel spacing related Attributes are encoded as the physical distance between the centers of each two-dimensional pixel, specified by two numeric values.The first value is the row spacing in mm, that is the spacing between the centers of adjacent rows, or vertical spacing.The second value is the column spacing in mm, that is the spacing between the centers of adjacent columns, or horizontal spacing.*\n\nConsequently it's related to the physical distance. Let's understand it given our extreme examples. The pixel spacing yields the mm of physical distance of one pixel. Hence we can compute the overall distance covered by one image width or height:"},{"metadata":{"trusted":true},"cell_type":"code","source":"np.min(pixelspacing_w) * 512","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is the true length in mm covered by our image. In my current case (can change as listdir choses file order at random) it's 195 mm meaning 19.5 cm. "},{"metadata":{"trusted":true},"cell_type":"code","source":"np.max(pixelspacing_w) * 512","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And in the maximum case it's ~500 mm, consequently 50 cm. \n\nWhat does this mean for us? It means that small heads do not automatically mean that this is a child. It could also be an adult but zoomed out. I'm expected that the dataset holds patients with varying true head sizes ranging from childrens to adults. But the resolution might differ from scan to scan even if patients show same head sizes. Nonetheless our model should still be able to detect the different types of hemorrhage. I would say it's worth to add zooming as image augmentation technique to our workflow later. What do you think? "},{"metadata":{},"cell_type":"markdown","source":"## The doctors windows <a class=\"anchor\" id=\"docwindows\"></a>\n\nTaking a look at the dicom dataset again we can see that there is already a window center and width given for us. But what does that mean? Was this done by a doctor who set the range to visualise the hemorrhage? Is this important for our algorithm? **We should be very careful now. If we would simply put 256 shades of grey into one window this would differ from patient to patient as the given window ranges are different. Consequently we would introduce a source of variation that is not given by original HU units per image.** What to do instead?\n\nI would like to collect window centers and width of ~1000 images to see the varity of doctos favorite windows. Then I would like to **setup a fixed window level and width that covers the majority of all window properties**. This way we can compare if a fixed custom window size is better suited that individual doctor window sizes. \n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_window_value(feature):\n    if type(feature) == pydicom.multival.MultiValue:\n        return np.int(feature[0])\n    else:\n        return np.int(feature)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"window_widths = []\nwindow_levels = []\nspacing_filenames = []\nfor file in train_files[0:1000]:\n    dataset = pydicom.dcmread(train_dir + file)\n    win_width = get_window_value(dataset.WindowWidth)\n    win_center = get_window_value(dataset.WindowCenter)\n    window_widths.append(win_width)\n    window_levels.append(win_center)\n    spacing_filenames.append(file)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,2,figsize=(20,15))\n\nsns.distplot(window_widths, kde=False, ax=ax[0,0], color=\"Tomato\")\nax[0,0].set_title(\"Window width distribution \\n of 1000 images\")\nax[0,0].set_xlabel(\"Window width\")\nax[0,0].set_ylabel(\"Frequency\")\n\nsns.distplot(window_levels, kde=False, ax=ax[0,1], color=\"Firebrick\")\nax[0,1].set_title(\"Window level distribution \\n of 1000 images\")\nax[0,1].set_xlabel(\"Window level\")\nax[0,1].set_ylabel(\"Frequency\")\n\nsns.distplot(np.log(window_widths), kde=False, ax=ax[1,0], color=\"Tomato\")\nax[1,0].set_title(\"Log window width distribution \\n of 1000 images\")\nax[1,0].set_xlabel(\"Log window width\")\nax[1,0].set_ylabel(\"Frequency\")\n\nsns.distplot(np.log(window_levels), kde=False, ax=ax[1,1], color=\"Firebrick\")\nax[1,1].set_title(\"Log window level distribution \\n of 1000 images\")\nax[1,1].set_xlabel(\"Log window level\")\nax[1,1].set_ylabel(\"Frequency\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Uhh! Very bad! Do you see it? There are extreme outliers in the window widths and centers. This does not make sense! I would have expected the doctos to focus on the brain tissue and we have already learnt that HU values are roughly between 8-70 in these cases. \n* Let's explore if we are better by focusing on our own window size! ;-)"},{"metadata":{},"cell_type":"markdown","source":"## My custom window <a class=\"anchor\" id=\"customwindow\"></a>\n\nOk, let's setup our own window width and center:"},{"metadata":{"trusted":true},"cell_type":"code","source":"doc_windows = pd.DataFrame(index=spacing_filenames, columns=[\"win_width\", \"win_level\"])\ndoc_windows[\"win_width\"] = window_widths\ndoc_windows[\"win_level\"] = window_levels\ndoc_windows.head(20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"doc_windows[doc_windows.win_width==doc_windows.win_width.median()]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"doc_windows.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.quantile(window_widths, 0.95)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.quantile(window_levels, 0.95)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* The majority of centers is located between 30 and 40 HUs.\n* In contrast the widths of windows vary between 70 to 150 HUs."},{"metadata":{},"cell_type":"markdown","source":"Now, let's setup a window that is centered at 30 and has a width of 80. This would be close to the median for both cases. Furthermore we have seen such combinations of width and center in our doc_windows dataframe."},{"metadata":{"trusted":true},"cell_type":"code","source":"def set_manual_window(hu_image, custom_center, custom_width):\n    min_value = custom_center - (custom_width/2)\n    max_value = custom_center + (custom_width/2)\n    hu_image[hu_image < min_value] = min_value\n    hu_image[hu_image > max_value] = max_value\n    return hu_image","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(3,4,figsize=(20,15))\n\ndocs_dict = {\"ID_352b300f9.dcm\": {\"width\": 4000, \"level\": 600},\n             \"ID_7e7d7633a.dcm\": {\"width\": 70, \"level\": 30},\n             \"ID_87e8b2528.dcm\": {\"width\": 80, \"level\": 40}}\nn = 0\nfor file in [\"ID_352b300f9.dcm\", \"ID_7e7d7633a.dcm\", \"ID_87e8b2528.dcm\"]:\n    dataset = pydicom.dcmread(train_dir + file)\n    pixelarray = dataset.pixel_array\n    ax[n,0].imshow(pixelarray, cmap=\"Spectral\")\n    ax[n,0].grid(False)\n    rescaled_image = rescale_pixelarray(dataset)\n    ax[n,1].imshow(rescaled_image, cmap=\"Spectral\")\n    ax[n,1].grid(False)\n    \n    org_windowed_image = set_manual_window(rescaled_image, docs_dict[file][\"level\"], docs_dict[file][\"width\"])\n    ax[n,2].imshow(org_windowed_image, cmap=\"Spectral\")\n    ax[n,2].grid(False)\n    \n    new_windowed_image = set_manual_window(rescaled_image, 40, 150)\n    ax[n,3].imshow(new_windowed_image, cmap=\"Spectral\")\n    ax[n,3].grid(False)\n    \n    n+=1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* The first extreme case seems to be faulty. After windowing to a center of 40 and width of 150 we can't see the same nice patterns as for the min and the median cases.\n* By windowing to the same window we can't see differences between the median and the min case. But pixel values differ! "},{"metadata":{},"cell_type":"markdown","source":"## The image shape <a class=\"anchor\" id=\"imageshape\"></a>"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"num_rows = []\nnum_cols = []\nspacing_filenames = []\nfor file in train_files[0:1000]:\n    dataset = pydicom.dcmread(train_dir + file)\n    num_rows.append(dataset.Rows)\n    num_cols.append(dataset.Columns)\n    spacing_filenames.append(file)\n\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(num_rows, ax=ax[0], color=\"Purple\", kde=False)\nax[0].set_title(\"Number of rows \\n distribution\")\nax[0].set_ylabel(\"Frequency given 1000 train images\")\nsns.distplot(num_cols, ax=ax[1], color=\"Violet\", kde=False)\nax[1].set_title(\"Number of columns \\n distribution\");\nax[1].set_ylabel(\"Frequency given 1000 train images\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok it seems that most images are of shape 512x512. But I don't know if we can be sure about it. Is it possible that this can change from stage 1 to stage 2? What about the test data?"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"num_rows = []\nnum_cols = []\nspacing_filenames = []\nfor file in test_files[0:1000]:\n    dataset = pydicom.dcmread(test_dir + file)\n    num_rows.append(dataset.Rows)\n    num_cols.append(dataset.Columns)\n    spacing_filenames.append(file)\n\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(num_rows, ax=ax[0], color=\"Purple\", kde=False)\nax[0].set_title(\"Number of rows \\n distribution\")\nax[0].set_ylabel(\"Frequency given 1000 test images\")\nsns.distplot(num_cols, ax=ax[1], color=\"Violet\", kde=False)\nax[1].set_title(\"Number of columns \\n distribution\");\nax[1].set_ylabel(\"Frequency given 1000 test images\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok for stage 1 this holds as well."},{"metadata":{},"cell_type":"markdown","source":"## Where to go next?\n\nThe kernel had been much longer in the past, but .... I reached the end of the kaggle kernel world by exceeding the maximum kernel length that can be displayed. Jackpot - This was the first time for me! :-) \n\n<img src=\"https://cdn.pixabay.com/photo/2017/05/26/11/37/worlds-end-2345861_1280.jpg\" width=\"600px\">\n\nI'm sorry that this happend. Have to learn how to keep complex stuff simple, more concise and light-weighted. If you like to continue with the modelling part - you can find it here:\n\nhttps://www.kaggle.com/allunia/rsna-ih-detection-baseline\n\nHappy kaggling!"}],"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":4,"nbformat_minor":1}