{"cells":[{"metadata":{},"cell_type":"markdown","source":"## This notebooks shows how to use Skimage method to segment Lung part only\n* Fast way to segment lung part only in CT by using Skimage"},{"metadata":{},"cell_type":"markdown","source":"### Load Module"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\n\nimport pydicom as dcm\nfrom pydicom.pixel_data_handlers.util import apply_modality_lut","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Get Segmented Lungs\n#### Step:\n* Covert into binary image (theshold set HU -400), [HU threshold information](https://en.wikipedia.org/wiki/Hounsfield_scale)\n* Remove the blobs connect to the border of the image\n* label the image\n* Keep the all labels.\n* Erosion operation with a disk of radius 2. This operation is seperate the lung nodules attached to the blood vessels.\n* Closure operation with a disk of radius 10. This operation is to keep nodules attached to the lung wall.\n* Fill in the small holes inside the binary mask of lungs.\n* Superimpose the binary mask on the input image."},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfrom skimage.segmentation import clear_border\nfrom skimage.measure import label, regionprops\nfrom skimage.morphology import disk, dilation, binary_erosion, binary_closing\nfrom skimage.filters import roberts, sobel\nimport cv2\nfrom scipy import ndimage as ndi\n\ndef get_segmented_lungs(im2, plot=False):\n    im = im2.copy()\n    # Step 1: Convert into a binary image.\n    binary = im < -400\n    \n    if plot:\n        plt.imshow(binary)\n        plt.show()\n        \n    # Step 2: Remove the blobs connected to the border of the image.\n    cleared = clear_border(binary)\n    \n    if plot:\n        plt.imshow(cleared)\n        plt.show()    \n        \n    # Step 3: Label the image.\n    label_image = label(cleared)\n    \n    if plot:\n        plt.imshow(label_image)\n        plt.show()    \n        \n    # Step 4: Keep the labels with 2 largest areas.\n    areas = [r.area for r in regionprops(label_image)]\n    areas.sort()\n    if len(areas) > 0:\n        for region in regionprops(label_image):\n            if region.area < areas[0]:\n                for coordinates in region.coords:\n                       label_image[coordinates[0], coordinates[1]] = 0\n    binary = label_image > 0\n    \n    if plot:\n        plt.imshow(binary)\n        plt.show()  \n        \n    # Step 5: Erosion operation with a disk of radius 2. This operation is seperate the lung nodules attached to the blood vessels.\n    selem = disk(2)\n    binary = binary_erosion(binary, selem)\n    \n    if plot:\n        plt.imshow(binary)\n        plt.show()  \n        \n    # Step 6: Closure operation with a disk of radius 10. This operation is to keep nodules attached to the lung wall.\n    selem = disk(10) # CHANGE BACK TO 10\n    binary = binary_closing(binary, selem)\n    \n    if plot:\n        plt.imshow(binary)\n        plt.show() \n        \n    # Step 7: Fill in the small holes inside the binary mask of lungs.\n    edges = roberts(binary)\n    \n    if plot:\n        plt.imshow(edges)\n        plt.show() \n        \n    binary = ndi.binary_fill_holes(edges)\n    \n    if plot:\n        plt.imshow(binary)\n        plt.show() \n        \n    # Step 8: Superimpose the binary mask on the input image.\n    selem = disk(4)\n    binary = dilation(binary, selem)\n    get_high_vals = binary == 0\n    im[get_high_vals] = -2000\n    \n    if plot:\n        plt.imshow(im)\n        plt.show()\n        \n    return im, binary","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Read the train file"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/rsna-str-pulmonary-embolism-detection/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Choose one CT image to check the segment result"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path = \"../input/rsna-str-pulmonary-embolism-detection/train/\"\n\nstudyID, SeriesID, SOPID = train.loc[50,['StudyInstanceUID','SeriesInstanceUID','SOPInstanceUID']].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dicom = data_path+studyID+\"/\"+SeriesID+\"/\"+SOPID+\".dcm\"\nimg = dcm.dcmread(dicom)\nimg_data = img.pixel_array # Read the pixel value\nhu = apply_modality_lut(img_data, img) # Transform to HU value\nlung_seg, _ = get_segmented_lungs(hu)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (10,10))\nplt.subplot(121)\nplt.imshow(hu)\nplt.subplot(122)\nplt.imshow(lung_seg)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Load one series and segment lung "},{"metadata":{"trusted":true},"cell_type":"code","source":"one_series_path = data_path+studyID+\"/\"+SeriesID+\"/\"\none_series = []\none_series_seg = []\n\nfor i in os.listdir(one_series_path):\n    dicom_path = one_series_path+\"/\"+i\n    img = dcm.dcmread(dicom_path)\n    img_data = img.pixel_array\n    hu = apply_modality_lut(img_data, img)\n    img_seg, _ = get_segmented_lungs(hu)\n    length = int(img.InstanceNumber)\n    one_series.append((length, img_data))\n    one_series_seg.append((length, img_seg))\n\none_series.sort()\none_series_seg.sort()\none_series_seg = [s[1] for s in one_series_seg]\none_series = [s[1] for s in one_series]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Animation\n* reference by @isaienkov greate [notebooks](https://www.kaggle.com/isaienkov/pulmonary-embolism-detection-eda/comments)"},{"metadata":{"trusted":true},"cell_type":"code","source":"from matplotlib import animation, rc\nrc('animation', html='jshtml')\n\ndef animate(ims,ims_seg):\n    fig , (ax1, ax2) = plt.subplots(1,2,figsize=(15,8))\n    ax1.axis('off')\n    ax2.axis('off')\n    im = ax1.imshow(ims[0])\n    im2 = ax2.imshow(ims_seg[0])\n\n    def animate_func(i):\n        im.set_data(ims[i])\n        im2.set_data(ims_seg[i])\n        return [im,im2]\n\n    anim = animation.FuncAnimation(fig, animate_func, frames = len(ims), interval = 1000//24)\n    \n    return anim\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"movie = animate(one_series,one_series_seg)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"movie","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Not sure lung segmentation is useful in this competition or not, just sharing how I processing lung CT in my work\n\n### Thanks for watch !!"}],"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":4,"nbformat_minor":4}