{"cells":[{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pydicom\nimport os\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly.figure_factory import create_trisurf\nfrom plotly.graph_objs import *\n\nfrom glob import glob\n\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input_path = '../input/osic-pulmonary-fibrosis-progression/test/ID00419637202311204720264'\n\noutput_path = working_path = \"/kaggle/working/\"\ng = glob(input_path + '/*.dcm')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print (\"Total of %d DICOM images.\\nFirst 5 filenames:\" % len(g))\nprint ('\\n'.join(g[:5]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load scan","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\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        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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    # Convert to Hounsfield units (HU)\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    \n    return np.array(image, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"id=2\npatient = load_scan(input_path)\nimgs = get_pixels_hu(patient)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.save(output_path + \"fullimages_%d.npy\" % (id), imgs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used=output_path+\"fullimages_%d.npy\" % id\nimgs_to_process = np.load(file_used).astype(np.float64) \n\nplt.hist(imgs_to_process.flatten(), bins=50, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"id = 2\nimgs_to_process = np.load(output_path+'fullimages_{}.npy'.format(id))\n\ndef sample_stack(stack, rows=2, cols=2, start_with=10, show_every=3):\n    fig,ax = plt.subplots(rows,cols,figsize=[12,12])\n    for i in range(rows*cols):\n        ind = start_with + i*show_every\n        ax[int(i/rows),int(i % rows)].set_title('slice %d' % ind)\n        ax[int(i/rows),int(i % rows)].imshow(stack[ind],cmap='gray')\n        ax[int(i/rows),int(i % rows)].axis('off')\n    plt.show()\n\nsample_stack(imgs_to_process)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print (\"Slice Thickness: %f\" % patient[0].SliceThickness)\nprint (\"Pixel Spacing (row, col): (%f, %f) \" % (patient[0].PixelSpacing[0], patient[0].PixelSpacing[1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"id = 2\nimgs_to_process = np.load(output_path+'fullimages_{}.npy'.format(id))\nprint (\"Shape before resampling\\t\", imgs_to_process.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    #spacing = map(float, ([scan[0].SliceThickness] + scan[0].PixelSpacing))\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    \n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"imgs_after_resamp, spacing = resample(imgs_to_process, patient)\nprint (\"Shape after resampling\\t\", imgs_after_resamp.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_mesh(image, threshold=-300, step_size=1):\n\n    print (\"Transposing surface\")\n    p = image.transpose(2,1,0)\n        \n    print(\"Calculating surface\")\n    verts, faces, norm, val = measure.marching_cubes_lewiner(p, threshold, step_size=step_size, allow_degenerate=True)\n    \n    return verts, faces\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plotly_3d(verts, faces):\n    x,y,z = zip(*verts) \n    \n    # Make the colormap single color since the axes are positional not intensity. \n    #colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']\n    colormap=['rgb(236, 236, 212)','rgb(236, 236, 212)']\n    \n    fig = create_trisurf(x=x,\n                        y=y, \n                        z=z, \n                        plot_edges=False,\n                        colormap=colormap,\n                        simplices=faces,\n                        backgroundcolor='rgb(64, 64, 64)',\n                        title=\"Interactive Visualization\")\n    iplot(fig)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plt_3d(verts, faces):\n    x,y,z = zip(*verts) \n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    print (\"Drawing\")\n    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n    mesh = Poly3DCollection(verts[faces], linewidths=0.05, alpha=1)\n    face_color = [1, 1, 0.9]\n    mesh.set_facecolor(face_color)\n    ax.add_collection3d(mesh)\n\n    ax.set_xlim(0, max(x))\n    ax.set_ylim(0, max(y))\n    ax.set_zlim(0, max(z))\n    ax.set_facecolor((0.7, 0.7, 0.7))\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"v, f = make_mesh(imgs_after_resamp, 350)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"v, f = make_mesh(imgs_after_resamp, 350, 5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plotly_3d(v, f)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt_3d(v, f)","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}