{"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>Lung segmentation without a CNN</h1>\n    - yet another chest x-ray processing notebook -\n</div>\n\n#### In this notebook, we'll use a combination of sklearn processing and filters to extract lung areas and bounding boxes.\n\nThe basic steps are:\n- Equalize the image with CLAHE\n- Create a threshold mask to separate tissue by pixel intensity\n- Find regions in the threshold\n- Remove borders\n- Fill small holes\n- Extract lung areas\n\nThe end result is a segment mask and a bounding boxes that represent the lung fields.","metadata":{}},{"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-05T21:40:53.603971Z","iopub.execute_input":"2021-07-05T21:40:53.604292Z","iopub.status.idle":"2021-07-05T21:40:53.607811Z","shell.execute_reply.started":"2021-07-05T21:40:53.604266Z","shell.execute_reply":"2021-07-05T21:40:53.606452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pydicom\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nfrom skimage.filters import threshold_otsu\nfrom skimage.exposure import equalize_adapthist\nfrom skimage.color import label2rgb\nfrom skimage.segmentation import clear_border\nfrom skimage.measure import label, regionprops\nfrom skimage.morphology import square,closing","metadata":{"execution":{"iopub.status.busy":"2021-07-05T21:40:53.94718Z","iopub.execute_input":"2021-07-05T21:40:53.947486Z","iopub.status.idle":"2021-07-05T21:40:53.952492Z","shell.execute_reply.started":"2021-07-05T21:40:53.94746Z","shell.execute_reply":"2021-07-05T21:40:53.951481Z"},"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.\n# I use this so I don't have to keep clicking down to get a test image path each time I try a new one.\n\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-05T21:40:53.990367Z","iopub.execute_input":"2021-07-05T21:40:53.990699Z","iopub.status.idle":"2021-07-05T21:40:53.996045Z","shell.execute_reply.started":"2021-07-05T21:40:53.990668Z","shell.execute_reply":"2021-07-05T21:40:53.994969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load a DICOM file and get the pixels\n\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-05T21:40:53.997818Z","iopub.execute_input":"2021-07-05T21:40:53.998092Z","iopub.status.idle":"2021-07-05T21:40:54.007305Z","shell.execute_reply.started":"2021-07-05T21:40:53.998066Z","shell.execute_reply":"2021-07-05T21:40:54.006167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This is the main function, it applies all the filters and plots the original and segmented image\n# Tweak params here to fine tune\n\ndef segment(img):\n\n    # Even out the contrast with CLAHE\n    img = equalize_adapthist(img, kernel_size=None, clip_limit=0.01, nbins=256)\n    \n    # Make a binary threshold mask and apply it to the image \n    thresh = threshold_otsu(image=img, nbins=256, hist=None)\n    thresh = img > thresh\n    bw = closing(img > thresh, square(3))\n\n    # clean up the borders\n    cleared = clear_border(bw)\n\n    # label image regions\n    label_image = label(cleared)\n    image_label_overlay = label2rgb(label_image, image=img, bg_label=0)\n    \n    # Plot the 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[1].imshow(image_label_overlay)\n    \n    # Iterate through the regions\n    for region in regionprops(label_image):\n        \n        # Only get large regions, 250,000 is a reasonably large enough area for lung fields, but this might take some tweaking\n        if region.area >= 250000:\n            \n            # draw a box around segments\n            minr, minc, maxr, maxc = region.bbox\n            rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr, fill=False, edgecolor='purple', linewidth=2)\n            ax[1].add_patch(rect)\n            print(rect)\n\n    ax[1].set_axis_off()\n    plt.tight_layout()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-07-05T21:40:54.009154Z","iopub.execute_input":"2021-07-05T21:40:54.009417Z","iopub.status.idle":"2021-07-05T21:40:54.020814Z","shell.execute_reply.started":"2021-07-05T21:40:54.00939Z","shell.execute_reply":"2021-07-05T21:40:54.019941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load and segment an image\n\n- I'm just randomly picking images from the train set.","metadata":{}},{"cell_type":"code","source":"# Use the studyID to get an image\nimage = load_image('00292f8c37bd')\nsegment(image)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T21:40:54.02186Z","iopub.execute_input":"2021-07-05T21:40:54.022136Z","iopub.status.idle":"2021-07-05T21:40:59.902882Z","shell.execute_reply.started":"2021-07-05T21:40:54.022108Z","shell.execute_reply":"2021-07-05T21:40:59.900718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- The results aren't bad. The colored lung pixel areas could be 'grown', or made into convex hull .. something to smooth them out.\n\n#### Let's look at a few more images.","metadata":{}},{"cell_type":"code","source":"image = load_image('013d698aeecb')\nsegment(image)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T21:40:59.905965Z","iopub.execute_input":"2021-07-05T21:40:59.906475Z","iopub.status.idle":"2021-07-05T21:41:07.110471Z","shell.execute_reply.started":"2021-07-05T21:40:59.906436Z","shell.execute_reply":"2021-07-05T21:41:07.109495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image = load_image('00086460a852')\nsegment(image)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T21:41:07.111631Z","iopub.execute_input":"2021-07-05T21:41:07.111874Z","iopub.status.idle":"2021-07-05T21:41:12.47641Z","shell.execute_reply.started":"2021-07-05T21:41:07.111851Z","shell.execute_reply":"2021-07-05T21:41:12.475574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image = load_image('00c74279c5b7')\nsegment(image)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T21:41:12.477455Z","iopub.execute_input":"2021-07-05T21:41:12.477743Z","iopub.status.idle":"2021-07-05T21:41:18.821146Z","shell.execute_reply.started":"2021-07-05T21:41:12.477715Z","shell.execute_reply":"2021-07-05T21:41:18.820262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image = load_image('015a2029ad0c')\nsegment(image)","metadata":{"execution":{"iopub.status.busy":"2021-07-05T21:41:18.822326Z","iopub.execute_input":"2021-07-05T21:41:18.822608Z","iopub.status.idle":"2021-07-05T21:41:25.623711Z","shell.execute_reply.started":"2021-07-05T21:41:18.822579Z","shell.execute_reply":"2021-07-05T21:41:25.622722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Conclusion\n\n- Bounding box data and segment maps can be extracted and used for further processing/modeling.\n- Using thresholding to segment images might actually be useful in finding lung tissue in chest radiographs due to the relatively contrasty periphery.\n- Adjustments to the input image's contrast might result in even better segmentation.\n- This is reasonably accurate on images that aren't diagnostically crappy.","metadata":{}},{"cell_type":"markdown","source":"**Here are some other processing notebooks I made:**\n- Applying filters to x-rays -> https://www.kaggle.com/davidbroberts/applying-filters-to-chest-x-rays\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- Standardizing Chest X-Ray Dataset Exports -> https://www.kaggle.com/davidbroberts/standardizing-cxr-datasets","metadata":{}}]}