{"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>Applying filters to Chest X-Rays</h1>- yet another chest x-ray processing notebook -</div>\n\n#### While exploring methods to better visualize lung tissue, I found some useful filters from sklearn.filters I thought I'd share.\n#### In this notebook we'll grab an image and apply a filter or two. Params can be tweaked to get some pretty unique views.\n\n- In this case, we'll use two modules from *sklearn.exposure* .. **equalize_hist** and **equalize_adapthist**\n- The filters are **unsharp_mask**, **meijering**, **sato**, **scharr**, and **hessian** from *sklearn.filters*","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-07-04T21:35:09.7136Z","iopub.execute_input":"2021-07-04T21:35:09.71427Z"}}},{"cell_type":"code","source":"# You will need to run this command if opening a JPG encoded DICOM file\n#!conda install gdcm -c conda-forge -y","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:05.52703Z","iopub.execute_input":"2021-07-05T16:25:05.527537Z","iopub.status.idle":"2021-07-05T16:25:05.530834Z","shell.execute_reply.started":"2021-07-05T16:25:05.527453Z","shell.execute_reply":"2021-07-05T16:25:05.529997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load stuff 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\nfrom skimage.filters import unsharp_mask, meijering, sato, scharr, hessian\nfrom skimage.exposure import exposure, equalize_hist, equalize_adapthist\nimport cv2 as cv","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:05.532075Z","iopub.execute_input":"2021-07-05T16:25:05.532414Z","iopub.status.idle":"2021-07-05T16:25:08.209674Z","shell.execute_reply.started":"2021-07-05T16:25:05.532389Z","shell.execute_reply":"2021-07-05T16:25:08.207886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This function gets the first image path in a StudyInstanceUID directory in the train set\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-07-05T16:25:08.212664Z","iopub.execute_input":"2021-07-05T16:25:08.213061Z","iopub.status.idle":"2021-07-05T16:25:08.219208Z","shell.execute_reply.started":"2021-07-05T16:25:08.213024Z","shell.execute_reply":"2021-07-05T16:25:08.218282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load a DICOM file and get the pixels\ndef load_image(study_id):\n    img_file = get_image_by_study_id(study_id)\n    image = pydicom.dcmread(img_file)\n    pixels = image.pixel_array\n\n    min_pixel = np.min(pixels)\n    max_pixel = np.max(pixels)\n\n    if image.PhotometricInterpretation == \"MONOCHROME1\":\n        pixels = max_pixel - pixels\n    else:\n        pixels = pixels\n\n    return pixels","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:08.220358Z","iopub.execute_input":"2021-07-05T16:25:08.220634Z","iopub.status.idle":"2021-07-05T16:25:08.238269Z","shell.execute_reply.started":"2021-07-05T16:25:08.22061Z","shell.execute_reply":"2021-07-05T16:25:08.237377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply filters. Tweak params here.\ndef apply_filter(f, img):\n    if f == 'equalize_hist':\n        img = equalize_hist(img, nbins=256, mask=None)\n        \n    if f == 'equalize_adapthist': \n        img = equalize_adapthist(img, kernel_size=None, clip_limit=0.01, nbins=256)\n        \n    if f == 'unsharp_mask':\n        img = unsharp_mask(img, radius=5, amount=2)\n        \n    if f == 'meijering':\n        img = meijering(img, sigmas=range(1, 10, 2), alpha=None, black_ridges=True, mode='reflect', cval=0)\n        \n    if f == 'sato':\n        img = sato(img, sigmas=range(1, 10, 2), black_ridges=True, mode='reflect', cval=0)\n        \n    if f == 'scharr':\n        img = scharr(img, mask=None, axis=None, mode='reflect', cval=0.0)\n        \n    if f == 'hessian':\n        img = hessian(img, sigmas=range(1, 10, 2), scale_range=None, scale_step=None, alpha=0.5, beta=0.5, gamma=15, black_ridges=True, mode='reflect', cval=0)\n    \n    if f == 'threshold_isodata':\n        img = threshold_isodata(image=img, nbins=256, return_all=False, hist=None)\n    \n    return img","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:08.239214Z","iopub.execute_input":"2021-07-05T16:25:08.239526Z","iopub.status.idle":"2021-07-05T16:25:08.262635Z","shell.execute_reply.started":"2021-07-05T16:25:08.239492Z","shell.execute_reply":"2021-07-05T16:25:08.261315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Display the original image and the processed image\ndef plot_images(title, image, image_processed):\n    # Plot both images\n    fig, axes = plt.subplots(nrows=1, ncols=2,sharex=True, sharey=True, figsize=(12, 12))\n    ax = axes.ravel()\n    ax[0].imshow(image, cmap=plt.cm.gray)\n    ax[0].set_title('Original image')\n    ax[1].imshow(image_processed, cmap=plt.cm.gray)\n    ax[1].set_title(title)\n    for a in ax:\n        a.axis('off')\n    fig.tight_layout()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:08.263926Z","iopub.execute_input":"2021-07-05T16:25:08.264186Z","iopub.status.idle":"2021-07-05T16:25:08.288164Z","shell.execute_reply.started":"2021-07-05T16:25:08.26416Z","shell.execute_reply":"2021-07-05T16:25:08.287023Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load an image and apply equalization\n- We'll use Histogram Equalization and Contrast Limited Adaptive Histogram Equalization (CLAHE)\n- Hist EQ gives a nice contrasty view of a chest x-ray.\n- CLAHE seems to be more balanced and 'wide' if the input image isn't poorly balanced.","metadata":{}},{"cell_type":"code","source":"# Load an image using the Study_ID\nimage = load_image('013d698aeecb')","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:08.289575Z","iopub.execute_input":"2021-07-05T16:25:08.289945Z","iopub.status.idle":"2021-07-05T16:25:09.035494Z","shell.execute_reply.started":"2021-07-05T16:25:08.289909Z","shell.execute_reply":"2021-07-05T16:25:09.034438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply Histogram Equalization\nimage_out = apply_filter('equalize_hist', image)\nplot_images('Histogram Equalization',image, image_out)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:09.036627Z","iopub.execute_input":"2021-07-05T16:25:09.037058Z","iopub.status.idle":"2021-07-05T16:25:11.315165Z","shell.execute_reply.started":"2021-07-05T16:25:09.037021Z","shell.execute_reply":"2021-07-05T16:25:11.314299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply Adaptive Histogram Equalization (CLAHE)\nimage_out = apply_filter('equalize_adapthist', image)\nplot_images('CLAHE',image, image_out)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:11.31686Z","iopub.execute_input":"2021-07-05T16:25:11.317227Z","iopub.status.idle":"2021-07-05T16:25:13.909093Z","shell.execute_reply.started":"2021-07-05T16:25:11.317198Z","shell.execute_reply":"2021-07-05T16:25:13.908087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Put an unsharp mask on the original image\n- Why is it called unsharp mask when it actually sharpens the image? (because the first step of unsharp is to blur, which is .. um .. unsharp)\n- The unsharpened image is more crisp and defined.","metadata":{}},{"cell_type":"code","source":"# Apply unsharp mask\nimage_out = apply_filter('unsharp_mask', image_out)\nplot_images('Unsharp',image, image_out)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:13.910331Z","iopub.execute_input":"2021-07-05T16:25:13.910568Z","iopub.status.idle":"2021-07-05T16:25:15.891216Z","shell.execute_reply.started":"2021-07-05T16:25:13.910541Z","shell.execute_reply":"2021-07-05T16:25:15.889847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sklearn has a lot of edge detection filters. Here's a few I found useful.\n- Play around with the default params in the apply_filter method to get some interesting results.","metadata":{}},{"cell_type":"code","source":"# Apply Meijering and then Hist EQ to brighten it up\nimage_out = apply_filter('meijering', image)\nimage_out = apply_filter('equalize_hist', image_out)\nplot_images('Meijering',image, image_out)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:15.892656Z","iopub.execute_input":"2021-07-05T16:25:15.893121Z","iopub.status.idle":"2021-07-05T16:25:26.793951Z","shell.execute_reply.started":"2021-07-05T16:25:15.893083Z","shell.execute_reply":"2021-07-05T16:25:26.793199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply Sato and Hist EQ\nimage_out = apply_filter('sato', image)\nimage_out = apply_filter('equalize_hist', image_out)\nplot_images('Sato',image, image_out)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:26.795081Z","iopub.execute_input":"2021-07-05T16:25:26.795476Z","iopub.status.idle":"2021-07-05T16:25:36.445951Z","shell.execute_reply.started":"2021-07-05T16:25:26.795431Z","shell.execute_reply":"2021-07-05T16:25:36.445011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Scharr is similar to Sobel, maybe a little faster?\nimage_out = apply_filter('scharr', image)\nimage_out = apply_filter('equalize_hist', image_out)\nplot_images('Scharr',image, image_out)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:36.447099Z","iopub.execute_input":"2021-07-05T16:25:36.44738Z","iopub.status.idle":"2021-07-05T16:25:39.035076Z","shell.execute_reply.started":"2021-07-05T16:25:36.44735Z","shell.execute_reply":"2021-07-05T16:25:39.034138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply Hessian and CLAHE\nimage_out = apply_filter('hessian', image)\nimage_out = apply_filter('equalize_adapthist', image_out)\nplot_images('Hessian',image, image_out)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T16:25:39.03616Z","iopub.execute_input":"2021-07-05T16:25:39.036432Z","iopub.status.idle":"2021-07-05T16:25:55.043662Z","shell.execute_reply.started":"2021-07-05T16:25:39.036403Z","shell.execute_reply":"2021-07-05T16:25:55.042812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### As you can see, there *are* methods to demonstrate tissue areas with more/less contrast. Some are noisy. Some are not.\n\n- Stacking filters and using various parameter settings can yield interesting results.\n\nHere are some other processing notebooks I made:\n- Lung Segmentation Without CNN -> https://www.kaggle.com/davidbroberts/lung-segmentation-without-cnn\n- Rib supression on Chest X-Rays -> https://www.kaggle.com/davidbroberts/rib-suppression-poc\n- Manual DICOM VOI LUT -> https://www.kaggle.com/davidbroberts/manual-dicom-voi-lut\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- Bounding Boxes on Cropped Images -> https://www.kaggle.com/davidbroberts/bounding-boxes-on-cropped-images\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\n- Standardizins Chest X-Ray Dataset Exports -> https://www.kaggle.com/davidbroberts/standardizing-cxr-datasets","metadata":{}}]}