{"cells":[{"metadata":{},"cell_type":"markdown","source":"# About Pulmonary Fibrosis\nPulmonary fibrosis is a lung disease that occurs when lung tissue becomes damaged and scarred. This thickened, stiff tissue makes it more difficult for your lungs to work properly. Over time, the scar tissue can destroy the normal lung and make it hard for oxygen to get into your blood.\n\n\n![](https://www.pulmonaryfibrosis.org/images/default-source/default-album/normal-and-impaired-gas-exchange.png?sfvrsn=c3b0918d_0)\n\n## Symptoms and Causes\n\nThere are five main categories of identifiable causes of pulmonary fibrosis: **Drug-induced, Radiation-induced, Environmental, Autoimmune, and Occupational**. Sometimes it can be challenging for doctors to figure out what causes PF. PF of unknown cause is called **“idiopathic”**.\n\n\"\"\"\"\n\n\n## Treatment Options\n\nAs given in the description, current methods make fibrotic lung diseases **difficult to treat**, even with access to a chest CT scan. In addition, the wide range of varied prognoses create issues organizing clinical trials.\n\n\n## Prognosis\n\n**Prognosis** is a medical term for predicting the likely or expected development of a disease, including whether the signs and symptoms will improve or worsen (and how quickly) or remain stable over time; expectations of quality of life, such as the ability to carry out daily activities; the potential for complications and associated health issues; and the likelihood of survival (including life expectancy). \n\n\n**Reference:** [https://www.pulmonaryfibrosis.org/life-with-pf/about-pf](http://).\n\n## Evaluation Metric\n\nFor each true FVC measurement, we will predict both an FVC and a confidence measure (standard deviation σ). The metric is computed as:\n\n\\begin{equation} \n\\sigma_{clipped} = max(\\sigma, 70)\\\\\n\\Delta = min ( |FVC_{true} - FVC_{predicted}|, 1000 )\\\\\nmetric = - \\frac{\\sqrt{2} \\Delta}{\\sigma_{clipped}} - \\ln ( \\sqrt{2} \\sigma_{clipped} )\n\\end{equation}","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"* [Understand the evaluation metric](https://www.kaggle.com/rohanrao/osic-understanding-laplace-log-likelihood).\n* [Understanding CT scans](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/167085).\n\nThanks to Twinkle Khanna for work https://www.kaggle.com/twinkle0705/your-starter-notebook-for-osic and \nGunes Evitan for work https://www.kaggle.com/gunesevitan/osic-pulmonary-fibrosis-progression-eda making this notebook possible.","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import 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\nimport matplotlib.lines as mlines\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,probplot, mode\nimport tqdm\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\")\n\n\n# Set Color Palettes for the notebook\ncustom_colors = ['#74a09e','#86c1b2','#98e2c6','#f3c969','#f2a553', '#d96548', '#c14953']\nsns.palplot(sns.color_palette(custom_colors))\n\n# Set Style\nsns.set_style(\"whitegrid\")\nsns.despine(left=True, bottom=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.rc('xtick',labelsize=11)\nplt.rc('ytick',labelsize=11)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# The Metadata \n\n### What we need to predict\n* FVC - final 3 values for each patient (only these will be used for the final score).\n* Confidence Value — [a thread about what it is here](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/166753). High value in confidence means that you are off a lot from the actual FVC, while a very low (or even 0) confidence means you're very sure on the FVC.\n\n> ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import train + test data\ntrain = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/train.csv\")\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":{"_kg_hide-input":true,"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":{},"cell_type":"markdown","source":"## Check for missing values","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Q: Are there any missing values?\", \"\\n\" +\n \"A: {}\".format(train.isnull().values.any()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## How many unique patients?","execution_count":null},{"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":{},"cell_type":"markdown","source":"## Patients bio: who are they?","execution_count":null},{"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 = (21, 7))\n\na = sns.distplot(data[\"Age\"], ax=ax1, color=custom_colors[1], hist=True, kde_kws=dict(lw=0))\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":{},"cell_type":"markdown","source":"## FVC & Percent\n\n\n

FVC the recorded lung capacity in ml (how much air you can exhale in a maximal forced expiration effort).

\n

Percent a computed field which approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics.

\n\n* **FVC**: most values lie between 1,000 and 5,000 ml. There are also some very high outliers above 5,000 and little values that lie below 1,000.\n* **Percent**: more than ~80% of the patients scored below 100%.","execution_count":null},{"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, ax3) = plt.subplots(1, 3, figsize = (21, 7))\n\na = sns.distplot(train[\"FVC\"], ax=ax1, color=custom_colors[6],kde_kws=dict(lw=3, ls=\"-\"))\nb = probplot(train['FVC'], plot=ax2)\nc = sns.distplot(train[\"Percent\"], ax=ax3, color=custom_colors[4], kde_kws=dict(lw=3, ls=\"-\"))\n\na.set_title(\"FVC Distribution\", fontsize=16)\nc.set_title(\"Percent Distribution\", fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Weeks","execution_count":null},{"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=True, kde_kws=dict(lw=3, 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\n\n\n* There is high correlation between FVC and Percent: when the volume of air increases, the Percent increases as well (when you exhale more, you get closer more to 100%).\n* There is no correlation between FVC/Percent and Age, meaning that Age has no influence on the volume of exhaled air.","execution_count":null},{"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 = (21, 7))\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)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> This is very weird: FVC and Percent are the highest for people that still smoke and the lowest for people that never smoked. However, we need to keep in mind that the percentage of people that still smoke is very low. So, we can't conclude that if a person smokes it's highly likely that will have a high FVC.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 6))\n\nsns.barplot(x=train.groupby('Patient')['SmokingStatus'].first().value_counts().index, y=train.groupby('Patient')['SmokingStatus'].first().value_counts(), palette=custom_colors[0:3])\npercentages = [(count / train.groupby('Patient')['SmokingStatus'].first().value_counts().sum() * 100).round(2) for count in train.groupby('Patient')['SmokingStatus'].first().value_counts()]\n\nplt.ylabel('')\nplt.xticks(np.arange(3), [f'Ex-smoker (%{percentages[0]})', f'Never smoked (%{percentages[1]})', f'Currently Smokes (%{percentages[2]})'])\nplt.tick_params(axis='x', labelsize=12)\nplt.tick_params(axis='y', labelsize=12)\nplt.title('SmokingStatus Counts in Training Set', size=15, pad=15)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Tabular Data Correlations","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(6, 6), dpi=150)\n\nsns.heatmap(train[['Weeks', 'FVC', 'Percent', 'Age']].corr(), annot=True, square=True, cmap='summer', annot_kws={'size': 12}, fmt='.2f') \n\nplt.tick_params(axis='x', labelsize=11, rotation=0)\nplt.tick_params(axis='y', labelsize=11, rotation=0)\nplt.title('Tabular Data Features Correlations')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# DICOM Data","execution_count":null},{"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":{"trusted":true},"cell_type":"code","source":"slice_counts = np.array([len(os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{directory}')) for directory in os.listdir('../input/osic-pulmonary-fibrosis-progression/train')])\n\nprint(f'Number of Image Slices in Training Set\\n{\"-\" * 38}')\nprint(f'Mean Slice Count: {slice_counts.mean():.6} - Median Slice Count: {int(np.median(slice_counts))} - Total Slice Count: {slice_counts.sum()}')\nprint(f'Min Slice Count: {slice_counts.min()} - Max Slice Count: {slice_counts.max()}')\n\nfig = plt.figure(figsize=(8, 2), dpi=150)\nax = sns.countplot(slice_counts, palette='autumn')\n\nfor idx, label in enumerate(ax.get_xticklabels()):\n if idx % 10 == 0:\n label.set_visible(True)\n else:\n label.set_visible(False)\n\nplt.ylabel('')\nplt.tick_params(axis='x', labelsize=8)\nplt.tick_params(axis='y', labelsize=8)\nplt.title('Number of Image Slices in Training Set', size=10)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Number of CT scans per Patient","execution_count":null},{"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":"## Metadata Extraction","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_metadata(patient_name):\n \n patient_directory = [pydicom.dcmread(f'../input/osic-pulmonary-fibrosis-progression/train/{patient_name}/{s}') for s in os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{patient_name}')]\n \n try:\n patient_directory.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n except AttributeError:\n patient_directory.sort(key=lambda s: int(s.InstanceNumber)) \n \n rows = patient_directory[0].Rows\n columns = patient_directory[0].Columns\n slices = len(patient_directory)\n \n slice_thicknesses = []\n slice_spacings = [] \n rescale_slopes = []\n rescale_intercepts = []\n \n for i, s in enumerate(patient_directory):\n slice_thicknesses.append(s.SliceThickness)\n rescale_slopes.append(s.RescaleSlope)\n rescale_intercepts.append(s.RescaleIntercept)\n try:\n slice_spacings.append(s.SpacingBetweenSlices)\n except AttributeError:\n slice_spacings.append(np.nan)\n \n train.loc[train['Patient'] == patient_name, 'Rows'] = rows\n train.loc[train['Patient'] == patient_name, 'Columns'] = columns\n train.loc[train['Patient'] == patient_name, 'Slices'] = slices\n train.loc[train['Patient'] == patient_name, 'SliceThickness'] = mode(slice_thicknesses)[0][0]\n train.loc[train['Patient'] == patient_name, 'SliceSpacing'] = mode(slice_spacings)[0][0] \n train.loc[train['Patient'] == patient_name, 'RescaleSlope'] = mode(rescale_slopes)[0][0]\n train.loc[train['Patient'] == patient_name, 'RescaleIntercept'] = mode(rescale_intercepts)[0][0]\n \nfor patient in tqdm.tqdm(train['Patient'].unique()):\n get_metadata(patient)\n \ntrain['Rows'] = train['Rows'].astype(np.uint16)\ntrain['Columns'] = train['Columns'].astype(np.uint16)\ntrain['Slices'] = train['Slices'].astype(np.uint16)\ntrain['SliceShape'] = train['Rows'].astype(str) + 'x' + train['Columns'].astype(str)\ntrain['SliceThickness'] = train['SliceThickness'].astype(np.float32)\ntrain['SliceSpacing'] = train['SliceSpacing'].astype(np.float32)\ntrain['RescaleSlope'] = train['RescaleSlope'].astype(np.float32)\ntrain['RescaleIntercept'] = train['RescaleIntercept'].astype(np.float32)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Slice Shapes\nPatients have different slice shapes but shapes are consistent in patients' directories. A patient can only have one unique slice shape. There are 10 unique 2D slice shapes in training set and the most common shapes are (512, 512) and (768, 768). There are some unusual slice shapes which are very rare. Those unusual slice shapes belong to 1 or 2 different patients.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(9, 3), dpi=100)\nsns.barplot(x=train.groupby('Patient')['SliceShape'].first().value_counts().values, y=train.groupby('Patient')['SliceShape'].first().value_counts().index, palette='summer')\n\nplt.xlabel('Patients', size=15)\nplt.ylabel('')\nplt.tick_params(axis='x', labelsize=10)\nplt.tick_params(axis='y', labelsize=10)\nplt.title(f'Training Set Slice Shapes', size=15)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Cropping images\nThe function crop_slice is used for cropping the frames. It drops all zero rows and columns in the given numpy array. After the cropping operation, all of the slices with unusual shapes are changed to 512x512.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def crop_slice(s):\n\n \"\"\"\n Crop frames from slices\n\n Parameters\n ----------\n s : numpy array, shape = (Rows, Columns)\n numpy array of slices with frame\n\n Returns\n -------\n s_cropped : numpy array, shape = (Rows - All Zero Rows, Columns - All Zero Columns)\n numpy array after the all zero rows and columns are dropped\n \"\"\"\n\n s_cropped = s[~np.all(s == 0, axis=1)]\n s_cropped = s_cropped[:, ~np.all(s_cropped == 0, axis=0)]\n return s_cropped","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(figsize=(20, 9), nrows=2, ncols=5)\n\nfor i, patient_name in enumerate(train.groupby('SliceShape')['Patient'].first()): \n \n patient_directory = [pydicom.dcmread(f'../input/osic-pulmonary-fibrosis-progression/train/{patient_name}/{s}') for s in os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{patient_name}')]\n patient_directory.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n first_slice_cropped = crop_slice(patient_directory[0].pixel_array)\n if first_slice_cropped.shape[0] != 512 and first_slice_cropped.shape[1] != 512:\n first_slice_cropped_resized = cv2.resize(first_slice_cropped, (512, 512), interpolation=cv2.INTER_NEAREST)\n else:\n first_slice_cropped_resized = first_slice_cropped\n old_aspect_ratio = train.groupby(\"SliceShape\")[\"Patient\"].first().index[i]\n new_aspect_ratio = f'{first_slice_cropped_resized.shape[0]}x{first_slice_cropped_resized.shape[1]}'\n \n if i > 4:\n i -= 5\n j = 1\n else:\n j = 0\n \n axes[j][i].imshow(first_slice_cropped_resized, cmap=plt.cm.bone)\n\n axes[j][i].tick_params(axis='x', labelsize=14)\n axes[j][i].tick_params(axis='y', labelsize=14)\n axes[j][i].set_title(f'{old_aspect_ratio} -> {new_aspect_ratio}', size=14)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## DICOM Lung Mask\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: \n* be sure you don't exclude important parts (like lesions)","execution_count":null},{"metadata":{"trusted":true},"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(imgrow_size/5 and B[2] **📌Remember**: Lungs are quite visible, but on images at the beginning and towards the end the lung dissapears completely.","execution_count":null},{"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":{"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}