{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# %reload_ext signature\n# %matplotlib inline\n\nimport numpy as np\nimport pydicom\nimport os\nimport matplotlib.pyplot as plt\nfrom glob import glob\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize\nfrom sklearn.cluster import KMeans\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nimport plotly.figure_factory as ff\nfrom plotly.graph_objs import *\ninit_notebook_mode(connected=True) ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"patient_id = '6897fa9de148'\npatient_folder = f'../input/rsna-str-pulmonary-embolism-detection/train/{patient_id}/'\ndata_paths = glob(patient_folder + '/*/*.dcm')\n\n# Print out the first 5 file names to verify we're in the right folder.\nprint (f'Total of {len(data_paths)} DICOM images.\\nFirst 5 filenames:' )\ndata_paths[:5]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Loop over the image files and store everything into a list.\n# \n\ndef load_scan(paths):\n    slices = [pydicom.read_file(path ) for path in paths]\n    slices.sort(key = lambda x: int(x.InstanceNumber), reverse = True)\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\n\ndef get_pixels_hu(scans):\n    image = np.stack([s.pixel_array for s in scans])\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 1\n    # The intercept is usually -1024, so air is approximately 0\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":"patient_id = '0003b3d648eb'\npatient = load_scan(data_paths)\nimgs = get_pixels_hu(patient)\n\n#This is a good time to save the new data set to disk so we don't have to reprocess the stack every time.\nnp.save(f'fullimages_{patient_id}.npy', imgs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used= f'./fullimages_{patient_id}.npy'\nimgs_to_process = np.load(file_used).astype(np.float64) \n\nplt.hist(imgs_to_process.flatten(), bins=100, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used= f'./fullimages_{patient_id}.npy'\nimgs_to_process = np.load(file_used).astype(np.float64) \n\ndef sample_stack(stack, rows=6, cols=6, start_with=10, show_every=3):\n    fig,ax = plt.subplots(rows,cols,figsize=[18,20])\n    for i in range(rows*cols):\n        ind = start_with + i*show_every\n        ax[int(i/rows),int(i % rows)].set_title(f'slice {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(f'Slice Thickness: {patient[0].SliceThickness}')\nprint(f'Pixel Spacing (row, col): ({patient[0].PixelSpacing[0]}, {patient[0].PixelSpacing[1]})')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.array([float(patient[0].SliceThickness), float(patient[0].PixelSpacing[0]), float(patient[0].PixelSpacing[0])])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_used= f'./fullimages_{patient_id}.npy'\nimgs_to_process = np.load(file_used).astype(np.float64)\n\ndef 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 = np.array([float(patient[0].SliceThickness), \n                        float(patient[0].PixelSpacing[0]), \n                        float(patient[0].PixelSpacing[0])])\n\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\n\nprint(f'Shape before resampling: {imgs_to_process.shape}')\nimgs_after_resamp, spacing = resample(imgs_to_process, patient, [1,1,1])\nprint(f'Shape after resampling: {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    return verts, faces\n\n\ndef plotly_3d(verts, faces):\n    x,y,z = zip(*verts) \n    \n    print(\"Drawing\")\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 = ff.create_trisurf(x=x, y=y, z=z, plot_edges=False,\n                        colormap=colormap,\n                        simplices=faces,\n                        backgroundcolor='rgb(64, 64, 64)',\n                        title=\"Interactive Visualization\")\n    iplot(fig)\n\ndef plt_3d(verts, faces):\n    print(\"Drawing\")\n    x,y,z = zip(*verts) \n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n\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_axis_bgcolor((0.7, 0.7, 0.7))\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, threshold = 350) #350 previously default value\nplt_3d(v, f)\n# plotly_3d(v, f)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def air_removal_mask(dilation):\n    labels = measure.label(dilation)\n    label_vals = np.unique(labels)\n    if labels[0,0] == labels[-1, -1]:\n        upper_cut = (labels==labels[0,0])\n        mask = np.abs(upper_cut*1 -1) \n    else:\n        upper_cut = (labels == labels[0,0])\n        lower_cut = (labels == labels[-1,-1])\n        mask = np.abs((upper_cut + lower_cut )*1 -1)      \n    return mask\n\n#Standardize the pixel values\ndef make_lungmask(img, display=False):\n    row_size= img.shape[0]\n    col_size = img.shape[1]\n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    # Find the average pixel value near the lungs\n    # to renormalize washed out images\n    middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n    mean = np.mean(middle)  \n    max = np.max(img)\n    min = np.min(img)\n    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    #\n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    #\n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n\n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n    #\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    #\n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation [10,10]\n    \n    mask = dilation.astype('int16')*air_removal_mask(dilation)\n    \n    if (display):\n        fig, ax = plt.subplots(3, 2, figsize=[12, 12])\n        ax[0, 0].set_title(\"Original\")\n        ax[0, 0].imshow(img, cmap='gray')\n        ax[0, 0].axis('off')\n        \n        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\n        \n        ax[1, 0].set_title(\"After Erosion and Dilation\")\n        ax[1, 0].imshow(dilation, cmap='gray')\n        ax[1, 0].axis('off')\n        \n        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        \n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\n        \n        ax[2, 1].set_title(\"Apply Mask on Original\")\n        ax[2, 1].imshow(mask*img, cmap='gray')\n        ax[2, 1].axis('off')\n        \n        plt.show()\n    return mask*img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img = imgs_after_resamp[120]\noutput = make_lungmask(img, display=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm import tqdm\nmasked_lung = []\n\nfor img in tqdm(imgs_after_resamp):\n    masked_lung.append(make_lungmask(img))\n\nsample_stack(masked_lung, show_every=6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.save(f'maskedimages_{patient_id}.npy', masked_lung)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\n\nPATH = \"../input/rsna-str-pulmonary-embolism-detection/\"\ntrain = pd.read_csv(PATH + \"train.csv\")\nsub = pd.read_csv(PATH + \"sample_submission.csv\")\n\nfeats = list(train.columns[3:5])+list(train.columns[8:12])+list(train.columns[13:17])\nmeans = train[feats].mean().to_dict()\n\n\nsub['label'] = 0.28\nfor feat in means.keys():\n    sub.loc[sub.id.str.contains(feat, regex=False), 'label'] = means[feat]\n    \nsub.to_csv('submission.csv', index = False)","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}