{"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":"# Introduction to DICOM & Voxels\n\n\n### What is a DICOM?\n**D**igital **I**maging and **Co**mmunications in **M**edicine (DICOM) - an international standard related to the exchange, storage and communication of digital medical images. Prior to this format, there was no standardized way to transfer medical scans. So loading up a single patient's study outside the hospital, in older formats took about 10-30 minutes for a single scan! \n\nWhile DICOM 16-bit images (with values ranging from -32768..32767), other 8-bit greyscale images store values 0 - 255. These value ranges in DICOM are useful, as they correlate with the [Hounsfield Scale](https://en.wikipedia.org/wiki/Hounsfield_scale). Each voxel can store a large amount of information.\n\nNB: If you want to get right to the code & image example, click [here](#example)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true}},{"cell_type":"markdown","source":"### Physics of CT Scans\n**CT (computed tomography)** uses X-rays to obtain images. A heated cathode releases high-energy electrons, which in turn release their energy as X-ray radiation. X-rays pass through tissues and hit a detector on the other side.**The more dense a tissue, the more X-rays it absorbs.**\n\n* Bone: X-rays absorbed = few X-rays reaching detector: White\n* Air: X-rays not absorbed = lots of X-rays reaching detector: Black\nCompared to plain film, **CT is able to distinguish more subtle density differences and there is no overlap of structures.\n\nCurrent CT machines use ‘Spiral CT’. This consists of a single radiation source with multiple detectors which rotates around the patient, obtaining a block of data as the patient is moved through.\n\n\"\"\n\nThe information obtained can be reconstructed by a computer to form a 3D “volume”, which can then be “re-sliced” digitally to obtain thinner slices as well as slices in different planes.\n\n### Phases of CT Scan:\nPhases of a scan refer to when the images are taken, relative to time of contrast administration.\n**Contrast Administration:** Intravenous contrast is used in CT to help highlight blood vessels and to enhance the tissue structure of various organs such as the brain, spine, liver and kidneys. \"Intravenous\" means that the contrast is injected into a vein using a small needle. Some imaging exams of the abdomen and gastrointestinal system use both the intravenous iodine and orally administered barium contrast for maximum sensitivity.\n\n$$(non or pre contrast) > (arterial) > (venous) > (delayed)$$\n\nThe arterial phase comes before the venous phase, because even though contrast is given into a vein, within approximately 30 seconds, the contrast has passed through the heart and into the arterial system.\n\n* The chest is usually imaged in the arterial phase.\n* The abdomen is usually imaged in the (portal) venous phase.\n* Liver lesions are usually imaged with a multiple phase scan.\n\nIn a CT chest/abdo/pelvis, the lung bases/liver may be imaged twice (overlap between chest - arterial phase, and abdomen - venous phase)","metadata":{}},{"cell_type":"markdown","source":"### Windowing: \nWindowing, also known as **grey-level mapping**, **contrast stretching**, **histogram modification** 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. The brightness of the image is adjusted via the window level. The contrast is adjusted via the window width.\n\nTissue density is measured in [Hounsfield units (HU)](https://en.wikipedia.org/wiki/Hounsfield_scale)\n\n* This is defined as **Air = −1000 HU**; **Water = 0 HU.**\n\nDensity of tissues in CT-Scans: $$Air < Fat < Fluid < Soft tissue < Bone < Metal$$\nThe easier way to remember this is (Fat floats on water, so is less dense than fluid; Soft tissue is mostly intracellular fluid with some connective tissue)\n\n* Air = −1000 HU\n* Lung ≈ −500 HU (partially air, partially soft tissue)\n* Fat ≈ −50 HU (slightly less dense than simple fluid)\n* Water = 0 HU\n* Soft tissue (& blood) ≈ +50 HU (slightly more dense than simple fluid)\n* Bone ≈ +1000 HU (much more dense)\n\n\nTo ascertain a window, a ‘level’ and a ‘width’ is defined. For example, a window with a level of 0 HU and a width of 400 HU will have a range of −200 HU to +200 HU. Any tissue with a density of −200 HU or less will be black, and any tissue with a density of +200 HU or more will be white. And values between -200 HU to +200 HU will be spread between the whole grayscale range. A **window** can be set to look at certain tissues of interest. **A small range of tissue density is represented by a full greyscale spectrum from black to white, thus making subtle density differences within the specified range easier to see.**\n\n\n\n\n\n\n### Typical window width and level values: \nAlthough this varies somewhat from institution to institution and vendor to vendor, window width and centers are generally fairly similar. **The values below are written as width and level (W:x L:y) in Hounsfield units (HU).**\n\n* head and neck\n* brain W:80 L:40\n* subdural W:130-300 L:50-100\n* stroke W:8 L:32 or W:40 L:40 3\n* temporal bones W:2800 L:600\n* soft tissues: W:350–400 L:20–60 4\n* chest\n* lungs W:1500 L:-600\n* mediastinum W:350 L:50\n* abdomen\n* soft tissues W:400 L:50\n* liver W:150 L:30\n* spine\n* soft tissues W:250 L:50\n* bone W:1800 L:400\n","metadata":{}},{"cell_type":"markdown","source":"\n### Window width\nThe window width (WW) as the name suggests is the measure of the range of CT numbers that an image contains. A wider window width (2000 HU), therefore, will display a wider range of CT numbers. Consequently, the transition of dark to light structures will occur over a larger transition area to that of a narrow window width (<1000 HU). Accordingly, it is important to note, that a significantly wide window displaying all the CT numbers will result in different attenuations between soft tissues to become obscured.\n\n**Wide window:** When you are looking at an area with predominantly different tissue density, a wide window is used. A good example is lungs or cortical tissue, where air and vessels will sit side by side.\n\n**Narrow window:** When you are looking at tissues with almost similar density, you should use narrow window. As a result subtle changes in tissued density (small window) is magnified over the whole grayscale range. \n\n### Window level/center\nThe window level (WL), often also referred to as window center, is the midpoint of the range of the CT numbers displayed. **When the window level is decreased the CT image will be brighter and vice versa.** ","metadata":{}},{"cell_type":"markdown","source":"### Upper and lower grey level calculation \nWhen presented with a Window width (WW) and Widnow Level (WL) one can calculate the upper and lower grey levels i.e. values over x will be white and values below y will be black. \n\n* the upper grey level (x) is calculated via WL + (WW ÷ 2)\n* the lower grey level (y) is calculated via WL - (WW ÷ 2)\n\nFor example, a brain is W:80 L:40. Therefore, all values above +80 will be set to maximum grayscale level (white) and all values below 0 will be set to lowest grayscale level in the display (black). And values between +0 to +80 will be spread between the whole grayscale range. \n\nExamples of commonly used windows are soft tissue, lung, and bone are given below: ","metadata":{}},{"cell_type":"markdown","source":"### Lung tissue window\nA lung window is used to view lung parenchyma. Lung parenchyma (−500 HU) would be within range, appearing grey. Air pockets (−1000 HU) around the lung, such as pneumothorax or bullae, would appear black, thus allowing clear differentiation.\n![](\"https://www.radiologycafe.com/images/basics/window-lung.png\")
\nLevel: −200 HU; Width: 2000 HU (Range: −1200 to +800),\n![](\"https://www.radiologycafe.com/images/basics/ct-window-lung.jpg\")
\nLung window on a chest CT","metadata":{}},{"cell_type":"markdown","source":"### Soft tissue window\nA soft tissue window is used to view most organs. A soft tissue window cannot be used for lung parenchyma, as lung density (−500 HU) is outside range and will appear completely black.\n![](\"https://www.radiologycafe.com/images/basics/window-soft.png\")
\n\nLevel: +50 HU; Width: 350 HU (Range: −125 to +225)\n![](\"https://www.radiologycafe.com/images/basics/ct-window-soft.jpg\")
\nSoft tissue window on a chest CT","metadata":{}},{"cell_type":"markdown","source":"### Bone window\nA bone window is used to view bone detail.\nThere is good differentiation within the range of high densities found in the bony cortex and medulla.\n![](\"https://www.radiologycafe.com/images/basics/window-bone.png\")
\nLevel: 300 HU; Width: 2000 HU (Range: −700 to +1300)\n![](\"https://www.radiologycafe.com/images/basics/ct-window-bone.jpg\")
\nBone window on a chest CT","metadata":{}},{"cell_type":"markdown","source":"\n# Workflow of a Radiologist\n\nThis website shows a typical workflow, thought process and professional approach that a radiologists takes, when given a task to detect any abnormalities on a CT scan of the lung. [radiopedia.org](https://radiopaedia.org/articles/ct-chest-non-contrast-technique?lang=gb) There are some details that I will skip over, as they are dealing with 3D scans in real life. In our dataset, for each patient ID we have only one slice (2D) of the lung to make our diagnosis from. \n\n> To be as good as a radiologist, you have to start thinking like one.\n\nThe issue that I want to highlight is how important windows are in a **radiologist's workflow**. Just looking at the bare images provided in the raw CT-scans provides a good looking image for non-med people however, from a radiologist's point of view, that provides little information about what he is actually looking for. In this case the importance of windowing comes in. \n","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0"}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd \nimport matplotlib.pyplot as plt\nimport os, glob\nimport pydicom as dcm\nprint('Loaded in libraries!')","metadata":{"execution":{"iopub.status.busy":"2022-06-13T07:53:27.066063Z","iopub.execute_input":"2022-06-13T07:53:27.066591Z","iopub.status.idle":"2022-06-13T07:53:27.293108Z","shell.execute_reply.started":"2022-06-13T07:53:27.066555Z","shell.execute_reply":"2022-06-13T07:53:27.292128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PATH = \"../input/rsna-str-pulmonary-embolism-detection/\"\n\ntrain_df = pd.read_csv(PATH + \"train.csv\")\ntest_df = pd.read_csv(PATH + \"test.csv\")\n\nTRAIN_PATH = PATH + \"train/\"\nTEST_PATH = PATH + \"test/\"\nsub = pd.read_csv(PATH + \"sample_submission.csv\")\ntrain_image_file_paths = glob.glob(TRAIN_PATH + '/*/*/*.dcm')\ntest_image_file_paths = glob.glob(TEST_PATH + '/*/*/*.dcm')\n\nprint(f'Train dataframe shape :{train_df.shape}')\nprint(f'Test dataframe shape :{test_df.shape}')\n\nprint(f'Number of train images : {len(train_image_file_paths)}')\nprint(f'Number of test images : {len(test_image_file_paths)}')","metadata":{"execution":{"iopub.status.busy":"2022-06-13T07:53:27.29497Z","iopub.execute_input":"2022-06-13T07:53:27.295289Z","iopub.status.idle":"2022-06-13T08:26:07.893315Z","shell.execute_reply.started":"2022-06-13T07:53:27.295255Z","shell.execute_reply":"2022-06-13T08:26:07.892145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Utility Functions. ","metadata":{}},{"cell_type":"code","source":"# Function to take care of teh translation and windowing. \ndef window_image(img, window_center,window_width, intercept, slope, rescale=True):\n img = (img*slope +intercept) #for translation adjustments given in the dicom file. \n img_min = window_center - window_width//2 #minimum HU level\n img_max = window_center + window_width//2 #maximum HU level\n img[imgimg_max] = img_max #set img_max for all HU levels higher than maximum HU level\n if rescale: \n img = (img - img_min) / (img_max - img_min)*255.0 \n return img\n \ndef get_first_of_dicom_field_as_int(x):\n #get x[0] as in int is x is a 'pydicom.multival.MultiValue', otherwise get int(x)\n if type(x) == dcm.multival.MultiValue: return int(x[0])\n else: return int(x)\n \ndef get_windowing(data):\n dicom_fields = [data[('0028','1050')].value, #window center\n data[('0028','1051')].value, #window width\n data[('0028','1052')].value, #intercept\n data[('0028','1053')].value] #slope\n return [get_first_of_dicom_field_as_int(x) for x in dicom_fields]","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:07.895047Z","iopub.execute_input":"2022-06-13T08:26:07.895457Z","iopub.status.idle":"2022-06-13T08:26:07.908683Z","shell.execute_reply.started":"2022-06-13T08:26:07.895412Z","shell.execute_reply":"2022-06-13T08:26:07.907695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def view_images(files, title = '', aug = None, windowing = True):\n width = 2\n height = 2\n fig, axs = plt.subplots(height, width, figsize=(15,15))\n \n for im in range(0, height * width):\n data = dcm.dcmread(files[im])\n image = data.pixel_array\n window_center , window_width, intercept, slope = get_windowing(data)\n if windowing:\n output = window_image(image, window_center, window_width, intercept, slope, rescale = False)\n else:\n output = image\n i = im // width\n j = im % width\n axs[i,j].imshow(output, cmap=plt.cm.gray) \n axs[i,j].axis('off')\n \n plt.suptitle(title)\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:07.910555Z","iopub.execute_input":"2022-06-13T08:26:07.911Z","iopub.status.idle":"2022-06-13T08:26:07.926688Z","shell.execute_reply.started":"2022-06-13T08:26:07.910955Z","shell.execute_reply":"2022-06-13T08:26:07.925504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_images(train_image_file_paths[3200:], 'Images with Windowing')","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:07.930318Z","iopub.execute_input":"2022-06-13T08:26:07.930652Z","iopub.status.idle":"2022-06-13T08:26:09.325834Z","shell.execute_reply.started":"2022-06-13T08:26:07.930618Z","shell.execute_reply":"2022-06-13T08:26:09.325053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_images(train_image_file_paths[3200:], title = 'Images with Windowing', windowing=False)","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:09.327109Z","iopub.execute_input":"2022-06-13T08:26:09.327388Z","iopub.status.idle":"2022-06-13T08:26:09.912887Z","shell.execute_reply.started":"2022-06-13T08:26:09.327358Z","shell.execute_reply":"2022-06-13T08:26:09.911959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = dcm.dcmread(train_image_file_paths[3203])\nimage = data.pixel_array\nwindow_center , window_width, intercept, slope = get_windowing(data)\noutput = window_image(image, window_center, window_width, intercept, slope, rescale = False)\nf, axarr = plt.subplots(1,2, figsize=(15,10))\naxarr[0].imshow(image, cmap='gray')\naxarr[1].imshow(output, cmap = 'gray')","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:09.914136Z","iopub.execute_input":"2022-06-13T08:26:09.914416Z","iopub.status.idle":"2022-06-13T08:26:10.275817Z","shell.execute_reply.started":"2022-06-13T08:26:09.914386Z","shell.execute_reply":"2022-06-13T08:26:10.27462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = dcm.dcmread(train_image_file_paths[3203])\nimage = data.pixel_array\nwindow_center , window_width, intercept, slope = get_windowing(data)\n\nprint(window_center , window_width, intercept, slope)\nfrom ipywidgets import interact\n","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:10.277306Z","iopub.execute_input":"2022-06-13T08:26:10.277636Z","iopub.status.idle":"2022-06-13T08:26:10.28778Z","shell.execute_reply.started":"2022-06-13T08:26:10.277603Z","shell.execute_reply":"2022-06-13T08:26:10.286968Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def int_print(window_center , window_width=500, intercept=-1024, slope=1):\n output = window_image(image, window_center, window_width, intercept, slope, rescale = False)\n f, axarr = plt.subplots(1,2, figsize=(15,10))\n axarr[0].imshow(image, cmap='gray')\n axarr[1].imshow(output, cmap = 'gray')\n ","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:10.288758Z","iopub.execute_input":"2022-06-13T08:26:10.28908Z","iopub.status.idle":"2022-06-13T08:26:10.297173Z","shell.execute_reply.started":"2022-06-13T08:26:10.289048Z","shell.execute_reply":"2022-06-13T08:26:10.295985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Interactively set the window parameter for the best view:\nThe following segment of code produces an interactive window where you are able to change the parameter of the function. \nFor viewing interactive features, please fork the code and then run it separately to view the interactive features. ","metadata":{}},{"cell_type":"code","source":"interact(int_print, window_center= 1000)","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:10.298023Z","iopub.execute_input":"2022-06-13T08:26:10.298245Z","iopub.status.idle":"2022-06-13T08:26:10.660522Z","shell.execute_reply.started":"2022-06-13T08:26:10.298221Z","shell.execute_reply":"2022-06-13T08:26:10.65956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\n\nPATH = \"../input/rsna-str-pulmonary-embolism-detection/\"\ntrain = pd.read_csv(PATH + \"train.csv\")\nsub = pd.read_csv(PATH + \"sample_submission.csv\")\n\nfeats = list(train.columns[3:5])+list(train.columns[8:12])+list(train.columns[13:17])\nmeans = train[feats].mean().to_dict()\n\n\nsub['label'] = 0.28\nfor feat in means.keys():\n sub.loc[sub.id.str.contains(feat, regex=False), 'label'] = means[feat]\n \nsub.to_csv('submission.csv', index = False)","metadata":{"execution":{"iopub.status.busy":"2022-06-13T08:26:10.661764Z","iopub.execute_input":"2022-06-13T08:26:10.662123Z","iopub.status.idle":"2022-06-13T08:26:16.532264Z","shell.execute_reply.started":"2022-06-13T08:26:10.662089Z","shell.execute_reply":"2022-06-13T08:26:16.531198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## References.\n* [See like a Radiologist with Systematic Windowing](https://www.kaggle.com/dcstang/see-like-a-radiologist-with-systematic-windowing) by [David Tang](https://www.kaggle.com/dcstang)\n* https://radiopaedia.org/articles/ct-chest-non-contrast-technique?lang=gb\n* https://radiopaedia.org/articles/ct-pulmonary-angiogram-technique?lang=gb\n* https://en.wikipedia.org/wiki/Hounsfield_scale\n* [EDA: View dicom images with correct windowing](https://www.kaggle.com/omission/eda-view-dicom-images-with-correct-windowing) by [Richard McKinley](https://www.kaggle.com/omission)\n","metadata":{}},{"cell_type":"markdown","source":"\n\n\n\n\n\n### In the meantime, check out my other ongoing works in this same competition: \n💥 [RSNA-STR Pulmonary Embolism [Dummy Sub]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-dummy-sub)
\n💥 [CT-Scans, DICOM files, Windowing Explained](https://www.kaggle.com/redwankarimsony/ct-scans-dicom-files-windowing-explained)
\n💥 [RSNA-STR-PE [Gradient & Sigmoid Windowing]](https://www.kaggle.com/redwankarimsony/rsna-str-pe-gradient-sigmoid-windowing)
\n💥 [RSNA-STR [✔️3D Stacking ✔️3D Plot ✔️Segmentation]](https://www.kaggle.com/redwankarimsony/rsna-str-3d-stacking-3d-plot-segmentation/edit/run/42517982)
\n💥 [RSNA-STR [DICOM 👉 GIF 👉 npy]](https://www.kaggle.com/redwankarimsony/rsna-str-dicom-gif-npy)
\n💥 [RSNA-STR Pulmonary Embolism [EDA]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-eda)
\n\n","metadata":{"trusted":true}}]}
