{"cells":[{"metadata":{},"cell_type":"markdown","source":"As we have noticed that each of the patient has almost 200+ slices of CT-Scan. While a single CT-Scan slice can carry a lot of information, however, using a single CT-Scan to identify any disease or condition undermines the true purpose of having a CT-Scan. With this view in mind, I intend to sort the slices in the correct order and then producing a 3D volumetric data so that it is possible to get a first hand look at the volumetric information. \n\nPlease upvote the notebook if you think it helped  you somehow.. :) "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nimport numpy as np, pandas as pd\nimport pydicom, imageio, os\nimport matplotlib.pyplot as plt\nfrom IPython import display\nfrom IPython.display import HTML\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":{},"cell_type":"markdown","source":"At first let's list all the patient IDs so that we can perform the same operation for each of the patient and save the result. "},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv('../input/rsna-str-pulmonary-embolism-detection/train.csv')\ntest_df = pd.read_csv('../input/rsna-str-pulmonary-embolism-detection/test.csv')\n\n\npatient_id_list = os.listdir('../input/rsna-str-pulmonary-embolism-detection/train/')\nprint(f'Total number of patient(experiment): {len(patient_id_list)}\\nFirst 5 patient IDs')\npatient_id_list[:5]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now I will scan all the files of a single patient(experiment). We will just produce a list of file-paths of that particular experiment. Here I have chosen id 6900 because it has both some positive slice and negative slices for Pulmonary Embolism (PE). "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"index = patient_id_list.index('6897fa9de148')\npatient_id = patient_id_list[index]\npatient_folder = f'../input/rsna-str-pulmonary-embolism-detection/train/{patient_id}/'\npatient_image_paths = glob(patient_folder + '/*/*.dcm')\n[neg, pos] = train_df[train_df.StudyInstanceUID == patient_id].pe_present_on_image.value_counts().values\n\n\n# Print out the first 5 file names to verify we're in the right folder.\nprint (f'Total of {len(patient_image_paths)} DICOM images.' )\nprint(f'The experiment {patient_id} has {pos} positive and {neg} negative examples for pe_present_on_image\\nFirst 5 filenames:')\n\npatient_image_paths[:5]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now I know that the experiment has both positive and negative samples. So it will be easier to visualize the image plotted. \n\n\nLet's make some utility functions to deal with the dicom imges. If you have problem understanding you can have a look at my other notebooks on this same competition. \n"},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_image_paths[0][-16:-4]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_slice(paths):\n    slices = [pydicom.read_file(path) for path in paths]\n#     labels = [train_df[train_df.SOPInstanceUID == path[-16:-4]].pe_present_on_image.values for path in patient_image_paths]\n#     labels = np.array(labels).squeeze()\n    slices.sort(key = lambda x: int(x.InstanceNumber), reverse = False)\n#     labels.sort(key = lambda x: int(x.InstanceNumber), reverse = False)\n    \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    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\ndef transform_to_hu(slices):\n    images = np.stack([file.pixel_array for file in slices])\n    images = images.astype(np.int16)\n\n    # convert ouside pixel-values to air:\n    # I'm using <= -1000 to be sure that other defaults are captured as well\n    images[images <= -1000] = 0\n    \n    # convert to HU\n    for n in range(len(slices)):    \n        intercept = slices[n].RescaleIntercept\n        slope = slices[n].RescaleSlope\n        if slope != 1:\n            images[n] = slope * images[n].astype(np.float64)\n            images[n] = images[n].astype(np.int16)      \n        images[n] += np.int16(intercept)\n    return np.array(images, dtype=np.int16)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stacked_dicoms = load_slice(patient_image_paths)\nstacked_patient_pixels = transform_to_hu(stacked_dicoms)\n\ndef sample_stack(stack, rows=6, cols=6, start_with=10, show_every=3):\n    fig,ax = plt.subplots(rows,cols,figsize=[20,22])\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\n     \nsample_stack(stacked_patient_pixels, show_every = 3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now let us have a look at the images with continuous volumetric information. "},{"metadata":{"trusted":true},"cell_type":"code","source":"imageio.mimsave(f'./{patient_id}.gif', stacked_patient_pixels, duration=0.1)\ndisplay.Image(f'./{patient_id}.gif', format='png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Resampling\nAlthough we have each individual slices, it is not immediately clear how thick each slice is. Fortunately, this is in the DICOM header."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Slice Thickness: {stacked_dicoms[0].SliceThickness}')\nprint(f'Pixel Spacing (row, col): ({stacked_dicoms[0].PixelSpacing[0]}, {stacked_dicoms[0].PixelSpacing[1]})')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This means we have 2.0 mm slices, and each voxel represents 0.7 mm.\n\nBecause a CT slice is typically reconstructed at 512 x 512 voxels, each slice represents approximately 370 mm of data in length and width.\n\nUsing the metadata from the DICOM we can figure out the size of each voxel as the slice thickness. In order to display the CT in 3D isometric form (which we will do below), and also to compare between different scans, it would be useful to ensure that each slice is resampled in 1x1x1 mm pixels and slices."},{"metadata":{"trusted":true},"cell_type":"code","source":"np.array([float(stacked_dicoms[0].SliceThickness), float(stacked_dicoms[0].PixelSpacing[0]), float(stacked_dicoms[0].PixelSpacing[0])])","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 = np.array([float(scan[0].SliceThickness), \n                        float(scan[0].PixelSpacing[0]), \n                        float(scan[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: {stacked_patient_pixels.shape}')\nimgs_after_resamp, spacing = resample(stacked_patient_pixels, stacked_dicoms, [1,1,1])\nprint(f'Shape after resampling: {imgs_after_resamp.shape}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let us make a mask that will remvoe the upper and lower part of the slices which are mostly air. "},{"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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now let us make another mask that will only the lung part of the slices. As a result it will be better for the ML models to concentrate on what's important. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#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\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's look at what we have done. "},{"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":"imageio.mimsave(f'segmented{patient_id}.gif', masked_lung, duration=0.1)\ndisplay.Image(f'segmented{patient_id}.gif', format='png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Saving to npy format\nThis is a good time to save the new data set to disk so we don't have to reprocess the stack every time. And this is not final. I will make a utility function from this and from that I will convert the whole dataset into this format. "},{"metadata":{"trusted":true},"cell_type":"code","source":"np.save(f'resampled_masked_lung_{patient_id}.npy', masked_lung)\nnp.save(f'stacked_slices_{patient_id}.npy', stacked_patient_pixels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"![](https://www.clipartmax.com/png/middle/265-2655834_work-in-progress-icon.png)\n\n\n\n\n\n### In the meantime check out my other ongoing works in this same competition: \n💥 [RSNA-STR Pulmonary Embolism [Dummy Sub]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-dummy-sub)<br>\n💥 [CT-Scans, DICOM files, Windowing Explained](https://www.kaggle.com/redwankarimsony/ct-scans-dicom-files-windowing-explained)<br>\n💥 [RSNA-STR-PE [Gradient & Sigmoid Windowing]](https://www.kaggle.com/redwankarimsony/rsna-str-pe-gradient-sigmoid-windowing)<br>\n💥 [RSNA-STR [✔️3D Stacking ✔️3D Plot ✔️Segmentation]](https://www.kaggle.com/redwankarimsony/rsna-str-3d-stacking-3d-plot-segmentation/edit/run/42517982)<br>\n💥 [RSNA-STR [DICOM 👉 GIF 👉 npy]](https://www.kaggle.com/redwankarimsony/rsna-str-dicom-gif-npy)<br>\n💥 [RSNA-STR Pulmonary Embolism [EDA]](https://www.kaggle.com/redwankarimsony/rsna-str-pulmonary-embolism-eda)<br>\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}