{"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":"## How to manually apply windowing to an x-ray image.\n\n- In this notebook, we'll grab an image from the SIIM Covid19 dataset and manually apply a Window Width and Window Level to it.\n\nNormally, the pixel ranges of x-rays are much larger than the pixel display range of consumer grade monitors. 10-16 bit images on 8 bit displays.\n\nSince we can only display 256 shades of gray (8 bit) on a most monitors, we must map *many-to-one* pixels here using a Look Up Table (LUT).\n\nAutomatically calculating the width/center does not always produce the best diagnostic quality image for the pathology in question.\n\nDICOM image viewers typically have a tool for windowing that allows radiologists to adjust the levels from the full image, rather than the 8 bit representation of it. This allows them to see subtle differences that they would not otherwise be able to see.\n\nSetting the Window Width and Window Level creates a VOI (Values of Interest) LUT and applies it to the image. This is a lossy operation.\n\nSome DICOM images contain VOI LUT arrays, reference external VOI LUTs, or use specified Width/Level tags to calculate the LUT.\n\n** VOI LUT should not be confused with Modality LUT, which maps modality unit specific values to usable pixel values. If a modality LUT is present in the DICOM image, it should always be applied first in the pipeline.*\n\n- Here's a couple DICOM pre-processing notebooks I made:\n- Rib supression on Chest X-Rays -> https://www.kaggle.com/davidbroberts/rib-suppression-poc\n- Apply Unsharp Mask to Chest X-Rays -> https://www.kaggle.com/davidbroberts/unsharp-masking-chest-x-rays\n- Cropping Chest X-Rays -> https://www.kaggle.com/davidbroberts/cropping-chest-x-rays\n- Visualizing Chest X-Ray bit planes -> https://www.kaggle.com/davidbroberts/visualizing-chest-x-ray-bitplanes\n- DICOM full range pixels as CNN input -> https://www.kaggle.com/davidbroberts/dicom-full-range-pixels-as-cnn-input","metadata":{}},{"cell_type":"code","source":"import os\nimport cv2\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport matplotlib.pyplot as plt\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:35:30.591647Z","iopub.execute_input":"2021-05-30T13:35:30.592298Z","iopub.status.idle":"2021-05-30T13:35:30.602138Z","shell.execute_reply.started":"2021-05-30T13:35:30.592229Z","shell.execute_reply":"2021-05-30T13:35:30.599796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Create a couple helper functions","metadata":{}},{"cell_type":"code","source":"# This function gets the first image path in a StudyInstanceUID directory\ndef get_image_by_study_id(study_id):\n    base_path = \"/kaggle/input/siim-covid19-detection/\"\n    study_path = base_path + \"train/\" + study_id + \"/\"\n    images = []\n    for subdir, dirs, files in os.walk(study_path):\n        for file in files:     \n            image = os.path.join(subdir, file)\n            if os.path.isfile(image):\n                return image\n    return \"none\"","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:35:30.604461Z","iopub.execute_input":"2021-05-30T13:35:30.604881Z","iopub.status.idle":"2021-05-30T13:35:30.630551Z","shell.execute_reply.started":"2021-05-30T13:35:30.604821Z","shell.execute_reply":"2021-05-30T13:35:30.629477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Make a simple linear VOI LUT from the raw (stored) pixel data\ndef make_lut(storedPixels, windowWidth, windowLevel, p_i):\n    \n    # Slope and Intercept set to 1 and 0 for X-ray. Get these from DICOM tags instead if using \n    # on a modality that requires them (CT, PT etc)\n    slope = 1.0\n    intercept = 0.0\n    minPixel = int(np.amin(storedPixels))\n    maxPixel = int(np.amax(storedPixels))\n\n    # Make an empty array for the LUT the size of the pixel 'width' in the raw pixel data\n    lut = [0] * (maxPixel + 1)\n    \n    # Invert pixels and windowLevel for MONOCHROME1. We invert the specified windowLevel so that \n    # increasing the level value makes the images brighter regardless of photometric intrepretation\n    invert = False\n    if p_i == \"MONOCHROME1\":\n        invert = True\n    else:\n        windowLevel = (maxPixel - minPixel) - windowLevel\n        \n    # Loop through the pixels and calculate each LUT value\n    for storedValue in range(minPixel, maxPixel):\n        modalityLutValue = storedValue * slope + intercept\n        voiLutValue = (((modalityLutValue - windowLevel) / windowWidth + 0.5) * 255.0)\n        clampedValue = min(max(voiLutValue, 0), 255)\n        if invert:\n            lut[storedValue] = round(255-clampedValue)\n        else:\n            lut[storedValue] = round(clampedValue)\n        \n    return lut","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:35:30.633313Z","iopub.execute_input":"2021-05-30T13:35:30.634086Z","iopub.status.idle":"2021-05-30T13:35:30.652888Z","shell.execute_reply.started":"2021-05-30T13:35:30.634033Z","shell.execute_reply":"2021-05-30T13:35:30.651242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply the LUT to a pixel array\ndef apply_lut(pixels_in, lut):\n    pixels_in = pixels_in.flatten()\n    pixels_out = [0] * len(pixels_in)\n    for i in range(0, len(pixels_in)):\n        pixel = pixels_in[i]\n        pixels_out[i] = int(lut[pixel])\n    return pixels_out","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:35:30.655937Z","iopub.execute_input":"2021-05-30T13:35:30.657344Z","iopub.status.idle":"2021-05-30T13:35:30.676481Z","shell.execute_reply.started":"2021-05-30T13:35:30.657271Z","shell.execute_reply":"2021-05-30T13:35:30.674643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Display an image with default windowing.","metadata":{}},{"cell_type":"code","source":"# Take a look at a random image and display it with imshow. The width and level are automagically calculated from the pixel data.\nimg_file = get_image_by_study_id(\"00c241c3fc0d\")\nprint(\"Loading image: \" + img_file)\n\nimg = pydicom.dcmread(img_file)\nplt.imshow(img.pixel_array,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:35:30.679195Z","iopub.execute_input":"2021-05-30T13:35:30.680268Z","iopub.status.idle":"2021-05-30T13:35:31.322152Z","shell.execute_reply.started":"2021-05-30T13:35:30.680222Z","shell.execute_reply":"2021-05-30T13:35:31.320769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Calculate the width and level from the pixel ranges manually.","metadata":{}},{"cell_type":"code","source":"# Set the width to be the distance between the pixels, and the level to be the middle of the pixel range.\n# The resulting image should look exactly like the one above.\npixels = img.pixel_array\nminPixel = np.min(pixels)\nmaxPixel = np.max(pixels)\nwindowWidth = maxPixel - minPixel\nwindowLevel = (minPixel + maxPixel) / 2\n\nlut = make_lut(pixels, windowWidth, windowLevel, img.PhotometricInterpretation)\npixels = apply_lut(pixels, lut)\n\n# Reshape the pixel array back into the image shape\nimg_out = np.reshape(pixels, (img.pixel_array.shape[0],img.pixel_array.shape[1]))\n\nprint(\"Calculated - Width: \" + str(windowWidth) + \" / Level  \" +  str(windowLevel))\nprint(\"Pixel range: \" + str(minPixel) + \" - \" + str(maxPixel))\nplt.imshow(img_out,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:42:55.405679Z","iopub.execute_input":"2021-05-30T13:42:55.406079Z","iopub.status.idle":"2021-05-30T13:42:59.305232Z","shell.execute_reply.started":"2021-05-30T13:42:55.406048Z","shell.execute_reply":"2021-05-30T13:42:59.302624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- In this case, the image pixel range is 0-4095, which indicates this is probably a 12 bit image.\n- Since there isn't a 12 bit data type, we have to store these in 16 bit arrays.\n- The Photometric Interpretation is MONOCHROME2, meaning the pixel intensities get brighter as the pixel value increases.\n- MNONOCHROME1 images, the intensities get brighter as the pixel value decreases.\n- We can verify by looking at a couple DICOM tags.","metadata":{}},{"cell_type":"code","source":"print(\"PhotometricInterpretation: \" + img.PhotometricInterpretation)\nprint(\"Bit stored: \" + str(img.BitsStored))","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:35:35.245402Z","iopub.execute_input":"2021-05-30T13:35:35.246149Z","iopub.status.idle":"2021-05-30T13:35:35.252859Z","shell.execute_reply.started":"2021-05-30T13:35:35.246109Z","shell.execute_reply":"2021-05-30T13:35:35.251504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Manually specify window width and level.","metadata":{}},{"cell_type":"code","source":"# If we use windowWidth = 1000, and keep the Level at center .. the image will no longer be 'full width' \n# Setting the width to less than the distance between the output pixels means the full range of pixels \n# is not used. The resulting image displays less than 255 shades of grey. If the width was set to '2',\n# the image would only display black or white pixels.\n\npixels = img.pixel_array\nwindowWidth = 1000\nwindowLevel = 2047\n\nlut = make_lut(pixels, windowWidth, windowLevel, img.PhotometricInterpretation)\npixels = apply_lut(pixels, lut)\n\nprint(\"Width: \" + str(windowWidth) + \" / Level  \" +  str(windowLevel))\nimg_out = np.reshape(pixels, (img.pixel_array.shape[0],img.pixel_array.shape[1]))\nplt.imshow(img_out,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:36:50.569894Z","iopub.execute_input":"2021-05-30T13:36:50.570384Z","iopub.status.idle":"2021-05-30T13:36:54.369589Z","shell.execute_reply.started":"2021-05-30T13:36:50.570346Z","shell.execute_reply":"2021-05-30T13:36:54.367948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's make this image full width, but a little brighter\n\npixels = img.pixel_array\nwindowWidth = 4095\nwindowLevel = 3500\n\nlut = make_lut(pixels, windowWidth, windowLevel, img.PhotometricInterpretation)\npixels = apply_lut(pixels, lut)\n\nprint(\"Manual Width: \" + str(windowWidth) + \" / Level  \" +  str(windowLevel))\nimg_out = np.reshape(pixels, (img.pixel_array.shape[0],img.pixel_array.shape[1]))\nplt.imshow(img_out,cmap=\"gray\");","metadata":{"execution":{"iopub.status.busy":"2021-05-30T13:38:59.201251Z","iopub.execute_input":"2021-05-30T13:38:59.201631Z","iopub.status.idle":"2021-05-30T13:39:03.059274Z","shell.execute_reply.started":"2021-05-30T13:38:59.2016Z","shell.execute_reply":"2021-05-30T13:39:03.058169Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Conclusion\n- Creating image sets with various Widths and Levels might offer more 'visible' information when exporting to 8 bit for model import.\n- You can try various width/level settings to best demonstrate the pathology in question.","metadata":{}}]}