{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"# heavily based on\nhttps://www.kaggle.com/sunpnwt12/osic-quick-eda-3d-plot-with-plotly","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import os\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\nimport 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\n\nINPUT_DIR = '/kaggle/input/osic-pulmonary-fibrosis-progression'\n\ntrainset = pd.read_csv(f'{INPUT_DIR}/train.csv')\ntestset = pd.read_csv(f'{INPUT_DIR}/test.csv')\nsample_sub = pd.read_csv(f'{INPUT_DIR}/sample_submission.csv')\n\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        \n    return slices\n\ndef dicom_file(idx_num, patient_id=None):\n    if patient_id:\n        return load_scan(dicom_dict[patient_id][0])\n    return load_scan(dicom_dict[p_id[idx_num]][0])\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\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)\n\n\ndef resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = np.array([scan[0].SliceThickness, scan[0].PixelSpacing[0], scan[0].PixelSpacing[1]], 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\n\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        \n    return slices\n\ndef dicom_file(idx_num, patient_id=None):\n    if patient_id:\n        return load_scan(dicom_dict[patient_id][0])\n    return load_scan(dicom_dict[p_id[idx_num]][0])\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\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)\n\ndef 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, angle=0):\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    return ax,\n\nDICOM_DIR = '/kaggle/input/osic-pulmonary-fibrosis-progression/train'\n\ndicom_dict = defaultdict(list)\n\n\nfor dirname in os.listdir(DICOM_DIR):\n    path = os.path.join(DICOM_DIR, dirname)\n    dicom_dict[dirname].append(path)\n    \np_id = sorted(trainset['Patient'].unique())\n\ntest2 = dicom_file(40)\n\ntest2_hu = get_pixels_hu(test2)\n\nresampled_test2_hu, spacing = resample(test2_hu, test2)\n\nfrom scipy.ndimage import rotate\n\nsub = resampled_test2_hu[::5,::5,::5]\n\nfrom matplotlib import animation\n\nTHRESHOLD = -300\n\nfig,_ = plt.subplots(figsize=(11,11))\nax = fig.add_subplot(111, projection='3d')\n\ndef make_mesh(image, threshold):\n    p = image.transpose(2, 1, 0)\n    verts, faces, normals, values = measure.marching_cubes_lewiner(p, threshold)\n    return verts, faces\n\ndef static_3d(image,ax, angle=0):\n    global x\n    global y\n    global z\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    return ax\n\n\n\ndef init_time_line(image=sub, axi=ax):\n    axi = static_3d(image, axi)\n    return axi\n\ndef animate_all(angle,axi=ax):\n    \n    image = (rotate(sub, angle=angle, axes=(1,2), reshape=False,cval=-1024) )\n    verts, faces = make_mesh(image, THRESHOLD)   \n    \n    axi.clear()\n    mesh = Poly3DCollection(verts[faces], alpha=0.1)\n    face_color = [0.5, 0.5, 1]\n    mesh.set_facecolor(face_color)\n    axi.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    return axi\n\nframes = np.arange(0,360, 6)\n\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"anim = animation.FuncAnimation(fig, animate_all, init_func=init_time_line, frames=frames, interval=41, blit=False)\n\nfrom matplotlib import animation, rc\nrc('animation', html='jshtml')\nanim","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}