{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"# Regular Imports\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 images\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\")\nfrom pathlib import Path","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import train + test data\ntrain = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/train.csv\")\nINPUT = Path(\"../input/osic-pulmonary-fibrosis-progression/\")\ntest = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/test.csv\")\n\n# Train len\nprint(\"Total Recordings in Train Data: {:,}\".format(len(train)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_n_dicom_df(train):\n    df_list = []\n    for patient_id in train['Patient'].unique():\n        patient_dir = INPUT / 'train' / patient_id\n        path_list = list(patient_dir.glob(\"*\"))\n        n_dicom = len(path_list)\n        n_list = [int(str(i).split('/')[-1].split('.')[0]) for i in path_list]\n        sort_n_list = sorted(n_list)\n        tmp_df = pd.DataFrame({'Patient': [patient_id],\n                               'n_dicom': [n_dicom],\n                               'n_list': [sort_n_list]})\n        df_list.append(tmp_df)\n        \n    n_dicom_df = pd.concat(df_list, sort=False)\n    return n_dicom_df.reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_dicom_df = get_n_dicom_df(train)\nn_dicom_df.to_csv('n_dicom_df.csv', index=False)\nn_dicom_df.head(15)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(n_dicom_df['n_dicom'], bins=20)\nplt.title('Number of dicom per patient');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df1_styler = train.head().style.set_table_attributes(\"style='display:inline'\").set_caption('Head Train Data')\ndf2_styler = test.style.set_table_attributes(\"style='display:inline'\").set_caption('Test Data (rest Hidden)')\n\ndisplay_html(df1_styler._repr_html_() + df2_styler._repr_html_(), raw=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"There are {} unique patients in Train Data.\".format(len(train[\"Patient\"].unique())), \"\\n\")\n\n# Recordings per Patient\ndata = train.groupby(by=\"Patient\")[\"Weeks\"].count().reset_index(drop=False)\n# Sort by Weeks\ndata = data.sort_values(['Weeks']).reset_index(drop=True)\nprint(\"Minimum number of entries are: {}\".format(data[\"Weeks\"].min()), \"\\n\" +\n      \"Maximum number of entries are: {}\".format(data[\"Weeks\"].max()))\n\n# Plot\nplt.figure(figsize = (16, 6))\np = sns.barplot(data[\"Patient\"], data[\"Weeks\"], color=custom_colors[2])\n\nplt.title(\"Number of Entries per Patient\", fontsize = 17)\nplt.xlabel('Patient', fontsize=14)\nplt.ylabel('Frequency', fontsize=14)\n\np.axes.get_xaxis().set_visible(False);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Select unique bio info for the patients\ndata = train.groupby(by=\"Patient\")[[\"Patient\", \"Age\", \"Sex\", \"SmokingStatus\"]].first().reset_index(drop=True)\n\n# Figure\nf, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (16, 6))\n\na = sns.distplot(data[\"Age\"], ax=ax1, color=custom_colors[1], hist=False, kde_kws=dict(lw=6, ls=\"--\"))\nb = sns.countplot(data[\"Sex\"], ax=ax2, palette=custom_colors[2:4])\nc = sns.countplot(data[\"SmokingStatus\"], ax=ax3, palette = custom_colors[4:7])\n\na.set_title(\"Patient Age Distribution\", fontsize=16)\nb.set_title(\"Sex Frequency\", fontsize=16)\nc.set_title(\"Smoking Status\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Min FVC value: {:,}\".format(train[\"FVC\"].min()), \"\\n\" +\n      \"Max FVC value: {:,}\".format(train[\"FVC\"].max()), \"\\n\" +\n      \"\\n\" +\n      \"Min Percent value: {:.4}%\".format(train[\"Percent\"].min()), \"\\n\" +\n      \"Max Percent value: {:.4}%\".format(train[\"Percent\"].max()))\n\n# Figure\nf, (ax1, ax2) = plt.subplots(1, 2, figsize = (16, 6))\n\na = sns.distplot(train[\"FVC\"], ax=ax1, color=custom_colors[6], hist=False, kde_kws=dict(lw=6, ls=\"--\"))\nb = sns.distplot(train[\"Percent\"], ax=ax2, color=custom_colors[4], hist=False, kde_kws=dict(lw=6, ls=\"-.\"))\n\na.set_title(\"FVC Distribution\", fontsize=16)\nb.set_title(\"Percent Distribution\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_dicom_df['n_dicom'].value_counts().head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Minimum no. weeks before CT: {}\".format(train['Weeks'].min()), \"\\n\" +\n      \"Maximum no. weeks after CT: {}\".format(train['Weeks'].max()))\n\nplt.figure(figsize = (16, 6))\n\na = sns.distplot(train['Weeks'], color=custom_colors[3], hist=False, kde_kws=dict(lw=8, ls=\"--\"))\nplt.title(\"Number of weeks before/after the CT scan\", fontsize = 16)\nplt.xlabel(\"Weeks\", fontsize=14);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Correlations between Variables"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Compute Correlation\ncorr1, _ = pearsonr(train[\"FVC\"], train[\"Percent\"])\ncorr2, _ = pearsonr(train[\"FVC\"], train[\"Age\"])\ncorr3, _ = pearsonr(train[\"Percent\"], train[\"Age\"])\nprint(\"Pearson Corr FVC x Percent: {:.4}\".format(corr1), \"\\n\" +\n      \"Pearson Corr FVC x Age: {:.0}\".format(corr2), \"\\n\" +\n      \"Pearson Corr Percent x Age: {:.2}\".format(corr3))\n\n# Figure\nf, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (16, 6))\n\na = sns.scatterplot(x = train[\"FVC\"], y = train[\"Percent\"], palette=[custom_colors[2], custom_colors[6]],\n                    hue = train[\"Sex\"], style = train[\"Sex\"], s=100, ax=ax1)\n\nb = sns.scatterplot(x = train[\"FVC\"], y = train[\"Age\"], palette=[custom_colors[2], custom_colors[6]],\n                    hue = train[\"Sex\"], style = train[\"Sex\"], s=100, ax=ax2)\n\nc = sns.scatterplot(x = train[\"Percent\"], y = train[\"Age\"], palette=[custom_colors[2], custom_colors[6]],\n                    hue = train[\"Sex\"], style = train[\"Sex\"], s=100, ax=ax3)\n\na.set_title(\"Correlation between FVC and Percent\", fontsize = 16)\na.set_xlabel(\"FVC\", fontsize = 14)\na.set_ylabel(\"Percent\", fontsize = 14)\n\nb.set_title(\"Correlation between FVC and Age\", fontsize = 16)\nb.set_xlabel(\"FVC\", fontsize = 14)\nb.set_ylabel(\"Age\", fontsize = 14)\n\nc.set_title(\"Correlation between Percent and Age\", fontsize = 16)\nc.set_xlabel(\"Percent\", fontsize = 14)\nc.set_ylabel(\"Age\", fontsize = 14);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Figure\nf, (ax1, ax2) = plt.subplots(1,2, figsize = (16, 6))\n\na = sns.barplot(x = train[\"SmokingStatus\"], y = train[\"FVC\"], ax=ax1, palette=custom_colors[0:4])\nb = sns.barplot(x = train[\"SmokingStatus\"], y = train[\"Percent\"], ax=ax2, palette=custom_colors[4:7])\n\na.set_title(\"Mean FVC per Smoking Status\", fontsize=16)\nb.set_title(\"Mean Perc per Smoking Status\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evolution of FVC in Patients"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create Time variable to count in ascending order the times the Patient has done a check in FVC\ndata_time = train.groupby(by=\"Patient\")[\"Weeks\"].count().reset_index()\ntrain[\"Time\"] = 0\n\nfor patient, times in zip(data_time[\"Patient\"], data_time[\"Weeks\"]):\n    train.loc[train[\"Patient\"] == patient, 'Time'] = range(1, times+1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"min_fvc = train[train[\"Time\"] == 1][[\"Patient\", \"FVC\"]].reset_index(drop=True)\n\nidx = train.groupby([\"Patient\"])[\"Weeks\"].transform(max) == train[\"Weeks\"]\nmax_fvc = train[idx][[\"Patient\", \"FVC\"]].reset_index(drop=True)\n\n# Compute difference and select only top patients with biggest difference\ndata = pd.merge(min_fvc, max_fvc, how=\"inner\", on=\"Patient\")\ndata[\"Dif\"] = data[\"FVC_x\"] - data[\"FVC_y\"]\n\n# Select only top n\nl = list(data.sort_values(\"Dif\", ascending=False).head(100)[\"Patient\"])\nx = train[train[\"Patient\"].isin(l)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\na = sns.lineplot(x = x[\"Time\"], y = x[\"FVC\"], hue = x[\"Patient\"], legend=False,\n                 palette=sns.color_palette(\"GnBu_d\", 100), size=1)\n\nplt.title(\"Patient FVC decrease on Weeks\", fontsize = 16)\nplt.xlabel(\"Weeks\", fontsize=14)\nplt.ylabel(\"FVC\", fontsize=14);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Number of CT scans per Patient"},{"metadata":{"trusted":true},"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);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualise DICOM"},{"metadata":{"trusted":true},"cell_type":"code","source":"class bcolors:\n    OKBLUE = '\\033[96m'\n    OKGREEN = '\\033[92m'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from PIL import Image\nfrom IPython.display import Image as show_gif\nimport scipy.misc\nimport matplotlib","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_gif(number_of_CT = 87):\n    \"\"\"Picks a patient at random and creates a GIF with their CT scans.\"\"\"\n    \n    # Select one of the patients\n    # patient = \"ID00007637202177411956430\"\n    patient = train[train[\"CT_number\"] == number_of_CT].sample(random_state=1)[\"Patient\"].values[0]\n    \n    # === READ IN .dcm FILES ===\n    patient_dir = \"../input/osic-pulmonary-fibrosis-progression/train/\" + patient\n    datasets = []\n\n    # First Order the files in the dataset\n    files = []\n    for dcm in list(os.listdir(patient_dir)):\n        files.append(dcm) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Read in the Dataset from the Patient path\n    for dcm in files:\n        path = patient_dir + \"/\" + dcm\n        datasets.append(pydicom.dcmread(path))\n        \n        \n    # === SAVE AS .png ===\n    # Create directory to save the png files\n    if os.path.isdir(f\"png_{patient}\") == False:\n        os.mkdir(f\"png_{patient}\")\n\n    # Save images to PNG\n    for i in range(len(datasets)):\n        img = datasets[i].pixel_array\n        matplotlib.image.imsave(f'png_{patient}/img_{i}.png', img)\n        \n        \n    # === CREATE GIF ===\n    # First Order the files in the dataset (again)\n    files = []\n    for png in list(os.listdir(f\"../working/png_{patient}\")):\n        files.append(png) \n    files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n\n    # Create the frames\n    frames = []\n\n    # Create frames\n    for file in files:\n    #     print(\"../working/png_images/\" + name)\n        new_frame = Image.open(f\"../working/png_{patient}/\" + file)\n        frames.append(new_frame)\n\n    # Save into a GIF file that loops forever\n    frames[0].save(f'gif_{patient}.gif', format='GIF',\n                   append_images=frames[1:],\n                   save_all=True,\n                   duration=200, loop=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"create_gif(number_of_CT=12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_gif(filename=\"./gif_ID00165637202237320314458.gif\", format='png', width=400, height=400)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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    \nimgs = []\nfor data in datasets:\n    img = data.pixel_array\n    imgs.append(img)\n    \n    \n# Show masks\nfig=plt.figure(figsize=(16, 6))\ncolumns = 10\nrows = 3\n\nfor i in range(1, columns*rows +1):\n    img = make_lungmask(datasets[i-1].pixel_array)\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(img, cmap=\"gray\")\n    plt.title(i, fontsize = 9)\n    plt.axis('off');","execution_count":null,"outputs":[]},{"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}