{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Gradient & Sigmoid Windowing\n\nI've been exploring a bunch of different ways to window the DICOM images and I thought I'd share a few ideas and results.\n\nHuge thanks to [David Tang](https://www.kaggle.com/dcstang/see-like-a-radiologist-with-systematic-windowing), [Marco](https://www.kaggle.com/marcovasquez/basic-eda-data-visualization), [Nanashi](https://www.kaggle.com/jesucristo/rsna-introduction-eda-models), and [Richard McKinley](https://www.kaggle.com/omission/eda-view-dicom-images-with-correct-windowing) and for their amazing Kernels that I totally borrowed code and ideas from.\n\n**Contents**\n* [No Windowing](#1)\n* [Brain Windowing](#2)\n* [Metadata Windowing](#3)\n* [One Window, Three Channels](#4)\n* [Gradient Windowing](#5)\n* [Brain + Subdural + Bone Windowing](#6)\n* [Exclusive Windowing](#7)\n* [Gradient (Brain + Subdural + Bone) Windowing](#8)\n* [Sigmoid Windowing](#9)\n* [Sigmoid (Brain + Subdural + Bone) Windowing](#10)\n* [Sigmoid Gradient (Brain + Subdural + Bone) Windowing](#11)\n* [Acknowledgements](#12)"},{"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\nimport pydicom\nimport os\n\nTRAIN_IMG_PATH = \"../input/rsna-intracranial-hemorrhage-detection/stage_1_train_images/\"\nTEST_IMG_PATH = \"../input/rsna-intracranial-hemorrhage-detection/stage_1_test_images/\"\nBASE_PATH = '/kaggle/input/rsna-intracranial-hemorrhage-detection/'\nTRAIN_DIR = 'stage_1_train_images/'\n\ntrain = pd.read_csv(\"../input/rsna-intracranial-hemorrhage-detection/stage_1_train.csv\")\ntrain_images = os.listdir(\"../input/rsna-intracranial-hemorrhage-detection/stage_1_train_images/\")\n\ntrain['filename'] = train['ID'].apply(lambda st: \"ID_\" + st.split('_')[1] + \".dcm\")\ntrain['type'] = train['ID'].apply(lambda st: st.split('_')[2])\ntrain = train[['Label', 'filename', 'type']].drop_duplicates().pivot(index='filename', columns='type', values='Label').reset_index()\n\nhem_types = ['epidural', 'intraparenchymal', 'intraventricular', 'subarachnoid', 'subdural']\n\n\ndef load_random_images():\n    image_names = [list(train[train[h_type] == 1].sample(1)['filename'])[0] for h_type in hem_types]\n    image_names += list(train[train['any'] == 0].sample(5)['filename'])\n    return [pydicom.read_file(os.path.join(TRAIN_IMG_PATH, img_name)) for img_name in image_names]\n\n\ndef view_images(images):\n    width = 5\n    height = 2\n    fig, axs = plt.subplots(height, width, figsize=(15,5))\n    \n    for im in range(0, height * width):\n        image = images[im]\n        i = im // width\n        j = im % width\n        axs[i,j].imshow(image, cmap=plt.cm.bone) \n        axs[i,j].axis('off')\n        title = hem_types[im] if im < len(hem_types) else 'normal'\n        axs[i,j].set_title(title)\n\n    plt.show()\n    \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) == pydicom.multival.MultiValue:\n        return int(x[0])\n    else:\n        return int(x)\n\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]\n\nprint('Loaded packages and setup utility functions')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's load a some images. We'll randomly pick one positive example from each class and 5 negative (normal) examples. "},{"metadata":{"trusted":true},"cell_type":"code","source":"imgs = load_random_images()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('DICOM Pixel Data Shape: ', imgs[0].pixel_array.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('Distribution of DICOM Pixel Values')\nax = plt.hist(np.array([img.pixel_array for img in imgs]).flatten(), bins=50, color='c')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see the DICOM data is 2-dimensional and has a range of values much wider than the typical png or jpg image. Let's remind ourselves what the scans look like if we include the full range of values... \n\n### No Windowing <a></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"view_images([img.pixel_array for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These don't look very useful for detecting hemorrhages.\n\n[David Tang's Kernel](https://www.kaggle.com/dcstang/see-like-a-radiologist-with-systematic-windowing) and [Richard McKinley's Kernel](https://www.kaggle.com/omission/eda-view-dicom-images-with-correct-windowing) both do a great job at explaining windowing and how it's used by radiologists. (Make sure to check them out!)\n\nThe way I've been thinking about it, is **how do we transform our 2D scan data into 3D image data in a way that makes it easy for our model to detect hemorrhages?**\n"},{"metadata":{},"cell_type":"markdown","source":"### Brain Windowing <a></a>\nLet's check out the range of values that corresponds with brain matter. We'll clip everything outside that range so that there's more contrast in the brain-range."},{"metadata":{"trusted":true},"cell_type":"code","source":"def brain_window(img):\n    window_min = 0\n    window_max = 80\n    _, _, intercept, slope = get_windowing(img)\n    img = img.pixel_array\n    img = img * slope + intercept\n    img[img < window_min] = window_min\n    img[img > window_max] = window_max\n    img = (img - np.min(img)) / (np.max(img) - np.min(img))\n    return img\n\nview_images([brain_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Metadata Windowing <a></a>\nThe DICOM images come with metadata specifying a window center and width. We could also use these values instead of the fixed range from above."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"def metadata_window(img, print_ranges=True):\n    # Get data from dcm\n    window_center, window_width, intercept, slope = get_windowing(img)\n    img = img.pixel_array\n    \n    # Window based on dcm metadata\n    img = img * slope + intercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    if print_ranges:\n        print(img_min, img_max)\n    img[img < img_min] = img_min\n    img[img > img_max] = img_max\n    \n    # Normalize\n    img = (img - img_min) / (img_max - img_min)\n    return img\n    \n\nprint('Metadata Window Ranges:')\nview_images([metadata_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looks like the metadata ranges are somewhat similar to the brain-range we initially used.\n\n### One Window, Three Channels <a></a>\nSince we'd like to eventually export the scans as png files, we have 3 channels (R,G,B) to work with. If we're only going to use one window setting, we can try to improve the contrast by spreading it out across all 3 channels."},{"metadata":{"trusted":true},"cell_type":"code","source":"def all_channels_window(img):\n    grey_img = brain_window(img) * 3.0\n    all_chan_img = np.zeros((grey_img.shape[0], grey_img.shape[1], 3))\n    all_chan_img[:, :, 2] = np.clip(grey_img, 0.0, 1.0)\n    all_chan_img[:, :, 0] = np.clip(grey_img - 1.0, 0.0, 1.0)\n    all_chan_img[:, :, 1] = np.clip(grey_img - 2.0, 0.0, 1.0)\n    return all_chan_img\n    \n\nview_images([all_channels_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Gradient Windowing <a></a>\nWe can spread our a single window across channels a different way, by mapping the pixel values to a gradient."},{"metadata":{"trusted":true},"cell_type":"code","source":"def map_to_gradient(grey_img):\n    rainbow_img = np.zeros((grey_img.shape[0], grey_img.shape[1], 3))\n    rainbow_img[:, :, 0] = np.clip(4 * grey_img - 2, 0, 1.0) * (grey_img > 0) * (grey_img <= 1.0)\n    rainbow_img[:, :, 1] =  np.clip(4 * grey_img * (grey_img <=0.75), 0,1) + np.clip((-4*grey_img + 4) * (grey_img > 0.75), 0, 1)\n    rainbow_img[:, :, 2] = np.clip(-4 * grey_img + 2, 0, 1.0) * (grey_img > 0) * (grey_img <= 1.0)\n    return rainbow_img\n\ndef rainbow_window(img):\n    grey_img = brain_window(img)\n    return map_to_gradient(grey_img)\n\nview_images([rainbow_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Brain + Subdural + Bone Windowing <a></a>\nAs David points out in his Kernel, different hemmorhages become more obvious at different window settings.. We can include more than one window in our training images by storing a different window in each channel."},{"metadata":{"trusted":true},"cell_type":"code","source":"def window_image(img, window_center, window_width):\n    _, _, intercept, slope = get_windowing(img)\n    img = img.pixel_array * slope + intercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img[img < img_min] = img_min\n    img[img > img_max] = img_max\n    img = (img - np.min(img)) / (np.max(img) - np.min(img))\n    return img\n\n\ndef bsb_window(img):\n    brain_img = window_image(img, 40, 80)\n    subdural_img = window_image(img, 80, 200)\n    bone_img = window_image(img, 600, 2000)\n    \n    bsb_img = np.zeros((brain_img.shape[0], brain_img.shape[1], 3))\n    bsb_img[:, :, 0] = brain_img\n    bsb_img[:, :, 1] = subdural_img\n    bsb_img[:, :, 2] = bone_img\n    return bsb_img\n\nview_images([bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Exclusive Windowing <a></a>\nSame idea as above, but removing values outside the window range by setting them to 0."},{"metadata":{"trusted":true},"cell_type":"code","source":"def window_image_bottom(img, window_center, window_width):\n    _, _, intercept, slope = get_windowing(img)\n    img = img.pixel_array * slope + intercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img[img < img_min] = img_min\n    img[img > img_max] = img_min\n    img = (img - np.min(img)) / (np.max(img) - np.min(img))\n    return img\n\n\ndef bsb_window(img):\n    brain_img = window_image_bottom(img, 40, 80)\n    subdural_img = window_image_bottom(img, 80, 200)\n    bone_img = window_image_bottom(img, 600, 2000)\n    \n    bsb_img = np.zeros((brain_img.shape[0], brain_img.shape[1], 3))\n    bsb_img[:, :, 0] = brain_img\n    bsb_img[:, :, 1] = subdural_img\n    bsb_img[:, :, 2] = bone_img\n    return bsb_img\n\nview_images([bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Gradient (Brain + Subdural + Bone) Windowing <a></a>\nWe can combine a few previous ideas by averaging 3 different window settings and then mapping the results to a gradient. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def rainbow_bsb_window(img):\n    brain_img = window_image(img, 40, 80)\n    subdural_img = window_image(img, 80, 200)\n    bone_img = window_image(img, 600, 2000)\n    combo = (brain_img*0.3 + subdural_img*0.5 + bone_img*0.2)\n    return map_to_gradient(combo)\n\nview_images([rainbow_bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Sigmoid Windowing <a></a>\nInstead of simply clipping values outside of our window, we can use a sigmoid function to increase the variance near the middle of the window, while limiting the variance at the extremes. (I didn't come up with this, I believe it's a pretty common way to window.)\n\nThis looks like it does a better job creating contrast near the center of the window and we aren't losing any data like we do when we clip to a min/max."},{"metadata":{"trusted":true},"cell_type":"code","source":"def sigmoid_window(img, window_center, window_width, U=1.0, eps=(1.0 / 255.0)):\n    _, _, intercept, slope = get_windowing(img)\n    img = img.pixel_array * slope + intercept\n    ue = np.log((U / eps) - 1.0)\n    W = (2 / window_width) * ue\n    b = ((-2 * window_center) / window_width) * ue\n    z = W * img + b\n    img = U / (1 + np.power(np.e, -1.0 * z))\n    img = (img - np.min(img)) / (np.max(img) - np.min(img))\n    return img\n\ndef sigmoid_brain_window(img):\n    return sigmoid_window(img, 40, 80)\n\nview_images([sigmoid_brain_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Sigmoid (Brain + Subdural + Bone) Windowing <a></a>\nAgain, combining two previous ideas."},{"metadata":{"trusted":true},"cell_type":"code","source":"def sigmoid_bsb_window(img):\n    brain_img = sigmoid_window(img, 40, 80)\n    subdural_img = sigmoid_window(img, 80, 200)\n    bone_img = sigmoid_window(img, 600, 2000)\n    \n    bsb_img = np.zeros((brain_img.shape[0], brain_img.shape[1], 3))\n    bsb_img[:, :, 0] = brain_img\n    bsb_img[:, :, 1] = subdural_img\n    bsb_img[:, :, 2] = bone_img\n    return bsb_img\n\nview_images([sigmoid_bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Sigmoid Gradient (Brain + Subdural + Bone) Windowing <a></a>\nAnd finally, putting it all together."},{"metadata":{"trusted":true},"cell_type":"code","source":"def map_to_gradient_sig(grey_img):\n    rainbow_img = np.zeros((grey_img.shape[0], grey_img.shape[1], 3))\n    rainbow_img[:, :, 0] = np.clip(4*grey_img - 2, 0, 1.0) * (grey_img > 0.01) * (grey_img <= 1.0)\n    rainbow_img[:, :, 1] =  np.clip(4*grey_img * (grey_img <=0.75), 0,1) + np.clip((-4*grey_img + 4) * (grey_img > 0.75), 0, 1)\n    rainbow_img[:, :, 2] = np.clip(-4*grey_img + 2, 0, 1.0) * (grey_img > 0.01) * (grey_img <= 1.0)\n    return rainbow_img\n\ndef sigmoid_rainbow_bsb_window(img):\n    brain_img = sigmoid_window(img, 40, 80)\n    subdural_img = sigmoid_window(img, 80, 200)\n    bone_img = sigmoid_window(img, 600, 2000)\n    combo = (brain_img*0.35 + subdural_img*0.5 + bone_img*0.15)\n    combo_norm = (combo - np.min(combo)) / (np.max(combo) - np.min(combo))\n    return map_to_gradient_sig(combo_norm)\n\nview_images([sigmoid_rainbow_bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I haven't done enough experimenting to know which one of these works best. (Although I have a hunch it's *Sigmoid (Brain + Subdural + Bone) Windowing*.) \n\nI was thinking these might help add some diversity to an ensemble?\n\nThanks for reading and let me know if you have any questions or comments! **Also, if you found any of this interesting, don't forget to upvote.** 😄\n\n## Acknowledgements <a></a>\nThanks again to [David Tang](https://www.kaggle.com/dcstang/see-like-a-radiologist-with-systematic-windowing), [Marco](https://www.kaggle.com/marcovasquez/basic-eda-data-visualization), [Nanashi](https://www.kaggle.com/jesucristo/rsna-introduction-eda-models), and [Richard McKinley](https://www.kaggle.com/omission/eda-view-dicom-images-with-correct-windowing) for sharing their excellent kernels.\n\nIdeas were also borrowed from [Practical Window Setting Optimization for Medical Image Deep Learning](https://arxiv.org/pdf/1812.00572.pdf) and [Precise diagnosis of intracranial hemorrhage and subtypes using a three-dimensional joint convolutional and recurrent neural network](https://rd.springer.com/content/pdf/10.1007%2Fs00330-019-06163-2.pdf)"},{"metadata":{},"cell_type":"markdown","source":"**Update 10/22: ** The `sigmoid_window` function can be pretty slow when processing the whole dataset on a cpu. You can get a nice speedup by running it on your gpu using cupy:"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"import cupy as cp\n\ndef sigmoid_window(dcm, img, window_center, window_width, U=1.0, eps=(1.0 / 255.0)):\n    img = cp.array(np.array(img))\n    _, _, intercept, slope = get_windowing(dcm)\n    img = img * slope + intercept\n    ue = cp.log((U / eps) - 1.0)\n    W = (2 / window_width) * ue\n    b = ((-2 * window_center) / window_width) * ue\n    z = W * img + b\n    img = U / (1 + cp.power(np.e, -1.0 * z))\n    img = (img - cp.min(img)) / (cp.max(img) - cp.min(img))\n    return cp.asnumpy(img)","execution_count":null,"outputs":[]}],"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}