{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 3D Interactive DICOM Viz : Pulmonary Embolism\n\n\n**Digital Imaging and Communications in Medicine (DICOM)** is the standard for the communication and management of medical imaging information and related data. I had a hard time understanding DICOM. \n\n\nThe 3D interactive images by @aravrs  helped me a lot!  Inspired by his work, I am trying the same for this competition. \n\nHope this helps you visualise DICOM as well!\n\n___"},{"metadata":{},"cell_type":"markdown","source":"# Import libraries"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import os\nimport plotly\nimport pydicom\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom glob import glob\nimport scipy.ndimage\nfrom skimage import measure\nimport matplotlib.pyplot as plt\nimport plotly.graph_objects as go\nfrom plotly.figure_factory import create_trisurf\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":" # Define functions"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def load_scan(path, reverse=True):\n    slices = [pydicom.read_file(path + \"/\" + s) for s in os.listdir(path)]\n    slices.sort(key=lambda x: int(x.InstanceNumber), reverse=reverse)\n\n    try:\n        slice_thickness = np.abs(\n            slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2]\n        )\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n\n    for s in slices:\n        s.SliceThickness = slice_thickness\n\n    return slices\n\n\ndef get_pixels_hu(scans):\n    image = np.stack([s.pixel_array for s in scans])\n    image = image.astype(np.int16)\n    image[image == -2000] = 0\n\n    intercept = scans[0].RescaleIntercept\n    slope = scans[0].RescaleSlope\n\n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n\n    image += np.int16(intercept)\n    return np.array(image, dtype=np.int16)\n\n\ndef resample(image, scan, new_spacing=[1, 1, 1]):\n    spacing = map(float, ([scan[0].SliceThickness] + list(scan[0].PixelSpacing)))\n    spacing = np.array(list(spacing))\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)\n    return image, new_spacing\n\n\ndef make_mesh(image, threshold=-300, step_size=1):\n    p = image.transpose(2, 1, 0)\n    verts, faces, _, _ = measure.marching_cubes_lewiner(p, threshold, step_size=step_size, allow_degenerate=True)\n    return verts, faces\n\n\ndef 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\n\ndef segment_lung_mask(image, fill_lung_structures=True):\n    binary_image = np.array(image >= -700, dtype=np.int8) + 1\n    labels = measure.label(binary_image)\n    background_label = labels[0, 0, 0]\n    binary_image[background_label == labels] = 2\n\n    if fill_lung_structures:\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:\n                binary_image[i][labeling != l_max] = 1\n    binary_image -= 1\n    binary_image = 1 - binary_image\n\n    labels = measure.label(binary_image, background=0)\n    l_max = largest_label_volume(labels, bg=0)\n    if l_max is not None:\n        binary_image[labels != l_max] = 0\n\n    return binary_image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Analysis"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"TRAIN_DIR = \"../input/rsna-str-pulmonary-embolism-detection/train/\"\nsub_folder_list = []\nfor x in os.listdir(TRAIN_DIR):\n    if os.path.isdir(TRAIN_DIR + '/' + x):\n        for y in os.listdir(TRAIN_DIR+ '/' + x):\n            if os.path.isdir(TRAIN_DIR + '/' +x  + '/'+ y):\n                sub_folder_list.append(x  + '/'+ y)\n\nn_dicom_dict = {\"Patient\":[],\"n_dicom\":[]}\nfor x in sub_folder_list:\n    g = glob(TRAIN_DIR+x + '/*.dcm')\n    n_dicom_dict[\"n_dicom\"].append(len(g))\n    n_dicom_dict[\"Patient\"].append(x)\n\ndicom_df = pd.DataFrame(n_dicom_dict)\ndicom_df.sort_values(['n_dicom'], inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(\"Minimum number of dicom files:\", min(dicom_df['n_dicom']))\nprint(\"Maximum number of dicom files:\", max(dicom_df['n_dicom']))\n\nplt.figure(figsize=(20,10))\nsns.distplot(dicom_df['n_dicom'], bins=20, color=\"#a55eea\")\nplt.title('Number of dicom files per patient');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3D Visualizations"},{"metadata":{},"cell_type":"markdown","source":"Let's consider patient `004b06506f54/a6862824bdd3/` from `/train` as an example, and visualize."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"data_path = \"../input/rsna-str-pulmonary-embolism-detection/train/004b06506f54/a6862824bdd3/\"\nprint(f\"{data_path.split('/')[-3].upper()} - {data_path.split('/')[-2]}\")\ng = glob(data_path + \"/*.dcm\")\nprint(f\"Total of {len(g)} DICOM images.\")\n\npatient = load_scan(data_path, False)\nprint(f\"Slice Thickness: {patient[0].SliceThickness}\")\nprint(f\"Pixel Spacing (row, col): ({patient[0].PixelSpacing[0]}, {patient[0].PixelSpacing[1]})\")\n\nimgs = get_pixels_hu(patient)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def plot_3d(data_path, reverse=False):\n    print(f\"{data_path.split('/')[-3].upper()} - {data_path.split('/')[-2]}\")\n    g = glob(data_path + \"/*.dcm\")\n    print(f\"Total of {len(g)} DICOM images.\")\n\n    patient = load_scan(data_path, reverse)\n    print(f\"Slice Thickness: {patient[0].SliceThickness}\")\n    print(f\"Pixel Spacing (row, col): ({patient[0].PixelSpacing[0]}, {patient[0].PixelSpacing[1]})\")\n\n    imgs = get_pixels_hu(patient)\n    print(f\"Shape resampling: {imgs.shape}\", end=\"\")\n    imgs_after_resamp, spacing = resample(imgs, patient, [1, 1, 1])\n    print(f\" -> {imgs_after_resamp.shape}\")\n\n    v1, f1 = make_mesh(imgs_after_resamp, 350, 2)\n\n    segmented_lungs = segment_lung_mask(imgs_after_resamp, fill_lung_structures=False)\n    segmented_lungs_fill = segment_lung_mask(imgs_after_resamp, fill_lung_structures=True)\n    internal_structures = segmented_lungs_fill - segmented_lungs\n    p = internal_structures.transpose(2, 1, 0)\n    v2, f2, _, _ = measure.marching_cubes_lewiner(p)\n\n    ### PLOTS\n    fig = plt.figure(figsize=(20, 10))\n    bg = np.array((30, 39, 46))/255.0\n    \n    # Ext\n    print(\".\", end=\"\")\n    x, y, z = zip(*v1)\n    ax1 = fig.add_subplot(121, projection=\"3d\")\n    mesh = Poly3DCollection(v1[f1], alpha=0.8)\n    face_color = (1, 1, 0.9)\n    mesh.set_facecolor(face_color)\n    ax1.add_collection3d(mesh)\n    ax1.set_xlim(0, max(x))\n    ax1.set_ylim(0, max(y))\n    ax1.set_zlim(0, max(z))\n    ax1.w_xaxis.set_pane_color((*bg, 1))\n    ax1.w_yaxis.set_pane_color((*bg, 1))\n    ax1.w_zaxis.set_pane_color((*bg, 1))\n\n    # Int\n    print(\".\", end=\"\")\n    x, y, z = zip(*v2)\n    ax2 = fig.add_subplot(122, projection=\"3d\")\n    mesh = Poly3DCollection(v2[f2], alpha=0.8)\n    face_color = np.array((255, 107, 107))/255.0\n    mesh.set_facecolor(face_color)\n    ax2.add_collection3d(mesh)\n    ax2.set_xlim(0, max(x))\n    ax2.set_ylim(0, max(y))\n    ax2.set_zlim(0, max(z))\n    ax2.w_xaxis.set_pane_color((*bg, 1))\n    ax2.w_yaxis.set_pane_color((*bg, 1))\n    ax2.w_zaxis.set_pane_color((*bg, 1))\n\n    print(\".\", end=\"\")\n    fig.tight_layout()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot_3d(data_path)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3D Interactive Plots"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"### 3D interactive ploting helper\ndef plotly_3d(verts, faces, ext=True):\n    x, y, z = zip(*verts)\n\n    fig = create_trisurf(\n        x=x,\n        y=y,\n        z=z,\n        plot_edges=False,\n        show_colorbar=False,\n        showbackground=False,\n        colormap=[\"rgb(236, 236, 212)\", \"rgb(236, 236, 212)\"] if ext else [\"rgb(255, 107, 107)\", \"rgb(255, 107, 107)\"],\n        simplices=faces,\n        backgroundcolor=\"rgb(30, 39, 46)\",\n        gridcolor=\"rgb(30, 39, 46)\",\n        title=\"<b>Interactive Visualization</b>\",\n    )\n    fig.layout.template = \"plotly_dark\"  # for dark theme \n    fig.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"### Ploting functions\n\n\ndef plot3d_interactive_ext(data_path, reverse=False):\n    print(f\"{data_path.split('/')[-3].upper()} - {data_path.split('/')[-2]}\")\n    g = glob(data_path + \"/*.dcm\")\n    patient = load_scan(data_path, reverse)\n    imgs = get_pixels_hu(patient)\n    imgs_after_resamp, spacing = resample(imgs, patient, [1, 1, 1])\n\n    v, f = make_mesh(imgs_after_resamp, 350, 2)\n    plotly_3d(v, f)\n\n\ndef plot3d_interactive_int(data_path, reverse=False):\n    print(f\"{data_path.split('/')[-3].upper()} - {data_path.split('/')[-2]}\")\n    g = glob(data_path + \"/*.dcm\")\n    patient = load_scan(data_path, reverse)\n    imgs = get_pixels_hu(patient)\n    imgs_after_resamp, spacing = resample(imgs, patient, [1, 1, 1])\n\n    segmented_lungs = segment_lung_mask(imgs_after_resamp, fill_lung_structures=False)\n    segmented_lungs_fill = segment_lung_mask(imgs_after_resamp, fill_lung_structures=True)\n    internal_structures = segmented_lungs_fill - segmented_lungs\n\n    p = internal_structures.transpose(2, 1, 0)\n    verts, faces, _, _ = measure.marching_cubes_lewiner(p)\n    plotly_3d(verts, faces, ext=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The below is an interactive plot. One can zoom in and out as well as rotate it!!"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot3d_interactive_ext(data_path)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot3d_interactive_int(data_path)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Overlapping scans"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot3d_interactive_ext(\"../input/rsna-str-pulmonary-embolism-detection/train/004b06506f54/a6862824bdd3/\", True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Acknowlegements and References\n\n- https://www.kaggle.com/aravrs/3d-dicom-visualizations-with-interactive-plots (main)\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\n- https://www.kaggle.com/wrrosa/advanced-dicom-ct-3d-visualizations-with-vtk\n"}],"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}