{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\nimport os\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport matplotlib.image as mpimg\nfrom tabulate import tabulate\nimport missingno as msno \nfrom IPython.display import display_html\nfrom PIL import Image\nimport gc\nimport cv2\nfrom scipy.stats import pearsonr\n\nimport pydicom # for DICOM \nfrom skimage.transform import resize\nimport copy\nimport re\n\n# Segmentation\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\nfrom plotly.tools import FigureFactory as FF\nfrom plotly.graph_objs import *\ninit_notebook_mode(connected=True) \n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\n\n","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# DICOM Data \n\n","metadata":{}},{"cell_type":"code","source":"# Create base director for Train .dcm files\ndirector = \"../input/osic-pulmonary-fibrosis-progression/train\"\n\n# Create path column with the path to each patient's CT\ntrain[\"Path\"] = director + \"/\" + train[\"Patient\"]\n\n# Create variable that shows how many CT scans each patient has\ntrain[\"CT_number\"] = 0\n\nfor k, path in enumerate(train[\"Path\"]):\n    train[\"CT_number\"][k] = len(os.listdir(path))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Minimum number of CT scans: {}\".format(train[\"CT_number\"].min()), \"\\n\" +\n      \"Maximum number of CT scans: {:,}\".format(train[\"CT_number\"].max()))\n\n# Scans per Patient\ndata = train.groupby(by=\"Patient\")[\"CT_number\"].first().reset_index(drop=False)\n# Sort by Weeks\ndata = data.sort_values(['CT_number']).reset_index(drop=True)\n\n# Plot\nplt.figure(figsize = (16, 6))\np = sns.barplot(data[\"Patient\"], data[\"CT_number\"], color=custom_colors[5])\nplt.axvline(x=85, color=custom_colors[2], linestyle='--', lw=3)\n\nplt.title(\"Number of CT Scans per Patient\", fontsize = 17)\nplt.xlabel('Patient', fontsize=14)\nplt.ylabel('Frequency', fontsize=14)\n\nplt.text(86, 850, \"Median=94\", fontsize=13)\n\np.axes.get_xaxis().set_visible(False);","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2 Visualise DICOM Info and Image\n\n> DICOM data can be extracted by using `pydicom.dcmread()`\n\n> [What color maps should you use in medical visualisation?](http://noeskasmit.com/colormaps-in-medical-visualization/)","metadata":{}},{"cell_type":"code","source":"path = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/19.dcm\"\ndataset = pydicom.dcmread(path)\n\nprint(bcolors.OKBLUE + \"Patient id.......:\", dataset.PatientID, \"\\n\" +\n      \"Modality.........:\", dataset.Modality, \"\\n\" +\n      \"Rows.............:\", dataset.Rows, \"\\n\" +\n      \"Columns..........:\", dataset.Columns)\n\nplt.figure(figsize = (7, 7))\nplt.imshow(dataset.pixel_array, cmap=\"plasma\")\nplt.axis('off');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430\"\ndatasets = []\n\n# First Order the files in the dataset\nfiles = []\nfor dcm in list(os.listdir(patient_dir)):\n    files.append(dcm) \nfiles.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n# Read in the Dataset\nfor dcm in files:\n    path = patient_dir + \"/\" + dcm\n    datasets.append(pydicom.dcmread(path))\n\n# Plot the images\nfig=plt.figure(figsize=(16, 6))\ncolumns = 10\nrows = 3\n\nfor i in range(1, columns*rows +1):\n    img = datasets[i-1].pixel_array\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(img, cmap=\"plasma\")\n    plt.title(i, fontsize = 9)\n    plt.axis('off');","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# show_gif(filename=\"./gif_ID00340637202287399835821.gif\", format='png', width=400, height=400)\n# show_gif(filename=\"./gif_ID00199637202248141386743.gif\", format='png', width=400, height=400)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.5 DICOM Lung Mask🎭\n\n<div class=\"alert alert-block alert-info\">\n<b>Reference:</b> <a link=\"https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\"> DICOM Processing Segmentation Visualization in Python</a>\n</div>\n\n**Mask on Lungs Purpose:**\n* Segmentation is part of the preprocessing method\n* Has the purpose of auto-detecting the boundaries surrounding a volume of interest (our case is the lungs)\n* Drawbacks: be sure you don't exclude important parts (like lesions)","metadata":{}},{"cell_type":"code","source":"# https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\n\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    \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    \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\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        ax[0, 1].set_title(\"Threshold\")\n        ax[0, 1].imshow(thresh_img, cmap='gray')\n        ax[0, 1].axis('off')\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        ax[1, 1].set_title(\"Color Labels\")\n        ax[1, 1].imshow(labels)\n        ax[1, 1].axis('off')\n        ax[2, 0].set_title(\"Final Mask\")\n        ax[2, 0].imshow(mask, cmap='gray')\n        ax[2, 0].axis('off')\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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### How does the mask work?","metadata":{}},{"cell_type":"code","source":"# Select a sample\npath = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/19.dcm\"\ndataset = pydicom.dcmread(path)\nimg = dataset.pixel_array\n\n# Masked image\nmask_img = make_lungmask(img, display=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4.1 Extracting Metadata from DICOM files\n\n> **📌Remember:** the `.get()` function is used overall for .dcm files that don't have all columns available. So, to not skip the entire file, we extract only what's available and fill in the rest with *None*.\n\n*Inspo from [Extract Metadata and Resize Notebook](https://www.kaggle.com/trsekhar123/nb-to-extract-metadata-and-resize-images-train)*","metadata":{}},{"cell_type":"code","source":"def get_observation_data(path):\n    \"\"\"Get information from the .dcm files.\n    path: complete path to the .dcm file\"\"\"\n\n    image_data = pydicom.read_file(path)\n    \n    # Dictionary to store the information from the image\n    observation_data = {\n        \"FileNumber\" : path.split(\"/\")[5],\n        \"Rows\" : image_data.get(\"Rows\"),\n        \"Columns\" : image_data.get(\"Columns\"),\n        \"PatientID\" : image_data.get(\"PatientID\"),\n        \"BodyPartExamined\" : image_data.get(\"BodyPartExamined\"),\n        \"RotationDirection\" : image_data.get(\"RotationDirection\"),\n        \"ConvolutionKernel\" : image_data.get(\"ConvolutionKernel\"),\n        \"PatientPosition\" : image_data.get(\"PatientPosition\"),\n        \"PhotometricInterpretation\" : image_data.get(\"PhotometricInterpretation\"),\n        \"Modality\" : image_data.get(\"Modality\"),\n        \"StudyInstanceUID\" : image_data.get(\"StudyInstanceUID\"),\n        \"PixelPaddingValue\" : image_data.get(\"PixelPaddingValue\"),\n        \"SamplesPerPixel\" : image_data.get(\"SamplesPerPixel\"),\n        \"BitsAllocated\" : image_data.get(\"BitsAllocated\"),\n        \"BitsStored\" : image_data.get(\"BitsStored\"),\n        \"HighBit\" : image_data.get(\"HighBit\"),\n        \"PixelRepresentation\" : image_data.get(\"PixelRepresentation\"),\n        \"RescaleType\" : image_data.get(\"RescaleType\"),\n    }\n\n    # Integer columns\n    int_columns = [\"SliceThickness\", \"KVP\", \"DistanceSourceToDetector\", \n        \"DistanceSourceToPatient\", \"GantryDetectorTilt\", \"TableHeight\", \n        \"XRayTubeCurrent\", \"GeneratorPower\", \"WindowCenter\", \"WindowWidth\", \n        \"SliceLocation\", \"RescaleIntercept\", \"RescaleSlope\"]\n    for k in int_columns:\n        observation_data[k] = int(image_data.get(k)) if k in image_data else None\n\n    # String columns\n    str_columns = [\"ImagePositionPatient\", \"ImageOrientationPatient\", \"ImageType\", \"PixelSpacing\"]\n    for k in str_columns:\n        observation_data[k] = str(image_data.get(k)) if k in image_data else None\n\n    \n    return observation_data","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test the function to see if it works","metadata":{}},{"cell_type":"code","source":"p = \"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/10.dcm\"\nexample = get_observation_data(p)\nexample\n# example","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get full paths for the images\npaths = []\nfor path in train[\"Path\"]:\n    for doc in os.listdir(path):\n        paths.append(path + \"/\" + doc)\n        ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create a dictionary with the data\nexceptions = 0\ndicts = []\n\nfor path in (paths):\n    # Get info in dict format\n    try:\n        d = get_observation_data(path)\n        dicts.append(d)\n    except Exception as e:\n        exceptions += 1\n        continue\n\n# Convert into a cudf dataframe\n# meta_train_data = cudf.DataFrame(data=dicts, columns=example.keys())\nmeta_train_data = pd.DataFrame(data=dicts, columns=example.keys())\n\n# Export information to a .csv\nmeta_train_data.to_csv(\"meta_train.csv\", index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Exceptions: {}\".format(exceptions))\nmeta_train_data.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}