{"cells":[{"metadata":{"_uuid":"81a2d52b-646a-4079-bfaa-56f87a85d8da","_cell_guid":"f007641c-6e13-42ce-81b2-a7d2ae737a4d","trusted":true},"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), [Richard McKinley](https://www.kaggle.com/omission/eda-view-dicom-images-with-correct-windowing) and [Ryan Epp](https://www.kaggle.com/reppic/gradient-sigmoid-windowing) for their amazing Kernels that I totally borrowed code and ideas from. Though the original notebook and idea of gradient and sigmoid windowing was applied for [RSNA Intracranial Hemorrhage Detection Challenge](https://www.kaggle.com/c/rsna-intracranial-hemorrhage-detection), in my view these techniques can put serious impact on the competition's LB. \n\n**Contents**<br>\n&#10687; [1. No Windowing](#1)<br>\n&#10687; [2. Lung Windowing](#2)<br>\n&#10687; [3. Metadata Windowing](#3)<br>\n&#10687; [4. One Window, Three Channels](#4)<br>\n&#10687; [5. Gradient Windowing](#5)<br>\n&#10687; [6. Mapping Multiple Windows in Multiple Channels](#6)<br>\n&#10687; [7. Exclusive Windowing](#7)<br>\n&#10687; [8. Gradient Ensemble Windowing](#8)<br>\n&#10687; [9. Sigmoid Windowing](#9)<br>\n&#10687; [10. Sigmoid Ensemble Windowing](#10)<br>\n&#10687; [11. Sigmoid Gradient Ensemble Windowing](#11)<br>\n&#10687; [12. Acknowledgements](#12)<br>"},{"metadata":{"_uuid":"edfdeb60-1867-47f1-a10b-f5d233abcacf","_cell_guid":"f9fe094e-06ac-495e-b67e-7aa94fb915cb","trusted":true},"cell_type":"markdown","source":"### Importing Packages and Setting Data Path"},{"metadata":{"_uuid":"e77cfb48-a3a6-4e9b-9fc3-55246cb7da77","_cell_guid":"96bddf81-e312-4553-9923-cd4e59a2706d","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 as dcm\nimport os\nimport glob\nimport random\n\nTRAIN_IMG_PATH = \"../input/rsna-str-pulmonary-embolism-detection/train/\"\nTEST_IMG_PATH = \"../input/rsna-str-pulmonary-embolism-detection/test/\"\n\ntrain_df = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")\ntest_df = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/test.csv\")\n\n# Uncomment this lines to get a list of paths of all training and testing files\n# train_file_paths = glob.glob(\"../input/rsna-str-pulmonary-embolism-detection/train/*/*/*.dcm\")\n# test_file_paths = glob.glob(\"../input/rsna-str-pulmonary-embolism-detection/test/*/*/*.dcm\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1cd8e093-7238-48cb-9ea3-5760b889427d","_cell_guid":"969d1e66-7209-4262-b472-62b6feac93e8","trusted":true},"cell_type":"markdown","source":"Let's select an experiment or examination which has both positive and negative sample. Here an experiment is handpicked so that it has half positive samples and half negative samples."},{"metadata":{"_uuid":"cecbf4a4-bbda-4a84-9c02-7a078832087a","_cell_guid":"33321e15-d889-4e64-a4cb-c638e41e1cd5","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"exam_UID = '6897fa9de148'\nexam = train_df[train_df[\"StudyInstanceUID\"] == exam_UID]\npositive_exam_paths = []\nnegative_exam_paths = []\n\nfor index, row in exam.iterrows():\n    if row[\"pe_present_on_image\"]:\n        positive_exam_paths.append('../input/rsna-str-pulmonary-embolism-detection/train/' + \n                                  row[\"StudyInstanceUID\"] + '/' +\n                                  row['SeriesInstanceUID'] + '/' + \n                                  row['SOPInstanceUID'] + '.dcm')\n    else:\n        negative_exam_paths.append('../input/rsna-str-pulmonary-embolism-detection/train/' + \n                                  row[\"StudyInstanceUID\"] + '/' +\n                                  row['SeriesInstanceUID'] + '/' + \n                                  row['SOPInstanceUID'] + '.dcm')\n\nprint(f'Statistics for examination UID: {exam_UID}')\nprint('pe_present_on_image -> (Positive): ', len(positive_exam_paths))\nprint('pe_present_on_image -> (Negative): ', len(negative_exam_paths))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d19dbc0c-70c9-4583-b740-ac50a67ec04e","_cell_guid":"c3836333-61f3-4ade-89b4-2b42064cc877","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def load_random_images(no_of_samples = 10):    \n    img_paths = random.sample(positive_exam_paths, no_of_samples//2)\n    img_paths += random.sample(negative_exam_paths, no_of_samples//2)\n    return [dcm.read_file(img_path) for img_path in img_paths]\n\n\ndef view_images(images, cmap = plt.cm.bone ):\n    width = len(images)//2\n    height = 2\n    fig, axs = plt.subplots(height, width, figsize=(20,8))\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=cmap) \n        axs[i,j].axis('off')\n        if i==0:\n            title = 'pe_present=1'\n        else:\n            title = 'pe_present=0'\n        axs[i,j].set_title(title)\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) == dcm.multival.MultiValue: return int(x[0])\n    else: 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 utility functions successfully')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d9b74ac1-7713-4608-bd65-94b53c8c419c","_cell_guid":"70206d36-e26a-4a49-a5dc-6ce2e9013987","trusted":true},"cell_type":"markdown","source":"Let's load a some images. We'll randomly pick five from positive example from each class and five negative examples. We will see the same examples through out the whole notebook so it is easier to distinguish them."},{"metadata":{"_uuid":"f9ca24e3-67fb-477b-896d-7c091e0a0d0f","_cell_guid":"b9e5fd2a-7949-4b90-b01e-eff11e32edaa","trusted":true},"cell_type":"code","source":"imgs = load_random_images()\nprint('DICOM Pixel Data Shape: ', [img_pixels.pixel_array.shape[:]  for img_pixels in imgs])\n\nplt.title('Distribution of DICOM Pixel Values')\nax = plt.hist(np.array([img.pixel_array for img in imgs]).flatten(), bins=300, color='c')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25c03676-194f-4caf-9de0-5a85edce91cd","_cell_guid":"9667a568-a4f6-4ec6-b12e-bb579bdf4bbd","trusted":true},"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<a id=\"1\"></a>\n<font color=\"black\" size=+2.5><b>1. No Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\n### Default color model <font color=\"blue\">plt.cm.bone</font>"},{"metadata":{"_uuid":"33d6e12e-5a8d-4318-9981-83446cae9917","_cell_guid":"0c2fb3a9-6fdf-41a7-943e-b7102cc3fe24","trusted":true},"cell_type":"code","source":"view_images([img.pixel_array for img in imgs] )","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25aaab31-f940-48c5-b3cd-250159682dfc","_cell_guid":"51be10f8-a953-4cb3-8a4e-dc86a2abd962","trusted":true},"cell_type":"markdown","source":"### Default color model <font color=\"blue\">plt.cm.gray</font>"},{"metadata":{"_uuid":"e2e79296-cef6-4709-a8d2-44c70abeb538","_cell_guid":"a25313ae-9fa8-4fa2-a459-6754014fe205","trusted":true},"cell_type":"code","source":"view_images([img.pixel_array for img in imgs], cmap = 'gray')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f445a088-766c-4de2-8ffc-f92673f7f12f","_cell_guid":"c1df24c0-b40e-429d-aebf-76532eb9a9a4","trusted":true},"cell_type":"markdown","source":"These don't look very useful for pulmonary embolism.\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 PE?**"},{"metadata":{"_uuid":"6bdce552-448d-4e8d-929b-50eb424d67d7","_cell_guid":"54cb0a0a-8d4b-4f25-a8cb-21544c0b96db","trusted":true},"cell_type":"markdown","source":"<a id=\"2\"></a>\n<font color=\"black\" size=+2.5><b>2. Lung Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\nLet's check out the range of values that corresponds with lung tissues. I take values from [radiopedia.org](https://radiopaedia.org/articles/windowing-ct) We'll clip everything outside that range so that there's more contrast in the lung-tissue range."},{"metadata":{"_uuid":"4d811eca-b0cc-4807-b608-8b4cf929469c","_cell_guid":"1bfefc83-6238-47da-bf06-ca7fb5d3a4c3","trusted":true},"cell_type":"code","source":"def lung_window(img):\n    window_min = -125\n    window_max = 225\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([lung_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"26daee85-4bfe-4285-b5e7-fd3424e1244b","_cell_guid":"4f02e8ba-e5bd-4f4a-ac01-12cde1d2495d","trusted":true},"cell_type":"markdown","source":"<a id=\"3\"></a>\n<font color=\"black\" size=+2.5><b>3. Metadata Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\nThe DICOM images come with metadata specifying a window center and width. Which was used when that CT-Scan was taken. We could also use these values instead of the fixed range from above. We can simply extract the metadata from the DICOM files. However for the first time, we will use that default windowing and later on we will use different windowing parameters based on our need."},{"metadata":{"_uuid":"a1a04e28-fb55-4a09-83e9-9082f506a0c2","_cell_guid":"d25ea471-dd53-499b-b8f8-b49763be96e1","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(f'Window Level: {window_center}, Window Width:{window_width}, Range:[{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":{"_uuid":"0110ab97-1dcb-4daa-a7d0-1855e0699919","_cell_guid":"19ce8af7-cc3c-44f1-a914-359d62606ce2","trusted":true},"cell_type":"markdown","source":"Looks like the metadata ranges are somewhat similar to the lung-range we initially used. So we will use that in our calculation throughout. \n\n<a id=\"4\"></a>\n<font color=\"black\" size=+2.5><b>4. One Window, Three Channels </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\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":{"_uuid":"815393e3-6c9b-4426-989f-066bed592ee4","_cell_guid":"4d12e109-8bbb-45f0-8515-4b1298924d0e","trusted":true},"cell_type":"code","source":"def all_channels_window(img):\n    grey_img = lung_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":{"_uuid":"d7754bde-8fc7-4727-860e-6dd44e47373f","_cell_guid":"439e80d0-4a67-46cc-bc19-4b2217121d38","trusted":true},"cell_type":"markdown","source":"<a id=\"5\"></a>\n<font color=\"black\" size=+2.5><b>5. Gradient Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\nWe can spread our a single window across channels a different way, by mapping the pixel values to a gradient."},{"metadata":{"_uuid":"e5052588-2b03-4e8a-af00-ab8ccca96781","_cell_guid":"195f91ad-9c24-4b55-9ce3-b77babbe9752","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 = lung_window(img)\n    return map_to_gradient(grey_img)\n\nview_images([rainbow_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"56d9e446-6ffc-4eb2-81aa-418d5ade2b7e","_cell_guid":"4ca5d72a-da8b-4fcb-9f28-86156d3b3078","trusted":true},"cell_type":"markdown","source":"<a id=\"6\"></a>\n<font color=\"black\" size=+2.5><b>6. Mapping Multiple Windows in Multiple Channels </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\n* As David points out in his Kernel, different density tissue 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":{"_uuid":"243da8f0-9778-4248-a046-1eed5bc27b59","_cell_guid":"256d6a33-878b-49aa-a8af-61bf4fe35cf9","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    lung_img1 = window_image(img, 40, 80)\n    lung_img2 = window_image(img, 80, 200)\n    lung_img3 = window_image(img, 600, 2000)\n    \n    bsb_img = np.zeros((lung_img1.shape[0], lung_img1.shape[1], 3))\n    bsb_img[:, :, 0] = lung_img1\n    bsb_img[:, :, 1] = lung_img2\n    bsb_img[:, :, 2] = lung_img3\n    return bsb_img\n\nview_images([bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ed426bd1-396c-440b-8422-2e290adac28b","_cell_guid":"e6056710-f764-4d12-bc76-1c8c1091aef8","trusted":true},"cell_type":"markdown","source":"<a id=\"7\"></a>\n<font color=\"black\" size=+2.5><b>7. Exclusive Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\nSame idea as above, but removing values outside the window range by setting them to 0."},{"metadata":{"_uuid":"e40a0346-d1ce-476e-998c-0492070d9d29","_cell_guid":"ef91e97d-2e93-4600-8b06-f9e6058de13e","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    lung_img1 = window_image_bottom(img, 40, 80)\n    lung_img2 = window_image_bottom(img, 80, 200)\n    lung_img3 = window_image_bottom(img, 600, 2000)\n    \n    bsb_img = np.zeros((lung_img1.shape[0], lung_img1.shape[1], 3))\n    bsb_img[:, :, 0] = lung_img1\n    bsb_img[:, :, 1] = lung_img2\n    bsb_img[:, :, 2] = lung_img3\n    return bsb_img\n\nview_images([bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c69e70ba-6896-430b-9000-13e7fffbba27","_cell_guid":"10062567-6509-4279-a23e-8047d1e62dcd","trusted":true},"cell_type":"markdown","source":"<a id=\"8\"></a>\n<font color=\"black\" size=+2.5><b>8. Gradient Ensemble Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\nWe can combine a few previous ideas by averaging 3 different window settings and then mapping the results to a gradient."},{"metadata":{"_uuid":"18954af3-596d-4b29-9e76-047331291574","_cell_guid":"378db12f-4eca-41ea-8c8a-92f867d3ece1","trusted":true},"cell_type":"code","source":"def rainbow_bsb_window(img):\n    lung_img1 = window_image(img, 40, 80)\n    lung_img2 = window_image(img, 80, 200)\n    lung_img3 = window_image(img, 600, 2000)\n    combo = (lung_img1*0.3 + lung_img2*0.5 + lung_img3*0.2)\n    return map_to_gradient(combo)\n\nview_images([rainbow_bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"073e5f83-c548-4c39-9d93-cf45e26753d4","_cell_guid":"3f90a9c7-ee11-488c-b3c4-bb52f8dc573b","trusted":true},"cell_type":"markdown","source":"<a id=\"9\"></a>\n<font color=\"black\" size=+2.5><b>9. Sigmoid Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\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":{"_uuid":"bd9ff063-8592-49f8-a8bd-80b245879053","_cell_guid":"c0e54085-0780-4ae0-a10d-e2bef72c660f","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":{"_uuid":"aee12558-cb3f-42cf-a300-5b8306481fd3","_cell_guid":"ec3b1a51-b99e-4fce-9c4e-15dd180b47e3","trusted":true},"cell_type":"markdown","source":"<a id=\"10\"></a>\n<font color=\"black\" size=+2.5><b>10. Sigmoid Ensemble Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\nAgain, combining two previous ideas."},{"metadata":{"_uuid":"bcb59e75-dcc6-4ec1-b853-78592c1fce96","_cell_guid":"910f61af-d10c-4f88-8b74-f28cdb8a2bcd","trusted":true},"cell_type":"code","source":"def sigmoid_bsb_window(img):\n    lung_img1 = sigmoid_window(img, 40, 80)\n    lung_img2 = sigmoid_window(img, 80, 200)\n    lung_img3 = sigmoid_window(img, 600, 2000)\n    \n    bsb_img = np.zeros((lung_img1.shape[0], lung_img1.shape[1], 3))\n    bsb_img[:, :, 0] = lung_img1\n    bsb_img[:, :, 1] = lung_img2\n    bsb_img[:, :, 2] = lung_img3\n    return bsb_img\n\nview_images([sigmoid_bsb_window(img) for img in imgs])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e315b77c-d511-437f-94b0-a6315adf5fa7","_cell_guid":"cb6dd011-729c-4d85-8661-3f6c6b0dcd63","trusted":true},"cell_type":"markdown","source":"<a id=\"11\"></a>\n<font color=\"black\" size=+2.5><b>11. Sigmoid Gradient Ensemble Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\n\nAnd finally, putting it all together."},{"metadata":{"_uuid":"3b7a7788-c808-413a-8038-05a625738e01","_cell_guid":"6338cbd7-1545-4b80-93ae-5c7678d85fd5","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    lung_img1 = sigmoid_window(img, 40, 80)\n    lung_img2 = sigmoid_window(img, 80, 200)\n    lung_img3 = sigmoid_window(img, 600, 2000)\n    combo = (lung_img1*0.35 + lung_img2*0.5 + lung_img3*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":{"_uuid":"b1793808-4836-4f5e-8991-64ba0c341a77","_cell_guid":"1715abac-d18f-45ee-9f4f-2e22df2732cd","trusted":true},"cell_type":"markdown","source":"<a id=\"12\"></a>\n<font color=\"black\" size=+2.5><b>12. Channelwise Windowing </b></font>\n\n<a href=\"#top\" class=\"btn btn-primary btn-sm\" role=\"button\" aria-pressed=\"true\" style=\"color:white\" data-toggle=\"popover\" title=\"go to Colors\">Go to TOC</a>\n\nThe values in CT scans tend to range from -1000 to 3000; we are used to 8-bit images with pixel values ranging from 0 to 255. Radiologists use a technique called windowing to visualize CT scans. Different types of tissues are better evaluated using different windows. Windows are defined by 2 numbers: window width and window level.\n\n\nNote that these images are single channel. I have provided them in 3-channel RGB format. Each channel is a different window.\n\n* RED channel / LUNG window / level=-600, width=1500\n* GREEN channel / PE window / level=100, width=700\n* BLUE channel / MEDIASTINAL window / level=40, width=400\n\n"},{"metadata":{"_uuid":"414c9c7d-737e-4a28-b494-adb6573c9170","_cell_guid":"1eb78635-28b6-446d-8cc7-73dd232e3406","trusted":true},"cell_type":"code","source":"def window(img, WL=50, WW=350):\n    upper, lower = WL+WW//2, WL-WW//2\n    X = np.clip(img.copy(), lower, upper)\n    X = X - np.min(X)\n    X = X / np.max(X)\n    X = (X*255.0).astype('uint8')\n    return X\n\nimage_raw = imgs[4].pixel_array\n\nimage_lung = np.expand_dims(window(image_raw, WL=-600, WW=1500), axis=3)\nimage_mediastinal = np.expand_dims(window(image_raw, WL=40, WW=400), axis=3)\nimage_pe_specific = np.expand_dims(window(image_raw, WL=100, WW=700), axis=3)\n# image = np.concatenate([image_mediastinal, image_pe_specific, image_lung], axis=3)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7f58262d-61c6-4c7e-aeb5-d1290a1a9b22","_cell_guid":"cd79158e-b716-4db2-84a0-dd822f62f829","trusted":true},"cell_type":"code","source":"rgb = np.dstack((image_lung,image_mediastinal,image_pe_specific))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0781f399-2c3b-4005-a13c-786fa61e301e","_cell_guid":"15393700-c37c-485d-9e26-336980c7d6b2","trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(1,4, figsize=(20, 10))\nax[0].imshow(image_lung.squeeze(), cmap='bone')\nax[0].axis('off')\nax[0].set_title('Window: Lung')\n\nax[1].imshow(image_mediastinal.squeeze(), cmap='bone')\nax[1].axis('off')\nax[1].set_title('Window: mediastinal ')\n\nax[2].imshow(image_pe_specific.squeeze(), cmap='bone')\nax[2].axis('off')\nax[2].set_title('Window: pe_specific')\n\nax[3].imshow(rgb )\nax[3].axis('off')\nax[3].set_title('Merged RGB')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c21cd74c-51cb-4426-9ec5-565349faecde","_cell_guid":"e6c3d857-0c36-44ed-b48a-b6defaeffdd2","trusted":true},"cell_type":"markdown","source":""},{"metadata":{"_uuid":"287b2de4-f373-4542-bb43-a83540d0e17c","_cell_guid":"1cdad06f-5407-445a-8d09-f4cccc205510","trusted":true},"cell_type":"markdown","source":"I haven't done enough experimenting to know which one of these works best because I haven't submitted any prediction to the competition.  I 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), [Richard McKinley](https://www.kaggle.com/omission/eda-view-dicom-images-with-correct-windowing) and [Ryan Epp](https://www.kaggle.com/reppic/gradient-sigmoid-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":{"_uuid":"358972c0-1fac-4201-aaad-85be9afa28cb","_cell_guid":"085800e0-56be-4b74-8755-da3529862858","trusted":true},"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":{"_uuid":"df22574d-23a0-4e92-926c-70965e639dc5","_cell_guid":"189ce664-eba1-498c-a2ff-ac3d735dc211","trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"# import cupy as cp\n\n# def 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":{"_uuid":"7ac7a40f-5973-4548-9a49-cb468ca556cf","_cell_guid":"5ca4137a-075e-41d6-ba62-cb12e4bb4218","trusted":true},"cell_type":"markdown","source":"![](https://www.clipartmax.com/png/middle/265-2655834_work-in-progress-icon.png)\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)<br>\n💥 [CT-Scans, DICOM files, Windowing Explained](https://www.kaggle.com/redwankarimsony/ct-scans-dicom-files-windowing-explained)<br>\n💥 [RSNA-STR-PE [Gradient & Sigmoid Windowing]](https://www.kaggle.com/redwankarimsony/rsna-str-pe-gradient-sigmoid-windowing)<br>\n💥 [RSNA-STR [✔️3D Stacking ✔️3D Plot ✔️Segmentation]](https://www.kaggle.com/redwankarimsony/rsna-str-3d-stacking-3d-plot-segmentation/edit/run/42517982)<br>\n💥 [RSNA-STR [DICOM 👉 GIF 👉 npy]](https://www.kaggle.com/redwankarimsony/rsna-str-dicom-gif-npy)<br>\n💥 [RSNA-STR Pulmonary Embolism [EDA]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-eda)<br>"},{"metadata":{"_uuid":"81079abe-7b84-4715-801f-14f0acd1f48a","_cell_guid":"ca711bd3-165f-415b-8c61-d332c69ebaaf","trusted":true},"cell_type":"code","source":"","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":4}