{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport scipy.ndimage as ndimage\nfrom skimage import measure, morphology, segmentation\nimport matplotlib.pyplot as plt\nimport os\nfrom pathlib import Path\nimport cv2\n\nimport time","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA_DIR = Path('/kaggle/input/osic-pulmonary-fibrosis-progression/')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Reference Kernels\n\n- https://www.kaggle.com/allunia/pulmonary-fibrosis-dicom-preprocessing\n- https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial\n- https://www.kaggle.com/aadhavvignesh/lung-segmentation-by-marker-controlled-watershed","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install ../input/lungmask/lungmask/dependencies/SimpleITK-1.2.4-cp37-cp37m-manylinux1_x86_64.whl\n!pip install ../input/lungmask/lungmask/dependencies/fastremap-1.10.2-cp37-cp37m-manylinux1_x86_64.whl\n!pip install ../input/lungmask/lungmask/dependencies/fill_voids-2.0.0-cp37-cp37m-manylinux1_x86_64.whl","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install ../input/lungmask/lungmask","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Dicom image utils","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Load the scans in given folder path\ndef load_scan(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n    for s in slices:\n        s.SliceThickness = slice_thickness\n    return slices\n\n#to HU values\ndef get_pixels_hu(slices):\n    image = np.stack([s.pixel_array for s in slices])\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slice_number] = image[slice_number].astype(np.int16)\n        image[slice_number] += np.int16(intercept)\n    return np.array(image, dtype=np.int16)\n\n# the scans can have different pixel size to real world mapping\n# resample to fix it to 1\ndef resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = np.array([float(scan[0].SliceThickness)] + list(scan[0].PixelSpacing), dtype=np.float32)\n\n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n    \n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Threshold based segmentation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def largest_label_volume(im, bg=-1):\n    vals, counts = np.unique(im, return_counts=True)\n    counts = counts[vals != bg]\n    vals = vals[vals != bg]\n    if len(counts) > 0:\n        return vals[np.argmax(counts)]\n    else:\n        return None\n\ndef segment_lung_mask(image, fill_lung_structures=True):\n    \n    # not actually binary, but 1 and 2. \n    # 0 is treated as background, which we do not want\n    binary_image = np.array(image > -320, dtype=np.int8)+1\n    labels = measure.label(binary_image)\n    \n    # Pick the pixel in the very corner to determine which label is air.\n    #   Improvement: Pick multiple background labels from around the patient\n    #   More resistant to \"trays\" on which the patient lays cutting the air \n    #   around the person in half\n    background_label = labels[0,0,0]\n    \n    #Fill the air around the person\n    binary_image[background_label == labels] = 2\n    \n    \n    # Method of filling the lung structures (that is superior to something like \n    # morphological closing)\n    if fill_lung_structures:\n        # For every slice we determine the largest solid structure\n        for i, axial_slice in enumerate(binary_image):\n            axial_slice = axial_slice - 1\n            labeling = measure.label(axial_slice)\n            l_max = largest_label_volume(labeling, bg=0)\n            \n            if l_max is not None: #This slice contains some lung\n                binary_image[i][labeling != l_max] = 1\n\n    \n    binary_image -= 1 #Make the image actual binary\n    binary_image = 1-binary_image # Invert it, lungs are now 1\n    \n    # Remove other air pockets insided body\n    labels = measure.label(binary_image, background=0)\n    l_max = largest_label_volume(labels, bg=0)\n    if l_max is not None: # There are air pockets\n        binary_image[labels != l_max] = 0\n \n    return binary_image\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Marker conrolled watershed segmentation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef generate_markers(image):\n    \"\"\"\n    Generates markers for a given image.\n    \n    Parameters: image\n    \n    Returns: Internal Marker, External Marker, Watershed Marker\n    \"\"\"\n    \n    #Creation of the internal Marker\n    marker_internal = image < -400\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    \n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    \n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    \n    marker_internal = marker_internal_labels > 0\n    \n    # Creation of the External Marker\n    external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n    \n    # Creation of the Watershed Marker\n    marker_watershed = np.zeros(image.shape, dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seperate_lungs(image, iterations = 1):\n    \"\"\"\n    Segments lungs using various techniques.\n    \n    Parameters: image (Scan image), iterations (more iterations, more accurate mask)\n    \n    Returns: \n        - Segmented Lung\n        - Lung Filter\n        - Outline Lung\n        - Watershed Lung\n        - Sobel Gradient\n    \"\"\"\n    \n    # Store the start time\n    start = time.time()\n    \n    marker_internal, marker_external, marker_watershed = generate_markers(image)\n    \n    \n    '''\n    Creation of Sobel Gradient\n    '''\n    \n    # Sobel-Gradient\n    sobel_filtered_dx = ndimage.sobel(image, 1)\n    sobel_filtered_dy = ndimage.sobel(image, 0)\n    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n    sobel_gradient *= 255.0 / np.max(sobel_gradient)\n    \n    \n    '''\n    Using the watershed algorithm\n    \n    \n    We pass the image convoluted by sobel operator and the watershed marker\n    to morphology.watershed and get a matrix matrix labeled using the \n    watershed segmentation algorithm.\n    '''\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n    \n    '''\n    Reducing the image to outlines after Watershed algorithm\n    '''\n    outline = ndimage.morphological_gradient(watershed, size=(3,3))\n    outline = outline.astype(bool)\n    \n    \n    '''\n    Black Top-hat Morphology:\n    \n    The black top hat of an image is defined as its morphological closing\n    minus the original image. This operation returns the dark spots of the\n    image that are smaller than the structuring element. Note that dark \n    spots in the original image are bright spots after the black top hat.\n    '''\n    \n    # Structuring element used for the filter\n    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [0, 0, 1, 1, 1, 0, 0]]\n    \n    blackhat_struct = ndimage.iterate_structure(blackhat_struct, iterations)\n    \n    # Perform Black Top-hat filter\n    outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n    \n    '''\n    Generate lung filter using internal marker and outline.\n    '''\n    lungfilter = np.bitwise_or(marker_internal, outline)\n    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n    \n    '''\n    Segment lung using lungfilter and the image.\n    '''\n    segmented = np.where(lungfilter == 1, image, -2000*np.ones(image.shape))\n    \n    return segmented, lungfilter, outline, watershed, sobel_gradient\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Segmentation using Lungmask","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from lungmask import mask\nimport SimpleITK as sitk","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mask_model = mask.load_model('unet','R231', '../input/lungmask/lungmask/models/unet_r231-d5d2fc3d.pth')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_img_sitk(path):\n    return sitk.ReadImage(path)\n\ndef generate_mask_image(path):\n    img = get_img_sitk(path)\n    segmentation = (mask.apply(img, mask_model)[0,:,:])\n    segmentation[segmentation > 0] = 1\n    segmentation = cv2.resize(segmentation, (512, 512))\n    img_array = cv2.resize(((sitk.GetArrayFromImage(img)[0,:,:]) - float(img.GetMetaData('0028|1052'))) / (float(img.GetMetaData('0028|1053')) * 1000), (512, 512))\n    masked_img = np.where(segmentation == 1, img_array, 0)\n    return segmentation, masked_img","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Image Normalization","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"MIN_BOUND = -1000.0\nMAX_BOUND = 400.0\n    \ndef normalize(image):\n    image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)\n    image[image>1] = 1.\n    image[image<0] = 0.\n    return image\n\n##Zero centering\nPIXEL_MEAN = 0.25\n\ndef zero_center(image):\n    image = image - PIXEL_MEAN\n    return image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visualizations","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_id = 'ID00007637202177411956430'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(DATA_DIR/'train.csv')\ntrain[train['Patient'] == patient_id]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slices = load_scan(str(DATA_DIR/f'train/{patient_id}'))\nhu_slices = get_pixels_hu(slices)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_hist(arr):\n    plt.hist(arr.flatten(), bins=80, color='c')\n    plt.xlabel(\"Hounsfield Units (HU)\")\n    plt.ylabel(\"Frequency\")\n    plt.show()\n\n# Show some slice in the middle\ndef plot_slice(arr, n=80):\n    plt.imshow(arr[n], cmap=plt.cm.gray)\n    plt.show()\n    \n#watershed plot markers\ndef plot_watershed_markers(hu_slice):\n    test_patient_internal, test_patient_external, test_patient_watershed = generate_markers(hu_slice)\n\n    f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(15,15))\n\n    ax1.imshow(test_patient_internal, cmap='gray')\n    ax1.set_title(\"Internal Marker\")\n    ax1.axis('off')\n\n    ax2.imshow(test_patient_external, cmap='gray')\n    ax2.set_title(\"External Marker\")\n    ax2.axis('off')\n\n    ax3.imshow(test_patient_watershed, cmap='gray')\n    ax3.set_title(\"Watershed Marker\")\n    ax3.axis('off')\n\n    plt.show()\n    \ndef plot_watershed_results(hu_slice, itrs=1):\n    test_segmented, test_lungfilter, test_outline, test_watershed, test_sobel_gradient = seperate_lungs(hu_slice, itrs)\n    f, ax = plt.subplots(3, 2, sharey=True, figsize = (12, 12))\n    ax[0][0].imshow(test_sobel_gradient, cmap='gray')\n    ax[0][0].set_title(\"Sobel Gradient\")\n    ax[0][0].axis('off')\n\n    ax[0][1].imshow(test_watershed, cmap='gray')\n    ax[0][1].set_title(\"Watershed\")\n    ax[0][1].axis('off')\n    \n    ax[1][0].imshow(test_segmented, cmap='gray')\n    ax[1][0].set_title('Segmented Lung')\n    ax[1][0].axis('off')\n    \n    ax[1][1].imshow(test_lungfilter, cmap='gray')\n    ax[1][1].set_title('Lungfilter')\n    ax[1][1].axis('off')\n    \n    ax[2][0].imshow(test_outline, cmap='gray')\n    ax[2][0].set_title('Outline')\n    ax[2][0].axis('off')\n\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_hist(hu_slices)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_slice(hu_slices, 20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# the images rearraged by instance numbers\nplt.imshow(pydicom.read_file(str(DATA_DIR/f'train/{patient_id}/20.dcm')).pixel_array, cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"st = time.time()\nresampled_slices, new_spacing = resample(hu_slices, slices)\nprint(new_spacing)\nprint(resampled_slices.shape)\nprint(time.time() - st)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_slice(resampled_slices,200)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Watershed segmentation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"st = time.time()\ntest_segmented, test_lungfilter, test_outline, test_watershed, test_sobel_gradient = seperate_lungs(resampled_slices[200], 1)\nprint(time.time() - st)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_watershed_results(resampled_slices[200])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_hist(test_segmented[test_segmented > -1000])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Lungmask library segmentation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"st = time.time()\nmask_filter, masked_lung = generate_mask_image(str(DATA_DIR/f'train/{patient_id}/20.dcm'))\nprint(time.time() - st)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(mask_filter,cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(masked_lung, cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Threshold based segmentation","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# this is fastest, but doesn't seem reliable, or maybe I'm missing something\n# There is no proper segmentation in the source kernel too\n# https://www.kaggle.com/allunia/pulmonary-fibrosis-dicom-preprocessing\nst = time.time()\nsegmented_lungs = segment_lung_mask(resampled_slices, False)\nprint((time.time() - st) / resampled_slices.shape[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(segmented_lungs[200, :, :], cmap='gray')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}