{"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":"<div class='alert alert-info' style='text-align: center'><h1>Manual VOI LUT (window) an MR Image</h1>\n- yet another MR processing notebook -</div>\n\n#### The idea is to adjust the contrast of MR images prior to conversion to 8 bit to better separate tissues of like-density.\n- This mimics what radiologists do with the WW/WL tool in a DICOM image viewer.\n- It helps to distinguish tumors from other brain tissue that surrounds them.\n- Since tumors are typically more dense than the stuff around them, but only slightly, we can make them more 'enhanced'.\n- The effect is more pronounced since the skull bones have been MIP'd out of the images.\n- Must be done before the pixels are crunched down to 8 bit.\n\n1. **Determine the pixel range of the image.**\n\n2. **Create a VOI LUT from raw pixels.**\n\n3. **Apply the LUT to the image.**\n\n4. **Rejoice!**\n\n### - Do imports and define some functions","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:48.689755Z","iopub.execute_input":"2021-09-13T02:16:48.690495Z","iopub.status.idle":"2021-09-13T02:16:48.697391Z","shell.execute_reply.started":"2021-09-13T02:16:48.690442Z","shell.execute_reply":"2021-09-13T02:16:48.695784Z"},"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(pixels, width, center, p_i):\n    \n    # Slope and Intercept set to 1 and 0 for MR. 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    min_pixel = int(np.amin(pixels))\n    max_pixel = int(np.amax(pixels))\n\n    # Make an empty array for the LUT the size of the pixel 'width' in the raw pixel data\n    lut = [0] * (max_pixel + 1)\n    \n    # Invert pixels and cent for MONOCHROME1. We invert the specified center so that \n    # increasing the center value makes the images brighter regardless of photometric intrepretation\n    invert = False\n    if p_i == \"MONOCHROME1\":\n        invert = True\n    else:\n        center = (max_pixel - min_pixel) - center\n        \n    # Loop through the pixels and calculate each LUT value\n    for pix_value in range(min_pixel, max_pixel):\n        lut_value = pix_value * slope + intercept\n        voi_value = (((lut_value - center) /  width + 0.5) * 255.0)\n        clamped_value = min(max(voi_value, 0), 255)\n        if invert:\n            lut[pix_value] = round(255 - clamped_value)\n        else:\n            lut[pix_value] = round(clamped_value)\n        \n    return lut","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:48.83962Z","iopub.execute_input":"2021-09-13T02:16:48.840082Z","iopub.status.idle":"2021-09-13T02:16:48.84873Z","shell.execute_reply.started":"2021-09-13T02:16:48.840044Z","shell.execute_reply":"2021-09-13T02:16:48.847877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply the LUT to a pixel array\ndef apply_lut(pixels_in, lut):\n    \n    pixels_in = pixels_in.flatten()\n    pixels_out = [0] * len(pixels_in)\n    \n    for i in range(0, len(pixels_in)):\n        pixel = pixels_in[i]\n        pixels_out[i] = int(lut[pixel])\n        \n    return pixels_out","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:48.850311Z","iopub.execute_input":"2021-09-13T02:16:48.850734Z","iopub.status.idle":"2021-09-13T02:16:48.863353Z","shell.execute_reply.started":"2021-09-13T02:16:48.850702Z","shell.execute_reply":"2021-09-13T02:16:48.862077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### - Load and display an image","metadata":{}},{"cell_type":"code","source":"from pydicom.pixel_data_handlers.util import apply_voi_lut\n# Load an image\nimage = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00148/T1wCE/Image-95.dcm')\npixels = image.pixel_array\n\n# Print out the pixel 'width'\nprint(\"Min pixel value: \" + str(np.min(pixels)))\nprint(\"Max pixel value: \" + str(np.max(pixels)))\n\nprint(image.PhotometricInterpretation)\n\n\nplt.figure(figsize= (6,6))\nplt.imshow(pixels, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:29:18.359564Z","iopub.execute_input":"2021-09-13T02:29:18.360019Z","iopub.status.idle":"2021-09-13T02:29:18.5954Z","shell.execute_reply.started":"2021-09-13T02:29:18.359983Z","shell.execute_reply":"2021-09-13T02:29:18.594341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot a histogram of the raw pixel data\nfig, axes = plt.subplots(nrows=1, ncols=1,sharex=False, sharey=False, figsize=(10,4))\nplt.title('Pixel Range: ' + str(np.min(pixels)) + '-' + str(np.max(pixels)))\nplt.hist(pixels.ravel(), np.max(pixels), (1, np.max(pixels)))\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:49.127664Z","iopub.execute_input":"2021-09-13T02:16:49.128014Z","iopub.status.idle":"2021-09-13T02:16:51.985516Z","shell.execute_reply.started":"2021-09-13T02:16:49.127981Z","shell.execute_reply":"2021-09-13T02:16:51.984114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- The pixel range or width, is from 0 to 2857. That's 2858 shades of gray.\n- Since we can't display such a 'wide' image on our consumer grade display adapters and monitors, we need to bin the pixels to a usable range (8 bit).\n\n### - Make a LUT\n- If we specify **window_center = window_width / 2** as the center, the image will appear as the default image above does.\n- Lowering the width has the effect of *decreasing the contrast* by binning pixels.\n- Lowering the center has the effect of *decreasing the brightness*.\n\n### - Now we'll tweak the window values to produce a more contrasty image.","metadata":{}},{"cell_type":"code","source":"# Apply three different WW/WL settings via LUT. We'll set the center slightly less than half to adjust for brightness.\nwindow_width_1 = np.max(image.pixel_array)\nwindow_center_1 = window_width_1 / 2\n\nlut = make_lut(image.pixel_array, window_width_1, window_center_1, image.PhotometricInterpretation)\nimage1 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0],pixels.shape[1]))\n\nwindow_width_2 = 450\nwindow_center_2 = 450\n\nlut = make_lut(image.pixel_array, window_width_2, window_center_2, image.PhotometricInterpretation)\nimage2 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0],pixels.shape[1]))\n\nwindow_width_3 = 900\nwindow_center_3 = 90\n\nlut = make_lut(image.pixel_array, window_width_3, window_center_3, image.PhotometricInterpretation)\nimage3 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0],pixels.shape[1]))","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:51.987239Z","iopub.execute_input":"2021-09-13T02:16:51.987636Z","iopub.status.idle":"2021-09-13T02:16:52.593019Z","shell.execute_reply.started":"2021-09-13T02:16:51.987601Z","shell.execute_reply":"2021-09-13T02:16:52.591665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2, ncols=2,sharex=True, sharey=True, figsize=(12, 12))\nax = axes.ravel()\nax[0].set_title('Default Image')\nax[0].imshow(image.pixel_array, cmap='gray')\nax[1].set_title(f'Width: {window_width_1} / Center: {window_center_1}')\nax[1].imshow(image1, cmap='gray')\nax[2].set_title(f'Width: {window_width_2} / Center: {window_center_2}')\nax[2].imshow(image2, cmap='gray')\nax[3].set_title(f'Width: {window_width_3} / Center: {window_center_3}')\nax[3].imshow(image3, cmap='gray')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:52.594773Z","iopub.execute_input":"2021-09-13T02:16:52.595125Z","iopub.status.idle":"2021-09-13T02:16:53.382205Z","shell.execute_reply.started":"2021-09-13T02:16:52.595088Z","shell.execute_reply":"2021-09-13T02:16:53.38084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### - Result: \n#### - The tumor area is more contrasted against the surrounding brain tissue. Especially in the bottom left image.\n#### - There are fine details visible in this image that are obscured on the other images.\n\n### - Let's load another image and try the same process.\n- Just copy the last 3 cells and change the filename.","metadata":{}},{"cell_type":"code","source":"# Load an image\nimage = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00014/FLAIR/Image-126.dcm')\npixels = image.pixel_array\n\n# Print out the pixel 'width'\nprint(\"Min pixel value: \" + str(np.min(pixels)))\nprint(\"Max pixel value: \" + str(np.max(pixels)))\n\nplt.figure(figsize= (6,6))\nplt.imshow(pixels, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:53.383815Z","iopub.execute_input":"2021-09-13T02:16:53.38416Z","iopub.status.idle":"2021-09-13T02:16:53.623813Z","shell.execute_reply.started":"2021-09-13T02:16:53.384125Z","shell.execute_reply":"2021-09-13T02:16:53.622753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### - This image doesn't have as wide of pixel range as the first image does. \n- We'll set it's level a little bit lower and try to make the tumor stand out better.","metadata":{}},{"cell_type":"code","source":"# Apply three different WW/WL settings via LUT. We'll set the level slightly less than half to adjust for brightness.\nwindow_width_1 = 1000\nwindow_center_1 = 900\n\nlut = make_lut(image.pixel_array, window_width_1, window_center_1, image.PhotometricInterpretation)\nimage1 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0],pixels.shape[1]))\n\nwindow_width_2 = 600\nwindow_width_2 = 900\n\nlut = make_lut(image.pixel_array, window_width_2, window_center_2, image.PhotometricInterpretation)\nimage2 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0],pixels.shape[1]))\n\nwindow_width_3 = 300\nwindow_center_3 = 900\n\nlut = make_lut(image.pixel_array, window_width_3, window_center_3, image.PhotometricInterpretation)\nimage3 = np.reshape(apply_lut(pixels, lut), (pixels.shape[0],pixels.shape[1]))","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:53.62533Z","iopub.execute_input":"2021-09-13T02:16:53.625807Z","iopub.status.idle":"2021-09-13T02:16:54.231143Z","shell.execute_reply.started":"2021-09-13T02:16:53.625753Z","shell.execute_reply":"2021-09-13T02:16:54.229973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2, ncols=2,sharex=True, sharey=True, figsize=(12, 12))\nax = axes.ravel()\nax[0].set_title('Default Image')\nax[0].imshow(image.pixel_array, cmap='gray')\nax[1].set_title(f'Width: {window_width_1} / Center: {window_center_1}')\nax[1].imshow(image1, cmap='gray')\nax[2].set_title(f'Width: {window_width_2} / Center: {window_center_2}')\nax[2].imshow(image2, cmap='gray')\nax[3].set_title(f'Width: {window_width_3} / Center: {window_center_3}')\nax[3].imshow(image3, cmap='gray')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-13T02:16:54.233503Z","iopub.execute_input":"2021-09-13T02:16:54.233861Z","iopub.status.idle":"2021-09-13T02:16:55.274128Z","shell.execute_reply.started":"2021-09-13T02:16:54.233811Z","shell.execute_reply":"2021-09-13T02:16:55.273048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Again, we can clearly see the tumor better. I suspect CNNs will too!\n\n#### Some of my other MR notebooks\n- Tumor Object Detection -> https://www.kaggle.com/davidbroberts/brain-tumor-object-detection\n- Determining MR image planes -> https://www.kaggle.com/davidbroberts/determining-mr-image-planes\n- Determining MR Slice Orientation -> https://www.kaggle.com/davidbroberts/determining-mr-slice-orientation\n- Determining DICOM image order -> https://www.kaggle.com/davidbroberts/determining-dicom-image-order\n- Reference Lines on MR images -> https://www.kaggle.com/davidbroberts/mr-reference-lines\n- Standardizing MR Images -> https://www.kaggle.com/davidbroberts/standardizing-mr-images\n- Export DICOM Images by Plane -> https://www.kaggle.com/davidbroberts/export-dicom-series-by-plane","metadata":{}}]}