{"cells":[{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')\nimport tqdm\nimport os\nimport gc\n\nimport numpy as np\nimport pandas as pd\npd.set_option('display.max_rows', 500)\npd.set_option('display.max_columns', 500)\npd.set_option('display.width', 1000)\n\nfrom scipy.stats import probplot, mode, linregress\n\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\nimport seaborn as sns\nfrom IPython.display import HTML\n\nimport gdcm\nimport pydicom\nimport cv2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ndf_test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\nprint(f'Training Set Shape = {df_train.shape} - Patients = {df_train[\"Patient\"].nunique()}')\nprint(f'Training Set Memory Usage = {df_train.memory_usage().sum() / 1024 ** 2:.2f} MB')\nprint(f'Test Set Shape = {df_test.shape} - Patients = {df_test[\"Patient\"].nunique()}')\nprint(f'Test Set Memory Usage = {df_test.memory_usage().sum() / 1024 ** 2:.2f} MB')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **1. Introduction**"},{"metadata":{},"cell_type":"markdown","source":"### **1.1. Training Set**\n\nThere are **1549** samples in `train.csv` and they belong to **176** different patients. All of the patients have their images acquired at Week 0, but all of the patients don't have FVC measured at Week 0 and their number of FVC measurements change from patient to patient.\n\nMost of the patients have FVC measurements at 9 different timesteps but this number can change between 6 and 10. Thus, number of FVC measurements are not consistent for different patients."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"training_sample_counts = df_train.rename(columns={'Weeks': 'Samples'}).groupby('Patient').agg('count')['Samples'].value_counts()\nprint(f'Training Set FVC Measurements Per Patient \\n{(\"-\") * 41}\\n{training_sample_counts}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **1.2. Test Set**\n\nThere are **5** samples in `test.csv` and they belong to **5** different patients. Those samples and patients also exist in `train.csv`. They are the last 5 patients of training set and the samples are first measurements of those patients. This is because `test.csv` is a placeholder and the real test set is hidden. The purpose of placeholder test set is showing the structure of real test set and testing submissions. When the notebook is submitted, it runs with the real test set.\n\nThe real test set will have more than 5 patients. However, it will only have a baseline CT scan and only the initial FVC measurement. You are asked to predict the final three FVC measurements for each patient, as well as a confidence value in your prediction. In order to avoid potential leakage, you are asked to predict every patient's FVC measurement for every possible week. Weeks prior to final three weeks are ignored in scoring."},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"df_test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **1.3. Sample Submission**\n\nThe simplest way of creating submissions is predicting `FVC` and `Confidence` for all test samples, then creating `Patient_Week` on test set and merging it to `sample_submission.csv`. This way submission file will have all of the predictions from test set regardless of their count. "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_submission = pd.read_csv( '../input/osic-pulmonary-fibrosis-progression/sample_submission.csv' )\ndf_submission","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **2. Laplace Log Likelihood**\n\nPredictions are evaluated with a modified version of the Laplace Log Likelihood. For each sample in test set, an `FVC` and a `Confidence` measure (standard deviation σ) has to be predicted.\n\n`Confidence` values smaller than 70 are clipped.\n\n$\\large \\sigma_{clipped} = max(\\sigma, 70),$\n\nErrors greater than 1000 are also clipped in order to avoid large errors.\n\n$\\large \\Delta = min ( |FVC_{true} - FVC_{predicted}|, 1000 ),$\n\nThe metric is defined as:\n\n$\\Large metric = -   \\frac{\\sqrt{2} \\Delta}{\\sigma_{clipped}} - \\ln ( \\sqrt{2} \\sigma_{clipped} ).$\n\nThe metric is implemented below and theoretical minimum score is calculated. If everything is predicted perfectly (`Delta = 0`) with the highest confidence (`Sigma = 70`), the score will be `-4.595068832329223`."},{"metadata":{"trusted":true},"cell_type":"code","source":"def laplace_log_likelihood(y_true, y_pred, sigma):\n    sigma_clipped = np.maximum(sigma, 70)\n    delta_clipped = np.minimum(np.abs(y_true - y_pred), 1000)\n    score = - np.sqrt(2) * delta_clipped / sigma_clipped - np.log(np.sqrt(2) * sigma_clipped)\n    return np.mean(score)\n\nlaplace_log_likelihood(df_train['FVC'], df_train['FVC'], 70)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **3. FVC (Forced Vital Capacity)**\n\n`FVC` measurement shows the amount of air a person can forcefully and quickly exhale after taking a deep breath. It is defined as `the recorded lung capacity in ml` under the Data tab. The change in `FVC` over the course of weeks is used for predicting the patients' lung function decline.\n\nEven though the `FVC` predictions smaller than 1000 are clipped, the minimum value in training set is **827**. The maximum `FVC` value in training set is **6399**. The distribution is heavily tailed on the right end because some patients have extremely high `FVC` measurements. However, most of the patients are close to mean `FVC`."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(f'FVC Statistical Summary\\n{\"-\" * 23}')\n\nprint(f'Mean: {df_train[\"FVC\"].mean():.6}  -  Median: {df_train[\"FVC\"].median():.6}  -  Std: {df_train[\"FVC\"].std():.6}')\nprint(f'Min: {df_train[\"FVC\"].min()}  -  25%: {df_train[\"FVC\"].quantile(0.25)}  -  50%: {df_train[\"FVC\"].quantile(0.5)}  -  75%: {df_train[\"FVC\"].quantile(0.75)}  -  Max: {df_train[\"FVC\"].max()}')\nprint(f'Skew: {df_train[\"FVC\"].skew():.6}  -  Kurtosis: {df_train[\"FVC\"].kurtosis():.6}')\nmissing_values_count = df_train[df_train[\"FVC\"].isnull()].shape[0]\ntraining_samples_count = df_train.shape[0]\nprint(f'Missing Values: {missing_values_count}/{training_samples_count} ({missing_values_count * 100 / training_samples_count:.4}%)')\n\nfig, axes = plt.subplots(ncols=2, figsize=(18, 6), dpi=150)\n\nsns.distplot(df_train['FVC'], label='FVC', ax=axes[0])\nprobplot(df_train['FVC'], plot=axes[1])\n\nfor i in range(2):\n    axes[i].tick_params(axis='x', labelsize=12)\n    axes[i].tick_params(axis='y', labelsize=12)\n    axes[i].set_xlabel('')\n    axes[i].set_ylabel('')\n    \naxes[0].set_title(f'FVC Distribution in Training Set', size=15, pad=15)\naxes[1].set_title(f'FVC Probability Plot', size=15, pad=15)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Every patients' `FVC` should be plotted individually as a function of time because the competition objective is predicting `FVC` values of different patients over multiple timesteps.\n\nMajority of the patients' conditions got worse over the course of weeks except 1-2%. Increase of `FVC` in that small group may be random because `FVC` fluctuates too much over time in some patients. Those patients may not be responding to treatment very well or they are not getting better for some reason. However, some patients are clearly getting better because their `FVC` are increasing linearly with very few fluctuations. Those patients are very rare.\n\nFluctuations don't appear to be predictable and that's why the score is calculated on the last three weeks of every patient and averaged. To avoid potential leakage, every patient's FVC measurement for every possible week must be predicted. That's why `submission.csv` has `Weeks` from **-12** to **133**. Those weeks which are not in the final three visits are ignored in scoring. Since the fluctuations don't appear to be predictable, this can be simplified into a slope prediction problem."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"df_train['FVC_diff-1'] = np.abs(df_train.groupby('Patient')['FVC'].diff(-1))\ndf_train['FVC_First'] = df_train.groupby('Patient')['FVC'].transform('first')\ndf_train['FVC_Last'] = df_train.groupby('Patient')['FVC'].transform('last')\ndf_train['Weeks_First'] = df_train.groupby('Patient')['Weeks'].transform('first')\ndf_train['Weeks_Last'] = df_train.groupby('Patient')['Weeks'].transform('last')\n\ndf_train['2_Data_Point_Slope'] =  (df_train['FVC_Last'] - df_train['FVC_First']) / (df_train['Weeks_Last'] - df_train['Weeks_First'])\n\nfor patient, df in df_train.groupby('Patient'):\n    slope, intercept, _, _, _ = linregress(df['Weeks'], df['FVC'])\n    df_train.loc[df_train['Patient'] == patient, 'All_Data_Point_Slope'] = slope\n    df_train.loc[df_train['Patient'] == patient, 'Intercept'] = intercept\n    \ndef plot_fvc(df, patient):\n    \n    xs = df['FVC_First'].unique().tolist() +  df['FVC_Last'].unique().tolist()\n    ys = df['Weeks_First'].unique().tolist() +  df['Weeks_Last'].unique().tolist()\n            \n    ax = df[(['Weeks', 'FVC',])].set_index('Weeks').plot(figsize=(30, 6), style=['b-'])\n    ax.plot(df['Weeks'], df['All_Data_Point_Slope'] * df['Weeks'] + df['Intercept'], label='Slope (All Data Points)', color='green', linestyle='--')\n    ax.plot(ys, xs, label='Slope (2 Data Points)', color='red', linestyle='--')\n    \n    plt.tick_params(axis='x', labelsize=20)\n    plt.tick_params(axis='y', labelsize=20)\n    plt.xlabel('')\n    plt.ylabel('')\n    plt.title(f'Patient: {patient} - {df[\"Age\"].tolist()[0]} - {df[\"Sex\"].tolist()[0]} - {df[\"SmokingStatus\"].tolist()[0]} ({len(df)} Measurements in {(df[\"Weeks\"].max() - df[\"Weeks\"].min())} Weeks Period)', size=25, pad=25)\n    plt.legend(prop={'size': 18})\n    plt.show()\n\nfor patient, df in list(df_train.groupby('Patient')):\n\n    print(f'Patient: {patient} FVC Statistical Summary\\n{\"-\" * 58}')\n    print(f'Mean: {df[\"FVC\"].mean():.6}  -  Median: {df[\"FVC\"].median():.6}  -  Std: {df[\"FVC\"].std():.6}')\n    print(f'Min: {df[\"FVC\"].min()} -  Max: {df[\"FVC\"].max()}')\n    print(f'Skew: {df[\"FVC\"].skew():.6}  -  Kurtosis: {df[\"FVC\"].kurtosis():.6}\\n{\"-\" * 58}')\n    print(f'Change Mean: {df[\"FVC_diff-1\"].mean():.6}  - Change Median: {df[\"FVC_diff-1\"].median():.6}  - Change Std: {df[\"FVC_diff-1\"].std():.6}')\n    print(f'Change Min: {df[\"FVC_diff-1\"].min()} -  Change Max: {df[\"FVC_diff-1\"].max()}')\n    print(f'Change Skew: {df[\"FVC_diff-1\"].skew():.6} -  Change Kurtosis: {df[\"FVC_diff-1\"].kurtosis():.6}\\n{\"-\" * 58}')\n    print(f'Slope (Calculated with All Data Points): {df[\"All_Data_Point_Slope\"].values[0]:.6}')\n    print(f'Slope (Calculated with 2 Data Points): {df[\"2_Data_Point_Slope\"].values[0]:.6}')\n    \n    plot_fvc(df, patient)\n    \ndf_train.drop(columns=['FVC_diff-1', 'FVC_First', 'FVC_Last', 'Weeks_First', 'Weeks_Last'], inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **4. Tabular Data**"},{"metadata":{},"cell_type":"markdown","source":"### **4.1. Continuous Features**\n\nThere are four continuous features along with `FVC` in tabular data. Those features are:\n\n* `Weeks`: The relative number of weeks pre/post the baseline CT (may be negative). It doesn't have any significant relationship with other features because patients got both better or worse over the course of time regardless of their `Age`.\n* `Percent`: A computed field which approximates the patient's `FVC` as a percent of the typical `FVC` for a person of similar characteristics. This feature has a strong relationship with `FVC` because it is derived from it, but it doesn't have any significant relationship with other features.\n* `Age`: Age of the patient. `Age` has a slight relationship with `FVC` and `Percent` since younger patients have higher lung capacity.\n\nDistributions of `FVC`, `Percent` and `Age` are very similar but `Weeks` is different than those features."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"g = sns.pairplot(df_train[['FVC', 'Weeks', 'Percent', 'Age']], aspect=1.4, height=5, diag_kind='kde', kind='reg')\n\ng.axes[3, 0].set_xlabel('FVC', fontsize=25)\ng.axes[3, 1].set_xlabel('Weeks', fontsize=25)\ng.axes[3, 2].set_xlabel('Percent', fontsize=25)\ng.axes[3, 3].set_xlabel('Age', fontsize=25)\ng.axes[0, 0].set_ylabel('FVC', fontsize=25)\ng.axes[1, 0].set_ylabel('Weeks', fontsize=25)\ng.axes[2, 0].set_ylabel('Percent', fontsize=25)\ng.axes[3, 0].set_ylabel('Age', fontsize=25)\n\ng.axes[3, 0].tick_params(axis='x', labelsize=20)\ng.axes[3, 1].tick_params(axis='x', labelsize=20)\ng.axes[3, 2].tick_params(axis='x', labelsize=20)\ng.axes[3, 3].tick_params(axis='x', labelsize=20)\ng.axes[0, 0].tick_params(axis='y', labelsize=20)\ng.axes[1, 0].tick_params(axis='y', labelsize=20)\ng.axes[2, 0].tick_params(axis='y', labelsize=20)\ng.axes[3, 0].tick_params(axis='y', labelsize=20)\n\ng.fig.suptitle('Tabular Data Feature Distributions and Interactions', fontsize=35, y=1.08)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **4.2. Sex**\n\nThe first categorical feature in tabular data is `Sex` which is basically gender of the patient. In training set, **139** (79%) patients are male and **37** (21%) patients are female."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 6), dpi=100)\n\nsns.barplot(x=df_train.groupby('Patient')['Sex'].first().value_counts().index, y=df_train.groupby('Patient')['Sex'].first().value_counts())\npercentages = [(count / df_train.groupby('Patient')['Sex'].first().value_counts().sum() * 100).round(2) for count in df_train.groupby('Patient')['Sex'].first().value_counts()]\n\nplt.ylabel('')\nplt.xticks(np.arange(2), [f'Male (%{percentages[0]})', f'Female (%{percentages[1]})'])\nplt.tick_params(axis='x', labelsize=12)\nplt.tick_params(axis='y', labelsize=12)\nplt.title('Sex Counts in Training Set', size=15, pad=15)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`FVC` distributions of males and females are very different from each other. Females have lower lung capacity compared to males due to genetics. `FVC` relationships with other features are also very different for males and females. `FVC` of males have a stronger relationship with `Percent` and `Age` compared to `FVC` of females.\n\nComparing `Weeks` for different genders is not logical but females have a decent `FVC` improvement over the course weeks compared to males.\n\n`Percent` distributions of males and females are very different from each other just like `FVC` distributions because `Percent` is derived from it.\n\n`Age` has no differences between males and females in terms of relationships and distributions except female's `Age` distribution have slightly longer tails and a shorter peak."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"g = sns.pairplot(df_train[['FVC', 'Weeks', 'Percent', 'Age', 'Sex']], hue='Sex', aspect=1.4, height=5, diag_kind='kde', kind='reg')\n\ng.axes[3, 0].set_xlabel('FVC', fontsize=25)\ng.axes[3, 1].set_xlabel('Weeks', fontsize=25)\ng.axes[3, 2].set_xlabel('Percent', fontsize=25)\ng.axes[3, 3].set_xlabel('Age', fontsize=25)\ng.axes[0, 0].set_ylabel('FVC', fontsize=25)\ng.axes[1, 0].set_ylabel('Weeks', fontsize=25)\ng.axes[2, 0].set_ylabel('Percent', fontsize=25)\ng.axes[3, 0].set_ylabel('Age', fontsize=25)\n\ng.axes[3, 0].tick_params(axis='x', labelsize=20)\ng.axes[3, 1].tick_params(axis='x', labelsize=20)\ng.axes[3, 2].tick_params(axis='x', labelsize=20)\ng.axes[3, 3].tick_params(axis='x', labelsize=20)\ng.axes[0, 0].tick_params(axis='y', labelsize=20)\ng.axes[1, 0].tick_params(axis='y', labelsize=20)\ng.axes[2, 0].tick_params(axis='y', labelsize=20)\ng.axes[3, 0].tick_params(axis='y', labelsize=20)\n\nplt.legend(prop={'size': 20})\ng._legend.remove()\ng.fig.suptitle('Tabular Data Feature Distributions and Interactions Between Sex Groups', fontsize=35, y=1.08)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **4.3. SmokingStatus**\n\nThe second categorical feature in tabular data is `SmokingStatus` which is also self-explanatory. In training set, **118** (67%) patients are ex-smokers, **49** (28%) patients had never smoked and **9** (5%) patients are smokers."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 6), dpi=100)\n\nsns.barplot(x=df_train.groupby('Patient')['SmokingStatus'].first().value_counts().index, y=df_train.groupby('Patient')['SmokingStatus'].first().value_counts())\npercentages = [(count / df_train.groupby('Patient')['SmokingStatus'].first().value_counts().sum() * 100).round(2) for count in df_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":"`FVC` distributions of `SmokingStatus` groups are quite unexpected. Mean `FVC` of smokers is  higher than mean `FVC` of ex-smokers and patients who had never smoked.\n\nDistribution of `Weeks` is similar for different `SmokingStatus`. Smokers have the strongest positive linear relationship between `FVC` and `Weeks` which is also another unexpected phenomenon.\n\n`Percent` distributions of different `SmokingStatus` groups is very similar to `FVC` distributions but peaks are taller. The linear relationship between `Percent` and `Weeks` is also stronger compared to `FVC` and `Weeks`.\n\n`Age` has no relationship with `SmokingStatus`."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"g = sns.pairplot(df_train[['FVC', 'Weeks', 'Percent', 'Age', 'SmokingStatus']], hue='SmokingStatus', aspect=1.4, height=5, diag_kind='kde', kind='reg')\n\ng.axes[3, 0].set_xlabel('FVC', fontsize=25)\ng.axes[3, 1].set_xlabel('Weeks', fontsize=25)\ng.axes[3, 2].set_xlabel('Percent', fontsize=25)\ng.axes[3, 3].set_xlabel('Age', fontsize=25)\ng.axes[0, 0].set_ylabel('FVC', fontsize=25)\ng.axes[1, 0].set_ylabel('Weeks', fontsize=25)\ng.axes[2, 0].set_ylabel('Percent', fontsize=25)\ng.axes[3, 0].set_ylabel('Age', fontsize=25)\n\ng.axes[3, 0].tick_params(axis='x', labelsize=20)\ng.axes[3, 1].tick_params(axis='x', labelsize=20)\ng.axes[3, 2].tick_params(axis='x', labelsize=20)\ng.axes[3, 3].tick_params(axis='x', labelsize=20)\ng.axes[0, 0].tick_params(axis='y', labelsize=20)\ng.axes[1, 0].tick_params(axis='y', labelsize=20)\ng.axes[2, 0].tick_params(axis='y', labelsize=20)\ng.axes[3, 0].tick_params(axis='y', labelsize=20)\n\nplt.legend(prop={'size': 20})\ng._legend.remove()\ng.fig.suptitle('Tabular Data Feature Distributions and Interactions Between SmokingStatus Groups', fontsize=35, y=1.08)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **4.4. Tabular Feature Correlations**\n\nAs seen from the plots above, the only strong correlation is between `FVC` and `Percent`. The other features' correlations are between -0.1 and 0.1."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 10), dpi=100)\n\nsns.heatmap(df_train[['Weeks', 'FVC', 'Percent', 'Age']].corr(), annot=True, square=True, cmap='coolwarm', annot_kws={'size': 15},  fmt='.2f')   \n\nplt.tick_params(axis='x', labelsize=18, rotation=0)\nplt.tick_params(axis='y', labelsize=18, rotation=0)\nplt.title('Tabular Data Feature Correlations', size=20, pad=20)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **5. DICOM Files**\n\nA `.dcm` file is an image file saved in the Digital Imaging and Communications in Medicine (DICOM) image format. It stores a medical image, such as a CT scan or ultrasound. DCM files may also include patient information to pair the image with the patient. There are **176** directories in `osic-pulmonary-fibrosis-progression/train` and **5** directories in `osic-pulmonary-fibrosis-progression/test`. Every directory has the `.dcm` files of the corresponding patient."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(f'Patients (directories) in osic-pulmonary-fibrosis-progression/train: {len(os.listdir(\"../input/osic-pulmonary-fibrosis-progression/train\"))}')\nprint(f'Patients (directories) in osic-pulmonary-fibrosis-progression/test: {len(os.listdir(\"../input/osic-pulmonary-fibrosis-progression/test\"))}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Multiple `.dcm` files represent different slices of a single CT scan which is acquired at Week 0. CT scans produce 3D volumes for each scan, those volumes consist of 2D slices and each slice is a `.dcm` file. Every directory has different number of slices in `osic-pulmonary-fibrosis-progression/train`. Those number of slices are between **12** and **1018** with a median of **98**. Total number of slices adds up to **33026** for 176 patients."},{"metadata":{"_kg_hide-input":true,"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=(20, 5), dpi=150)\nax = sns.countplot(slice_counts)\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=15)\nplt.tick_params(axis='y', labelsize=15)\nplt.title('Number of Image Slices in Training Set', size=20, pad=20)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"DICOM files can be read and processed easily with `pydicom` package. DICOM files allow to store metadata along with pixel data inside them. The first patient's (`ID00228637202259965313869`) first slice (`1.dcm`) is read with `dcmread` method in the cell below.\n\nReading the file creates a `pydicom.dataset.FileDataset` object. Dataset can be displayed by simply printing its string (`str()` or `repr()`) value. `FileDataset` object wraps `dict` and contains `DataElement` instances. The value of each element can be one of a regular numeric, string or text value, a `list` of regular values or a `Sequence` instance, where `Sequence` is a `list` of `Dataset` instances."},{"metadata":{"trusted":true},"cell_type":"code","source":"file_path = '../input/osic-pulmonary-fibrosis-progression/train/ID00228637202259965313869/1.dcm'\ndicom_file = pydicom.dcmread(file_path)\n\nprint(f'Patient: ID00228637202259965313869 Image: 1.dcm Dataset\\n{\"-\" * 55}\\n\\n{dicom_file}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`pydicom` uses dictionary interface so data inside the files can be accessed with `keys()` and `values()` methods. `.keys()` method shows that keys are tuple pairs and they are called tag numbers. Specific elements can be accessed by their DICOM keywords or tag numbers, but using DICOM keywords is the recommended way. If data is accessed with the tag number, `.value` should be added to end of square brackets."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Accessing Patient Name with DICOM Keyword (PatientName): {dicom_file.PatientName}')\nprint(f'Accessing Patient Name with Tag Number ((0x10, 0x10)): {dicom_file[(0x10, 0x10)].value}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If you don't remember the exact keywords, `Dataset.dir()` can be used to return all non-private element keywords in the dataset. All of the data associated with the keywords below, can be accessed by `Dataset.Keyword`."},{"metadata":{"trusted":true},"cell_type":"code","source":"dicom_file.dir()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`PixelData` keyword contains the raw bytes of the slice. It is much more convenient to use `Dataset.pixel_array` in order to return a `numpy.ndarray`. As seen from below, `ID00228637202259965313869/1.dcm` slice is a 512x512 greyscale image."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(f'Patient: ID00228637202259965313869 Image: 1.dcm Pixel Array {dicom_file.pixel_array.shape}\\n{\"-\" * 70}\\n\\n{dicom_file.pixel_array}')\n\nfig = plt.figure(figsize=(6, 6), dpi=100)\n\nplt.imshow(dicom_file.pixel_array, cmap=plt.cm.bone)\n\nplt.tick_params(axis='x', labelsize=12)\nplt.tick_params(axis='y', labelsize=12)\nplt.title(f'ID00228637202259965313869/1.dcm', size=15, pad=15)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As the images are acquired at Week 0 in one CT scan, a single slice by itself isn't meaningful. All of the slices make up the 3D volume, so they have to be analyzed together. `load_scan` function loads every slice in the given patient directory and stacks them up. However, slice names may not be correctly ordered so it is safer to sort them by `ImagePositionPatient[2]`, that is defined as z coordinates of the upper left hand corner of the slice. The output is a 3D volume with the shape of `(n_slices, n_rows, n_columns)`."},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(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    patient_directory.sort(key=lambda s: float(s.ImagePositionPatient[2]))\n    patient_slices = np.zeros((len(patient_directory), patient_directory[0].Rows, patient_directory[0].Columns))\n\n    for i, s in enumerate(patient_directory):\n        patient_slices[i] = s.pixel_array\n            \n    return patient_slices\n\npatient = 'ID00228637202259965313869'\npatient_slices = load_scan(patient)\nprint(f'Patient {patient} CT scan is loaded - Volume Shape: {patient_slices.shape}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This animation shows the slices in an order. The slices are taken from bottom to top or top to bottom while the patients are holding their breath. That is the reason why area of the lungs is changing between different slices. The lung area is small at bottom and top, but it is larger at middle part."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"%%capture\n\nfig = plt.figure(figsize=(7, 7))\n\nims = []\nfor i in patient_slices:\n    im = plt.imshow(i, animated=True, cmap=plt.cm.bone)\n    plt.axis('off')\n    ims.append([im])\n\nani = animation.ArtistAnimation(fig, ims, interval=25, blit=False, repeat_delay=1000)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"HTML(ani.to_html5_video())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## **6. Image Metadata**"},{"metadata":{},"cell_type":"markdown","source":"### **6.1. Metadata Extraction**\n\nDICOM files consist of lots of different metadata and some of those fields are very useful.\n\nFirst group of useful metadata features are `Rows` (X), `Columns` (Y) and `Slices` (Z). They are the three dimensions that make up the 3D volume of a single CT scan. `Rows` and `Columns` can be accessed from DICOM files but `Slices` doesn't exist in them as a field. It is the count of DICOM files in patient's directories.\n\nSecond group of useful metadata features are `SpacingBetweenSlices`, `PixelSpacing` and `SliceThickness`. `SpacingBetweenSlices` gives the distance between two adjacent slices in mm. This distance can also be calculated with `ImagePositionPatient Z` difference between slices. The other metadata feature is `PixelSpacing` that gives the distance between pixels for both X and Y axis in mm. `SliceThickness` is the nominal slice thickness in mm and it shouldn't be used for volume calculation."},{"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        slice_positions = np.round([s.ImagePositionPatient[2] for s in patient_directory], 4)\n        non_duplicate_idx = np.unique([np.where(slice_position == slice_positions)[0][0] for slice_position in slice_positions])\n    except AttributeError:\n        patient_directory.sort(key=lambda s: int(s.InstanceNumber))\n        instance_numbers = np.array([int(s.InstanceNumber) for s in patient_directory])\n        non_duplicate_idx = np.unique([np.where(instance_number == instance_numbers)[0][0] for instance_number in instance_numbers])\n\n    patient_directory = list(np.array(patient_directory)[non_duplicate_idx])   \n    \n    rows = patient_directory[0].Rows\n    columns = patient_directory[0].Columns\n    slices = len(patient_directory)\n        \n    pixel_spacings = np.zeros((len(patient_directory), 2))\n    slice_positions = np.zeros((len(patient_directory)))\n    slice_thicknesses = []\n        \n    for i, s in enumerate(patient_directory):\n        slice_thicknesses.append(s.SliceThickness)\n        try:\n            pixel_spacings[i, :] = np.array(s.PixelSpacing)\n        except AttributeError:\n            pixel_spacings[i, :] = np.nan\n\n        try:\n            slice_positions[i] = s.ImagePositionPatient[2]\n        except AttributeError:\n            pass\n                \n    df_train.loc[df_train['Patient'] == patient_name, 'Rows'] = rows\n    df_train.loc[df_train['Patient'] == patient_name, 'Columns'] = columns\n    df_train.loc[df_train['Patient'] == patient_name, 'Slices'] = slices\n    df_train.loc[df_train['Patient'] == patient_name, 'PixelSpacingX'] = list(np.round(pixel_spacings.mean(axis=0), 3))[0]\n    df_train.loc[df_train['Patient'] == patient_name, 'PixelSpacingY'] = list(np.round(pixel_spacings.mean(axis=0), 3))[1]\n    df_train.loc[df_train['Patient'] == patient_name, 'SliceSpacing'] = mode(np.abs(np.diff(np.round(slice_positions, 3))))[0][0]\n    df_train.loc[df_train['Patient'] == patient_name, 'SliceThickness'] = mode(slice_thicknesses)[0][0]\n\n        \nfor patient in tqdm.tqdm(df_train['Patient'].unique()):\n    get_metadata(patient)\n    \ndf_train['Rows'] = df_train['Rows'].astype(np.uint16)\ndf_train['Columns'] = df_train['Columns'].astype(np.uint16)\ndf_train['Slices'] = df_train['Slices'].astype(np.uint16)\ndf_train['SliceShape'] = df_train['Rows'].astype(str) + 'x' + df_train['Columns'].astype(str)\ndf_train['SliceThickness'] = df_train['SliceThickness'].astype(np.float32)\ndf_train['SliceSpacing'] = df_train['SliceSpacing'].astype(np.float32)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Two of the patients (`ID00132637202222178761324` and `ID00128637202219474716089`) don't have `SpacingBetweenSlices` and `ImagePositionPatient` fields in their DICOM files, so there is no way to retrieve their `SliceSpacing` feature. However, those missing values can be filled with a value from a similar CT scan. The similarity is mainly based on slice count and other metadata features.\n\nPatient `ID00132637202222178761324` have 407 slices and most of the CT scans in that slice count range have `0.7` mm `SliceSpacing` so it is filled with `0.7`. The other patient is `ID00128637202219474716089` and it is harder to fill because `SliceSpacing` values are very diverse in that slice count range. `5.0` is used for filling the missing `SliceSpacing` in this case because it looks like it is the median value of the range."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Slice Counts between 400 and 410\\n{\"-\" * 32}\\n\\n', df_train.groupby('Patient').first().query('400 < Slices < 410')[['Rows', 'Columns', 'Slices', 'PixelSpacingX', 'PixelSpacingY', 'SliceSpacing', 'SliceThickness']], '\\n')\ndf_train.loc[df_train['Patient'] == 'ID00132637202222178761324', 'SliceSpacing'] = 0.7\nprint(f'Slice Counts between 45 and 55\\n{\"-\" * 30}\\n\\n', df_train.groupby('Patient').first().query('45 < Slices < 55')[['Rows', 'Columns', 'Slices', 'PixelSpacingX', 'PixelSpacingY', 'SliceSpacing', 'SliceThickness']])\ndf_train.loc[df_train['Patient'] == 'ID00128637202219474716089', 'SliceSpacing'] = 5.0","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **6.2. Slice Shapes**\n\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."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12, 9))\nsns.barplot(x=df_train.groupby('Patient')['SliceShape'].first().value_counts().values, y=df_train.groupby('Patient')['SliceShape'].first().value_counts().index)\n\nplt.xlabel('Patients', size=18)\nplt.ylabel('')\nplt.tick_params(axis='x', labelsize=15)\nplt.tick_params(axis='y', labelsize=15)\nplt.title(f'Training Set Slice Shapes', size=18, pad=18)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A random slice from every unique shape is drawn below. It looks like some of the shapes have a frame around the lungs. Those frames exist in only slices in which X is not equal to Y (`1100x888`, `733x888`, `734x888`, `752x888`, `788x888` and `843x888`), and it is related to LightSpeed VCT and LightSpeed Pro 16 CT scanner models. Other shapes (`512x512`, `632x632`, `768x768` and `1302x1302`) don't have those frames."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axes = plt.subplots(figsize=(32, 14), nrows=2, ncols=5)\n\nfor i, patient_name in enumerate(df_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 = patient_directory[0].pixel_array\n    aspect_ratio = df_train.groupby(\"SliceShape\")[\"Patient\"].first().index[i]\n    \n    if i > 4:\n        i -= 5\n        j = 1\n    else:\n        j = 0\n        \n    axes[j][i].imshow(first_slice, cmap=plt.cm.bone)\n\n    axes[j][i].tick_params(axis='x', labelsize=18)\n    axes[j][i].tick_params(axis='y', labelsize=18)\n    axes[j][i].set_title(f'{aspect_ratio}', size=20, pad=20)\n\nplt.show()","execution_count":null,"outputs":[]},{"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":{},"cell_type":"markdown","source":"The 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`. The other shapes aren't changed because they don't any all zero rows or columns. This function can be used safely on the private test set."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(figsize=(32, 14), nrows=2, ncols=5)\n\nfor i, patient_name in enumerate(df_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    old_aspect_ratio = df_train.groupby(\"SliceShape\")[\"Patient\"].first().index[i]\n    new_aspect_ratio = f'{first_slice_cropped.shape[0]}x{first_slice_cropped.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, cmap=plt.cm.bone)\n\n    axes[j][i].tick_params(axis='x', labelsize=18)\n    axes[j][i].tick_params(axis='y', labelsize=18)\n    axes[j][i].set_title(f'{old_aspect_ratio} -> {new_aspect_ratio}', size=20, pad=20)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After the cropping operation, there are only 4 unique slice shapes left and they are `512x512`, `632x632`, `768x768` and `1302x1302`. Those shapes are resized to `512x512` because it is the most common shape in dataset. This process is done with `cv2.resize` and nearest neighbor interpolation is used for preserving the pixel values."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(figsize=(32, 14), nrows=2, ncols=5)\n\nfor i, patient_name in enumerate(df_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 = df_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=18)\n    axes[j][i].tick_params(axis='y', labelsize=18)\n    axes[j][i].set_title(f'{old_aspect_ratio} -> {new_aspect_ratio}', size=20, pad=20)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **6.3. Volume Shapes**\n\n`SliceThickness` is the nominal slice thickness in millimeters. There are 15 unique `SliceThickness` values in training set. There are two unusual values which are `0.80000001` and `0.8999999` because all of the other values are multipliers of `0.125`. Those unusual values are rounded to `0.8` and `0.9`. The most common `SliceThickness` values are `1.0` and `1.25`.\n\n`SliceSpacing` is the distance between slices in millimeters. It has more unique values and less unusual values than `SliceThickness`. It should be used for the volume calculation and it shouldn't be mistaken with `SliceThickness`."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(figsize=(36, 36), nrows=3)\n\nsns.countplot(df_train.groupby('Patient')['Slices'].first(), ax=axes[0])\nsns.barplot(y=df_train.groupby('Patient')['SliceThickness'].first().value_counts().values, x=np.round(df_train.groupby('Patient')['SliceThickness'].first().value_counts().index, 4), ax=axes[1])\nsns.countplot(df_train.groupby('Patient')['SliceSpacing'].first(), ax=axes[2])\n\nfor idx, label in enumerate(axes[0].get_xticklabels()):\n    if idx % 10 == 0:\n        label.set_visible(True)\n    else:\n        label.set_visible(False)\n\nfor i in range(3):\n    axes[i].set_xlabel('')\n    axes[i].set_ylabel('Patient Count', size=20)\n    if i == 0:\n        axes[i].tick_params(axis='x', labelsize=20)\n    else:\n        axes[i].tick_params(axis='x', labelsize=20)\n    \n    axes[i].tick_params(axis='y', labelsize=20)\n    \naxes[0].set_title(f'Training Set Slice Counts', size=25, pad=25)\naxes[1].set_title(f'Training Set Slice Thickness Counts', size=25, pad=25)\naxes[2].set_title(f'Training Set SliceSpacing Counts', size=25, pad=25)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **6.4. Volume Calculation**\n\nAfter the metadata are extracted, volume of the CT scans can be calculated. Features created with segmentation should always be implemented in a normalized unit because shapes of CT scans are different, that's the reason why volume (mm^3) is required.\n\nVolume of a single voxel is calculated with `PixelSpacing X * PixelSpacing Y * SliceSpacing` and voxel count is calculated with `Rows * Columns * Slices`. These two features are enough to calculate volume of selected areas."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train['VoxelVolume'] = df_train['PixelSpacingX'] * df_train['PixelSpacingY'] * df_train['SliceSpacing']\ndf_train['VoxelCount'] = df_train['Rows'].astype(int) * df_train['Columns'].astype(int) * df_train['Slices'].astype(int)\n\ndf_train.drop(columns=['SliceShape', 'PixelSpacingX', 'PixelSpacingY', 'SliceSpacing', 'Rows', 'Columns', 'Slices', 'SliceThickness'], inplace=True)\n\ndf_train.groupby('Patient').first()[['VoxelVolume', 'VoxelCount']].head(10)","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}