{"cells":[{"metadata":{},"cell_type":"markdown","source":"# ABOUT\n\nthis kernel is to preprocess DICOM images and 3D plot with ploty from the competition [OSIC Pulmonary Fibrosis Progression](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression).\n\nDICON images are cropped, resized, segmentated, converted into Poly3DCollection and then saved into npy files for later use in deep learning.\n\nI will keep working on it and improve preprocessings\n\nNote that it is highly based on [OSIC |quick EDA + 3D plot with plotly](https://www.kaggle.com/sunpnwt12/osic-quick-eda-3d-plot-with-plotly)"},{"metadata":{},"cell_type":"markdown","source":"# Import"},{"metadata":{"trusted":true,"scrolled":true,"_kg_hide-input":true},"cell_type":"code","source":"# !conda install -c conda-forge gdcm -y # run this code for the first time\nimport os\n# import gdcm\nfrom collections import defaultdict\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# import plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly import figure_factory as FF\n\nimport scipy.ndimage\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\n\nimport random\nimport pydicom","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Read files"},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"# load train and test image\n\nDICOM_DIR = '/kaggle/input/osic-pulmonary-fibrosis-progression/train'\nDICOM_DIR_TEST = '/kaggle/input/osic-pulmonary-fibrosis-progression/test'\n\ndicom_dict = defaultdict(list)\ndicom_dict_test = defaultdict(list)\n\ndefault_image_size = 512\n\nfor dirname in os.listdir(DICOM_DIR):\n    path = os.path.join(DICOM_DIR, dirname)\n    dicom_dict[dirname].append(path)\n    \nfor dirname in os.listdir(DICOM_DIR_TEST):\n    path = os.path.join(DICOM_DIR_TEST, dirname)\n    dicom_dict_test[dirname].append(path)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"# ### load_scan:\n# since there are a couple of dicom files that don't have 'ImagePositionPatient' attribute so instead,\n# i will use 'InstanceNumber' attribute for those\n\n# ### dicom_file:\n# 1. take index number of patient which stored in dict_dicom earlier\n# 2. this might be useful when you need to pick some random patient\n# 3. It also takes specific patient Id in case you need.\n# 4. Note that this function is going to read all file in taken path.\n\n# ### get_pixels_hu\n# 1. take dicom file which had called through dicom_file function\n# 2. It stacks up all the load slices of certain patient\n# 3. stacked slices will be calculated into Hounsfield Units\n\n\n\ndef load_scan(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    length = len(slices)\n    actualLength = sum([1 if hasattr(x, 'ImagePositionPatient') else 0 for x in slices])\n    if length == actualLength:\n        slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    else:\n        slices.sort(key = lambda x: float(x.InstanceNumber))\n        \n    return slices\n\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    \n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        \n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slbice_number] = image[slice_number].astype(np.int16)\n            \n        image[slice_number] += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# 3D Visualization"},{"metadata":{},"cell_type":"markdown","source":"Useful links: \n\n* https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\n* https://medium.com/@hengloose/a-comprehensive-starter-guide-to-visualizing-and-analyzing-dicom-images-in-python-7a8430fcb7ed"},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"test = load_scan(dicom_dict['ID00047637202184938901501'][0])\n\ntest_hu = get_pixels_hu(test)\nprint('Patient {}'.format(test[0].PatientName))\nprint('Slices : {}\\nPixels : ({} x {})'.format(test_hu.shape[0], test_hu.shape[1], test_hu.shape[2]))\n\nplt.figure(figsize=(12, 8))\nax = sns.distplot(test_hu.flatten(), bins=80, norm_hist=True)\n# ax.set_title('Hounsfield Units of patient {}'.format(test_hu[0].PatientName), fontsize=25)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ref) https://stackoverflow.com/questions/48121916/numpy-resize-rescale-image\n# INTER_NEAREST - a nearest-neighbor interpolation\n# INTER_LINEAR - a bilinear interpolation (used by default)\n# INTER_AREA - resampling using pixel area relation. It may be a preferred method for image decimation, as it gives moire’-free results. But when the image is zoomed, it is similar to the INTER_NEAREST method.\n# INTER_CUBIC - a bicubic interpolation over 4x4 pixel neighborhood\n# INTER_LANCZOS4 - a Lanczos interpolation over 8x8 pixel neighborhood\n\nimport cv2\ndef resize(slices):\n    new_slice = []\n    for slice_number in range(len(slices)):\n        new_slice.append(cv2.resize(slices[slice_number], dsize=(default_image_size, default_image_size), interpolation=cv2.INTER_AREA))\n    return np.array(new_slice)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # do crop here\ndef crop(slices):\n    slice_height = len(slices[0])\n    slice_width = len(slices[0][0])\n    diff_height_half = int((slice_height - default_image_size) / 2)\n    diff_width_half = int((slice_width - default_image_size) / 2)\n    new_slice = []\n    for slice_number in range(len(slices)):\n        new_slice.append(slices[slice_number][diff_width_half: slice_width - diff_width_half, diff_height_half: slice_height - diff_height_half])\n    return np.array(new_slice)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"scrolled":true},"cell_type":"code","source":"def resample(image, slice_tickness, x_pixel_spacing, y_pixel_spacing, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = np.array([slice_tickness, x_pixel_spacing, y_pixel_spacing], 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 = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n    \n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def make_mesh(image, threshold):\n    p = image.transpose(2, 1, 0)\n    \n    verts, faces, normals, values = measure.marching_cubes_lewiner(p, threshold)\n    return verts, faces\n\ndef static_3d(image, threshold=-300):\n    \n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    \n    verts, faces = make_mesh(image, threshold)\n    x, y, z = zip(*verts)\n    \n    mesh = Poly3DCollection(verts[faces], alpha=0.1)\n    face_color = [0.5, 0.5, 1]\n    mesh.set_facecolor(face_color)\n    \n    ax.add_collection3d(mesh)\n    ax.set_xlim(0, max(x))\n    ax.set_ylim(0, max(y))\n    ax.set_zlim(0, max(z))\n    plt.show()\n    \ndef interactive_3d(image, threshold=-300):\n    verts, faces = make_mesh(image, threshold)\n    x, y, z = zip(*verts)\n    fig = FF.create_trisurf(x=x,\n                            y=y,\n                            z=z,\n                            plot_edges=False,\n                            simplices=faces)\n    iplot(fig)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"x_size = test[0].Columns\ny_size = test[0].Rows\nslice_thinkness = test[0].SliceThickness\nx_pixel_spacing = test[0].PixelSpacing[0] * (x_size / default_image_size)\ny_pixel_spacing = test[0].PixelSpacing[1] * (y_size / default_image_size)\n\n# if ratio is not 1:1 then crop\nif (x_size != y_size):\n    test_hu = crop(test_hu)\n\n# if size is not 512 then resize\nif (x_size != default_image_size or y_size != default_image_size):\n    test_hu = resize(test_hu)\n\n# resample test\nresampled_test_hu, spacing = resample(test_hu, slice_thinkness, x_pixel_spacing, y_pixel_spacing)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The ratio of images that are not 1:1 are the ones with margins outside of images. crop those margins and resize to 512 resolution\n\nAs total resolution gets smaller, I muliplied pixel spacings by the ration too ( I guess as total number of pixel gets smaller, single pixel covers larger space thus bigger millimetres ) "},{"metadata":{"trusted":true},"cell_type":"code","source":"static_3d(resampled_test_hu)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"scrolled":true},"cell_type":"code","source":"def largest_label_volume(im, bg=-1):\n    vals, counts = np.unique(im, return_counts=True)\n\n    counts = counts[vals != bg]\n    vals = vals[vals != bg]\n\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    # 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    z, y, x = labels.shape\n    for slice_number in range(len(image)):\n        #Fill the air around the person\n        binary_image[slice_number][labels[slice_number, 0, 0] == labels[slice_number]] = 2\n        binary_image[slice_number][labels[slice_number, 0, x - 1] == labels[slice_number]] = 2\n        binary_image[slice_number][labels[slice_number, y - 1, 0] == labels[slice_number]] = 2\n        binary_image[slice_number][labels[slice_number, y - 1, x - 1] == labels[slice_number]] = 2\n        binary_image[slice_number][labels[slice_number, 0, int((x - 1) / 2)] == labels[slice_number]] = 2\n        binary_image[slice_number][labels[slice_number, y - 1, int((x - 1) / 2)] == labels[slice_number]] = 2\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    return binary_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"segmented_lungs = segment_lung_mask(resampled_test_hu, False)\nsegmented_lungs_fill = segment_lung_mask(resampled_test_hu, True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Segmented lungs"},{"metadata":{"trusted":true},"cell_type":"code","source":"static_3d(segmented_lungs, 1.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Segmented lungs filled"},{"metadata":{"trusted":true},"cell_type":"code","source":"static_3d(segmented_lungs_fill, 1.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Difference between both"},{"metadata":{"trusted":true},"cell_type":"code","source":"static_3d(segmented_lungs_fill - segmented_lungs, -0.5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"interactive_3d(segmented_lungs_fill - segmented_lungs, -0.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":" # Save to npy\n"},{"metadata":{},"cell_type":"markdown","source":"the codes below is to save segmented lung data into npy files."},{"metadata":{"trusted":true},"cell_type":"code","source":"# create directories\n# os.removedirs(\"/kaggle/working/train\")\n# os.removedirs(\"/kaggle/working/test\")\n# os.makedirs('train')\n# os.makedirs('test')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load all dicom\n# import gdcm\n\n# PATH = '/kaggle/working/train'\n# PATH_TEST = '/kaggle/working/test'\n\n# # create train data npy\n# for k, v in dicom_dict.items():\n#     print(k ,v)\n#     scan = load_scan(v[0])\n#     scan_hu = get_pixels_hu(scan)\n#     resample_scan_hu, spacing = resample(scan_hu, scan)\n#     segmented_lungs = segment_lung_mask(resample_scan_hu, False)\n#     segmented_lungs_fill = segment_lung_mask(resample_scan_hu, True)\n#     np.save(f'{PATH}/{k}' , segmented_lungs_fill - segmented_lungs)\n    \n# create test data npy\n# for k, v in dicom_dict_test.items():\n#     print(k ,v)\n#     scan = load_scan(v[0])\n#     scan_hu = get_pixels_hu(scan)\n#     resample_scan_hu, spacing = resample(scan_hu, scan)\n#     segmented_lungs = segment_lung_mask(resample_scan_hu, False)\n#     segmented_lungs_fill = segment_lung_mask(resample_scan_hu, True)\n#     np.save(f'{PATH_TEST}/{k}' , segmented_lungs_fill - segmented_lungs)","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}