{"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>Standardizing MR Images</h1>\n- yet another MR processing notebook -</div>\n\n#### The goal of this notebook is to demonstrate how to normalize a DICOM MR image with minimal loss of information.\n#### - We'll load and normalize an MR image and attempt to maintain it's original pixel distribution.","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport matplotlib.pyplot as plt\nfrom skimage import exposure\nimport cv2\nfrom scipy import stats\nfrom IPython.display import clear_output","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:26.266827Z","iopub.execute_input":"2021-09-03T14:53:26.267222Z","iopub.status.idle":"2021-09-03T14:53:26.272952Z","shell.execute_reply.started":"2021-09-03T14:53:26.267192Z","shell.execute_reply":"2021-09-03T14:53:26.272151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### - Define a LUT function","metadata":{}},{"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-03T14:53:26.280121Z","iopub.execute_input":"2021-09-03T14:53:26.280731Z","iopub.status.idle":"2021-09-03T14:53:26.292549Z","shell.execute_reply.started":"2021-09-03T14:53:26.280695Z","shell.execute_reply":"2021-09-03T14:53:26.291275Z"},"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 = pixels_in.flatten()\n    pixels_out = [0] * len(pixels)\n    for i in range(0, len(pixels)):\n        pixel = pixels[i]\n        if pixel > 0:\n            pixels_out[i] = int(lut[pixel])\n    return np.reshape(pixels_out, (pixels_in.shape[0],pixels_in.shape[1]))","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:26.294899Z","iopub.execute_input":"2021-09-03T14:53:26.295769Z","iopub.status.idle":"2021-09-03T14:53:26.307479Z","shell.execute_reply.started":"2021-09-03T14:53:26.295715Z","shell.execute_reply":"2021-09-03T14:53:26.306372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### - Load an image","metadata":{}},{"cell_type":"code","source":"# Load an image\nimage = pydicom.dcmread('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00216/T1wCE/Image-98.dcm')\nraw_pixels = image.pixel_array\n\n# Plot the image\nplt.figure(figsize= (6,6))\nplt.imshow(raw_pixels, cmap='gray');","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:26.309604Z","iopub.execute_input":"2021-09-03T14:53:26.310339Z","iopub.status.idle":"2021-09-03T14:53:26.560479Z","shell.execute_reply.started":"2021-09-03T14:53:26.310292Z","shell.execute_reply":"2021-09-03T14:53:26.559302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot a histogram of the raw pixels\nfig, axes = plt.subplots(nrows=1, ncols=1,sharex=False, sharey=False, figsize=(10,4))\nplt.title('Pixel Range of raw pixels: ' + str(np.min(raw_pixels)) + '-' + str(np.max(raw_pixels)))\nplt.hist(raw_pixels.ravel(), np.max(raw_pixels), (1, np.max(raw_pixels)))\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:26.562483Z","iopub.execute_input":"2021-09-03T14:53:26.562801Z","iopub.status.idle":"2021-09-03T14:53:28.395116Z","shell.execute_reply.started":"2021-09-03T14:53:26.562771Z","shell.execute_reply":"2021-09-03T14:53:28.394184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- The raw pixels have a range of 0 - 737. This image has more than 8 bits of data.\n- We can't display a 16 bit image, so we have to bin the pixels .. either with a LUT, or through normalization.\n\n#### - Let's normalize this image five different ways.","metadata":{}},{"cell_type":"markdown","source":"#### 1. Create an 'Auto-LUT' image\n- Use the raw pixel range to calculate a 'centered' image.\n- This is what pydicom and matplotlib do.","metadata":{}},{"cell_type":"code","source":"# Calculate the width and center of the pixels to make a LUT\nauto_lut_window_width = np.max(raw_pixels)\nauto_lut_window_center = (np.max(raw_pixels) - np.min(raw_pixels)) / 2 + np.min(raw_pixels)\n\nlut = make_lut(raw_pixels, auto_lut_window_width, auto_lut_window_center, image.PhotometricInterpretation)\nimage_autolut = apply_lut(raw_pixels, lut)","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:28.396898Z","iopub.execute_input":"2021-09-03T14:53:28.397303Z","iopub.status.idle":"2021-09-03T14:53:28.620522Z","shell.execute_reply.started":"2021-09-03T14:53:28.397268Z","shell.execute_reply":"2021-09-03T14:53:28.619302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 2. Use default WindowWidth and WindowCenter DICOM tags\n- This is how the modality unit (MR machine) thinks the image should be displayed according to a predefined protocol.\n- The problem with these is that they're not consistent and some appear to have been set prior to bone stripping. So they're just wrong in some cases. But they are great in others.","metadata":{}},{"cell_type":"code","source":"# Get the DICOM tags for ww/wc and use them to create the LUT\nif image.WindowWidth != '' and image.WindowCenter != '':\n    window_width = image.WindowWidth\n    window_center = image.WindowCenter\n\nlut = make_lut(raw_pixels, window_width, window_center, image.PhotometricInterpretation)\nimage_default_window = apply_lut(raw_pixels, lut)","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:28.624278Z","iopub.execute_input":"2021-09-03T14:53:28.624979Z","iopub.status.idle":"2021-09-03T14:53:28.834657Z","shell.execute_reply.started":"2021-09-03T14:53:28.624908Z","shell.execute_reply":"2021-09-03T14:53:28.833337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3. Normalize with code\n- This is how we routinely do it, a simple linear transform.","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:17:29.854354Z","iopub.execute_input":"2021-09-03T14:17:29.854939Z","iopub.status.idle":"2021-09-03T14:17:29.867745Z","shell.execute_reply.started":"2021-09-03T14:17:29.854899Z","shell.execute_reply":"2021-09-03T14:17:29.866968Z"}}},{"cell_type":"code","source":"pixels = raw_pixels - np.min(raw_pixels)\npixels = pixels / np.max(pixels)\nimage_manual_norm = (pixels * 255).astype(np.uint8)","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:28.836116Z","iopub.execute_input":"2021-09-03T14:53:28.836417Z","iopub.status.idle":"2021-09-03T14:53:28.842111Z","shell.execute_reply.started":"2021-09-03T14:53:28.836389Z","shell.execute_reply":"2021-09-03T14:53:28.841071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 4. Apply CLAHE\n- Use CLAHE with default settings on the raw pixels","metadata":{}},{"cell_type":"code","source":"image_clahe = exposure.equalize_adapthist(raw_pixels)","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:28.843616Z","iopub.execute_input":"2021-09-03T14:53:28.843904Z","iopub.status.idle":"2021-09-03T14:53:28.870149Z","shell.execute_reply.started":"2021-09-03T14:53:28.843877Z","shell.execute_reply":"2021-09-03T14:53:28.868781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 5. Apply a manual VOI LUT\n- Make the image wider and maintain contrast.","metadata":{}},{"cell_type":"code","source":"lut_window_width = 350\nlut_window_center = 500\n\nlut = make_lut(raw_pixels, lut_window_width, lut_window_center, image.PhotometricInterpretation)\nimage_lut = apply_lut(raw_pixels, lut)","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:28.871191Z","iopub.execute_input":"2021-09-03T14:53:28.871481Z","iopub.status.idle":"2021-09-03T14:53:29.07788Z","shell.execute_reply.started":"2021-09-03T14:53:28.871453Z","shell.execute_reply":"2021-09-03T14:53:29.076656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### - Display the images and histograms after normalization","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=5, ncols=2,sharex=False, sharey=False, figsize=(16, 16))\nax = axes.ravel()\n\nax[0].set_title('Default WW/WC ' + str(image.WindowWidth) + \"/\" + str(image.WindowCenter))\nax[0].imshow(image_default_window, cmap='gray')\n\nax[1].set_title('Pixel Range: ' + str(np.min(image_default_window)) + '-' + str(np.max(image_default_window)))\nax[1].hist(image_default_window.ravel(), np.max(image_default_window), (1, np.max(image_default_window)))\n\nax[2].set_title('Manual Norm ')\nax[2].imshow(image_manual_norm, cmap='gray')\n\nax[3].set_title('Pixel Range: ' + str(np.min(image_manual_norm)) + '-' + str(np.max(image_manual_norm)))\nax[3].hist(image_manual_norm.ravel(), np.max(image_manual_norm), (1, np.max(image_manual_norm)))\n\nax[4].set_title('Auto-LUT ' + str(auto_lut_window_width) + '/' + str(auto_lut_window_center))\nax[4].imshow(image_autolut, cmap='gray')\n\nax[5].set_title('Pixel Range: ' + str(np.min(image_autolut)) + '-' + str(np.max(image_autolut)))\nax[5].hist(image_autolut.ravel(), 254, (1, 255))\n\nax[6].set_title('CLAHE ')\nax[6].imshow(image_clahe, cmap='gray')\n\nax[7].set_title('Pixel Range: ' + str(np.min(image_clahe)) + '-' + str(np.max(image_clahe)))\nax[7].hist(image_clahe.ravel(), 254, (0.1,1))\n\n\nax[8].set_title('Manual LUT: ' + str(lut_window_width) + '/' + str(lut_window_center))\nax[8].imshow(image_lut, cmap='gray')\n\nax[9].set_title('Pixel Range: ' + str(np.min(image_lut)) + '-' + str(np.max(image_lut)))\nax[9].hist(image_lut.ravel(), 254, (1, 254))\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:29.079656Z","iopub.execute_input":"2021-09-03T14:53:29.079936Z","iopub.status.idle":"2021-09-03T14:53:32.974714Z","shell.execute_reply.started":"2021-09-03T14:53:29.079909Z","shell.execute_reply":"2021-09-03T14:53:32.973345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- The raw pixels have values within the 0 to 737 range. That's 738 shades of gray\n- In the rest of the images, the same distribution shape occurs, but the pixel value ranges are reduced to 0-255 range. 256 shades of gray (this is Lossy)\n- Notice the Default WW/WC image is a little brighter .. it's difficult to see on the picture, but the hist is shifted toward the right.\n- Manual Normalization and Auto-LUT result in the same image (maybe need some rounding or slight tweaks, but it's almost identical).\n\n#### IMPORTANT !!\n- You can see the tumor in the histograms, it's the bump on the left .. a large group of darker pixels that are a different shade then most of the other pix.\n- CLAHE removes this by equalizing contrast between nearby pixels, making the soft tumor tissue look more like the normal tissue .. exactly what we don't want.\n- By applying a VOI LUT, we've 'widened' the image (more shades of gray), but not lost the contrast. The distribution is basically the same shape, just fatter.\n- Notice the bump on the left of the histogram. We didn't change the tumor to match the brain. You can see, the tumor stands out a little more in the LUT image.\n\nIt's hard to visualize, let's enlarge the images a little bit.","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=1, ncols=3,sharex=False, sharey=False, figsize=(16, 10))\nax = axes.ravel()\n\nax[0].set_title('Norm')\nax[0].imshow(image_manual_norm, cmap='gray')\n\nax[1].set_title('CLAHE')\nax[1].imshow(image_clahe, cmap='gray')\n\nax[2].set_title('LUT')\nax[2].imshow(image_lut, cmap='gray')\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-03T14:53:32.97625Z","iopub.execute_input":"2021-09-03T14:53:32.976636Z","iopub.status.idle":"2021-09-03T14:53:33.528903Z","shell.execute_reply.started":"2021-09-03T14:53:32.976539Z","shell.execute_reply":"2021-09-03T14:53:33.525206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Conclusion:\n\n- Any normalization techniques applied to MR images must maintain the pixel distribution of the raw pixel data.\n- Applying equalization or filters that disturb the inherent contrast of a 16 bit image is counter-productive.","metadata":{}},{"cell_type":"markdown","source":"- 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- Manual VOI LUT on MR images -> https://www.kaggle.com/davidbroberts/manual-voi-lut-on-mr-images\n- Export DICOM Images by Plane -> https://www.kaggle.com/davidbroberts/export-dicom-series-by-plane","metadata":{}}]}