{"cells":[{"metadata":{},"cell_type":"markdown","source":"I am trying to make this notebook easy to read, please excuse me for any un-'pythonic' or inefficient code. Any feedback is greatly appreciated!\nFor modelling part of this notebook, visit https://www.kaggle.com/gunawanmarbun/osic-pulmonary-modelling/."},{"metadata":{},"cell_type":"markdown","source":"# Import Packages"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport glob\nimport math\nimport random\nimport functools\nimport warnings\n\nimport scipy\nimport numpy as np\nimport pandas as pd\n\nimport pydicom\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\n\nimport ipywidgets as widgets\nfrom ipywidgets import interact, fixed\nfrom tqdm.auto import tqdm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visualization Utilities\n\n* `distplot_numerical`: distribution plot for numerical columns\n* `distplot_categorical`: distribution plot for categorical columns (max 6 disctinct)\n* `plot_slices_data`: plot DICOM images with gray cmap"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"def distplot_numerical(data, cols_num, col_target=None, grid_c=3, w=15, h_factor=3, **kwargs):\n    \"\"\"\n    Distplot numerical column attributes in small multiple grid.\n\n    Parameters\n    ----------\n    data : pandas.DataFrame\n        dataframe without infinite values. will drop null values while plotting.\n    cols_num : list of str\n        interval or ratio column in data\n    col_target : str, optional\n        the target variable we want to distingusih the cols_num distributino\n    grid_c : int, default=3\n        number of grid columns\n    w : int, default=15\n        figsize witdh arguments\n    h_factor : float, default=3.5\n        height of small plot\n    \"\"\"\n    n = math.ceil(len(cols_num) / grid_c)\n    fig, ax = plt.subplots(n, grid_c, figsize=(w, h_factor*n))\n    if type(ax) != np.ndarray:\n        ax = np.array([ax])\n    sorted_cols_num = sorted(cols_num)  # we wnat it sorted for easier search\n\n    if col_target is None:\n        for col, a in zip(sorted_cols_num, ax.reshape(-1)):\n            sns.distplot(data[col], ax=a, **kwargs)\n            a.set_xlabel(col)\n    else:\n        sorted_cols_target = sorted(data[col_target].unique())\n        if len(sorted_cols_target) > 1 and len(sorted_cols_target) <= 5:  # > 5 will be too crowded\n            for col, a in zip(sorted_cols_num, ax.reshape(-1)):\n                for t in sorted_cols_target:\n                    sns.distplot(data[data[col_target] == t][col].dropna(), ax=a, **kwargs)\n                a.legend(sorted_cols_target)\n                a.set_xlabel(col)\n        else:  # most probably regression analysis\n            for col, a in zip(sorted_cols_num, ax.reshape(-1)):\n                sns.distplot(data[col], ax=a, **kwargs)\n                a.set_xlabel(col)\n    plt.tight_layout()\n        \ndef distplot_categorical(data, cols_cat, col_target=None, normalize=True, grid_c=3, w=15,\n                         h_factor=3, sort=False, kind='bar'):\n    \"\"\"\n    Distplot categorical column attributes in small multiple grid.\n\n    Parameters\n    ----------\n    data : pandas.DataFrame\n        dataframe without infinite values. will drop null values while plotting.\n    cols_cat : list of str\n        categorical column in data\n    col_target : str, optional\n        the target variable we want to distingusih the cols_num distributino\n    normalize : bool, default=True\n        wether to normalize the count or not\n    grid_c : int, default=3\n        number of grid columns\n    w : int, default=15\n        figsize witdh arguments\n    h_factor : float, default=3.5\n        height of small plot\n    sort : bool, default=False\n        prevent sorting based on counts, will fallback to .cat.categories if the series is having\n        category dtype\n    kind : str, default='bar'\n        matplotlib plot kind, really recommend to do bar plot, alternative would be 'barh'\n    \"\"\"\n    n = math.ceil(len(cols_cat) / grid_c)\n    fig, ax = plt.subplots(n, grid_c, figsize=(w, h_factor*n))\n    if type(ax) != np.ndarray:\n        ax = np.array([ax])\n    sorted_cols_cat = sorted(cols_cat)  # we want it sorted for easier search\n\n    if col_target is None:\n        for col, a in zip(sorted_cols_cat, ax.reshape(-1)):\n            data[col].value_counts(normalize=normalize, sort=sort).plot(ax=a, kind=kind)\n            xlabels = [x.get_text()[:15]+'...' if (len(x.get_text()) > 15) else x for x in a.get_xticklabels()]\n            a.set_xticklabels(xlabels, rotation=30, ha='right')\n            a.set_xlabel(col)\n    else:\n        sorted_cols_target = sorted(data[col_target].unique())\n        if len(sorted_cols_target) > 1 and len(sorted_cols_target) <= 6:  # > 5 will be too crowded\n            for col, a in zip(sorted_cols_cat, ax.reshape(-1)):\n                data.groupby(col_target)[col].value_counts(normalize=normalize, sort=sort).unstack(0).plot(ax=a, kind=kind)\n                xlabels = [x.get_text()[:15]+'...' if (len(x.get_text()) > 15) else x for x in a.get_xticklabels()]\n                a.set_xticklabels(xlabels, rotation=30, ha='right')\n        else:  # most probably regression analysis\n            for col, a in zip(sorted_cols_cat, ax.reshape(-1)):\n                data[col].value_counts(normalize=normalize, sort=sort).plot(ax=a, kind=kind)\n                xlabels = [x.get_text()[:15]+'...' if (len(x.get_text()) > 15) else x for x in a.get_xticklabels()]\n                a.set_xticklabels(xlabels, rotation=30, ha='right')\n                a.set_xlabel(col)\n    plt.tight_layout()\n    \n    \ndef plot_slices_data(slices_data, n_cols=10, cmap='gray', **kwargs):\n    n_rows = math.ceil(slices_data.shape[0] / n_cols)\n    fig, axes = plt.subplots(n_rows, n_cols, figsize=(14, n_rows*1.5))\n    for img, ax in tqdm(zip(slices_data, axes.reshape(-1)), leave=False, total=slices_data.shape[0]):\n        ax.imshow(img, cmap=cmap, **kwargs)\n        ax.axis('off')\n    \n    missing_image_cnt = (n_rows * n_cols) - slices_data.shape[0]\n    if missing_image_cnt > 0:\n        for ax in axes.reshape(-1)[::-1][:-missing_image_cnt]:\n            ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Q & A Utils\n\n* `add_first_last_FVC`\n* `add_first_last_weeks`"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def add_first_last_FVC(df_groupby):\n    min_obs, max_obs = df_groupby['Weeks'].agg([\"min\", \"max\"])\n    min_obs_FVC = df_groupby.loc[df_groupby['Weeks'] == min_obs, 'FVC'].values[0]\n    max_obs_FVC = df_groupby.loc[df_groupby['Weeks'] == max_obs, 'FVC'].values[0]\n    is_decline = max_obs_FVC < min_obs_FVC\n    df_groupby['is_decline'] = is_decline\n    df_groupby['min_obs_FVC'] = min_obs_FVC\n    df_groupby['max_obs_FVC'] = max_obs_FVC\n    df_groupby['diff_FVC'] = (max_obs_FVC-min_obs_FVC)\n    df_groupby['diff_pct_FVC'] = (max_obs_FVC-min_obs_FVC) / min_obs_FVC\n    return df_groupby\n\ndef filter_min_max_obs(df_groupby):\n    min_obs, max_obs = df_groupby['Weeks'].agg(['min', 'max'])\n    return df_groupby[df_groupby['Weeks'].isin([min_obs, max_obs])]\n\ndef add_first_last_weeks(df_groupby):\n    min_obs, max_obs = df_groupby['Weeks'].agg([\"min\", \"max\"])\n    df_groupby['min_week'] = min_obs\n    df_groupby['max_week'] = max_obs\n    df_groupby['diff_weeks'] = max_obs - min_obs\n    df_groupby['num_obs'] = df_groupby['Patient'].count()\n    df_groupby['rate_obs'] = df_groupby['num_obs'] / df_groupby['diff_weeks']\n    return df_groupby","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Utilities\n\n* `sort_nicely`: sort filepaths numerically"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import re\n\ndef tryint(s):\n    try:\n        return int(s)\n    except ValueError:\n        return s\n    \ndef alphanum_key(s):\n    \"\"\" Turn a string into a list of string and number chunks.\n        \"z23a\" -> [\"z\", 23, \"a\"]\n    \"\"\"\n    return [ tryint(c) for c in re.split('([0-9]+)', s) ]\n\ndef sort_nicely(l):\n    \"\"\" Sort the given list in the way that humans expect.\n    \"\"\"\n    l.sort(key=alphanum_key)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# References\n\n## Introductory Notebook on DICOM Image by allunia\nI am a fan of her EDA notebook, it includes what I'd say must watch introductory video on how to properly see DICOM images: https://www.kaggle.com/allunia/pulmonary-fibrosis-dicom-preprocessing/notebook\n\n\n## Introductory Notebook on Patient's Data by Heroseo\nAnother recommended read to walk through all patient's data EDA processes: https://www.kaggle.com/piantic/osic-pulmonary-fibrosis-progression-basic-eda\n\n## Watershed Segmentation Algorithm by Ankasor\nhttps://www.kaggle.com/ankasor/improved-lung-segmentation-using-watershed"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"basepath = \"../input/osic-pulmonary-fibrosis-progression/\"\ntrain_df = pd.read_csv(f\"{basepath}train.csv\")\ntest_df = pd.read_csv(f\"{basepath}test.csv\")\nsubmission_df = pd.read_csv(f\"{basepath}sample_submission.csv\")\nprint(train_df.shape, test_df.shape, submission_df.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Patient's Data EDA\n\nThis EDA will use a simple questions & answers format for I believe this is easier to understand. Well, isn't the purpose of EDA is to ask better questions? As this is my first time analyzing medical data, I find it useful to think ahead of all questions ('dumb' or not) and then work my understanding from there. We can skecth all questions beforehand, then work on a simple visualization describing the answer.\n\nLet's first have a quick start from univariate distribution."},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_num = ['Weeks', 'FVC', 'Percent', 'Age']\ncols_cat = ['Sex', 'SmokingStatus']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Univariate Analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"distplot_numerical(train_df, cols_num, grid_c=4)\ndistplot_categorical(train_df, cols_cat, grid_c=2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some evaluation:\n* `Patient`- a unique Id for each patient (also the name of the patient's DICOM folder)\n* `Age`- we see most of patient are in 65-75 range\n* `FVC` - the recorded lung capacity in ml, we see normal distribution centered around 3000 mL, which considering most of our patient is male, is a rather low FVC values\n* `Percent`- a computed field which approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics, we can see from the plot above almost 90% of all patient are below the approximated FVC value\n* `Weeks`- the relative number of weeks pre/post the baseline CT (may be negative). Most of observation is less than a year, we can see later how this is elaborated for each patients.\n* `Sex`- most of the patients are male, which is summarized in the FVC value as we know distribution for male FVC value is around 1000 mL higher than female\n* `SmokingStatus`- could infer some intereseting relationship between smoking behaviour and FVC value development over time."},{"metadata":{},"cell_type":"markdown","source":"## Questions & Answers\n\nExhaustive list of all explorative question."},{"metadata":{},"cell_type":"markdown","source":"### 1. How many patient have increasing/declining FVC measurement over time? What is the rate of change?"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"temp_df = train_df\\\n    .drop_duplicates(subset=['Patient', 'Weeks'], keep='first')\\\n    .groupby('Patient')\\\n    .apply(add_first_last_FVC)\\\n    .drop_duplicates(subset='Patient', keep='first')\\\n    .loc[:, ['Patient', 'is_decline', 'min_obs_FVC', 'max_obs_FVC', 'diff_FVC', 'diff_pct_FVC']]\ntemp_df\\\n    .groupby('is_decline')[['Patient']].agg(['count'])\\\n    .join(temp_df.groupby('is_decline').agg(['mean', 'std']))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2. How many observation we have for each patient? What is the spread of time vs number of observations for each patient?"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"train_df\\\n    .drop_duplicates(subset=['Patient', 'Weeks'], keep='first')\\\n    .groupby('Patient')\\\n    .apply(add_first_last_weeks)\\\n    .drop_duplicates(subset=['Patient'])\\\n    .loc[:, ['diff_weeks', 'num_obs', 'rate_obs']]\\\n    .agg(['min', 'max', 'mean', 'std'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From here we can see that the longest observation we had for a patient is 63 weeks, which is roughly 1.2 year. Later we see in test dataset that we are to predict up to 133 weeks test period, which is more than double the observation period of train dataset. It seems that we could gather domain knowledge of the disease development to patient over the years. From [here](https://www.medscape.com/answers/301226-95979/what-is-the-mortality-rate-of-idiopathic-pulmonary-fibrosis-ipf#:~:text=Idiopathic%20pulmonary%20fibrosis%20(IPF)%20portends,deaths%20per%20million%20in%20women.) for example we know that mean survival of diopathic pulmonary fibrosis (IPF) is 2-5 years from the time of diagnosis. We maybe even can predict the mortality rate over the years."},{"metadata":{},"cell_type":"markdown","source":"# CT-scan EDA\n\nWe create a class `DICOMImages` that wrap around a patient dicom slices."},{"metadata":{"trusted":true},"cell_type":"code","source":"class DICOMImages:\n    DOUBLE_IDS = ['ID00078637202199415319443']\n    \"\"\"Wrapper for multiple slices of a patient CT-Scan results.\"\"\"\n    def __init__(self, id, dirpath='../input/osic-pulmonary-fibrosis-progression/train/'):\n        self.id = id\n        self.basepath = os.path.join(dirpath, self.id)\n        self.filepaths = glob.glob(os.path.join(self.basepath, \"*.dcm\"))\n        if self.id in self.DOUBLE_IDS:\n            self.filepaths = self.filepaths[:len(self.filepaths)//2]\n        sort_nicely(self.filepaths)\n        \n    def __iter__(self):\n        for filepath in self.filepaths:\n            yield pydicom.dcmread(filepath)\n\n    def __len__(self):\n        return len(self.filepaths)\n    \n    @property\n    def image_type(self):\n        \"\"\"\n        Infer dicom image type by its first slice metadata.\n        Categories:\n            - 'zero' : Rescale Intercept value is 0\n            - 'not-zero': Rescale Intercept value is either -1000 or -1024\n        \"\"\"\n        mapper = {0: 'zero'}\n        rescale_intercept = self.get_dicom_metadata(self.get_slice(index=0))['Rescale Intercept']\n        return {\n            'name': mapper.get(rescale_intercept, 'not-zero'),\n            'rescale_intercept': rescale_intercept\n        }\n        \n    @property\n    def slices(self):\n        return list(self.__iter__())\n    \n    def get_slice(self, index):\n        return pydicom.dcmread(self.filepaths[index])\n    \n    @property\n    def df(self):\n        return pd.DataFrame(\n            [self.get_dicom_metadata(slice) for slice in self.__iter__()]\n        )\n    \n    @staticmethod\n    def get_dicom_metadata(slice):\n        dict_res = {}\n        for x in slice.values():\n            if isinstance(x, pydicom.dataelem.RawDataElement):\n                metadata = pydicom.dataelem.DataElement_from_raw(x)\n            else:\n                metadata = x\n            if metadata.name == 'Pixel Data':\n                continue\n            dict_res.update({\n                f\"{metadata.name}\": metadata.value\n            })\n        return dict_res\n    \n    @property\n    def slices_data(self):\n        return np.stack([self._to_HU(slice) for slice in self.__iter__()])\n    \n    @property\n    def middle_slice_data(self):\n        mid_slice_index = (len(self.filepaths)-1) // 2\n        return self._to_HU(pydicom.dcmread(self.filepaths[mid_slice_index]))\n        \n    def sampled_slices_data(self, n_samples=30):\n        if len(self.filepaths) < n_samples:\n            msg = f\"Total slices is less than number of samples: {len(self.filepaths)} < {n_samples}.\"\n            msg += \" Number of samples default to total slices.\"\n            warnings.warn(msg, UserWarning)\n            n_samples = len(self.filepaths)\n        sample_indexes = np.linspace(0, len(self.slices)-1, n_samples).astype(int)\n        sampled_slices = np.array(self.slices)[sample_indexes]\n        return np.stack([self._to_HU(slice) for slice in sampled_slices])\n\n    @staticmethod\n    def _to_HU(slice):\n        intercept, slope = slice.RescaleIntercept, slice.RescaleSlope\n        \n        slice_data = slice.pixel_array.astype(np.int16)\n        slice_data[slice_data <= -1000] = 0\n        \n        if slope != 1:\n            slice_data = slope * slice_data.astype(np.float64)\n            slice_data = slice_data.astype(np.int16)\n            \n        slice_data += np.int16(intercept)\n        return slice_data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Image Metadata EDA\n\nCheck metadata of patient DICOMImage and see if we can found anything interesting. The definition of metadata is taken from https://dicom.innolitics.com/ciods.\n\nWe will specifically see for these image-related columns:\n* `Rows` : total number of pixel in vertical axis\n* `Columns` : total number of pixel in horizontal axis\n* `Pixel Spacing` : The first value is the row spacing in mm, that is the spacing between the centers of adjacent rows, or vertical spacing. The second value is the column spacing in mm, that is the spacing between the centers of adjacent columns, or horizontal spacing.\n* `Bits Allocated` : Number of bits allocated for each pixel sample\n* `High Bit` : Most significant bit for pixel sample data. Each sample shall have the same high bit\n* `Pixel Representation` : Data representation of the pixel samples. Each sample shall have the same pixel representation\n* `Rescale Intercept`: The value b in the relationship between stored values (SV) in Pixel Data (7FE0,0010) and the output units specified in Rescale Type (0028,1054). Output units = m*SV + b.\n* `Rescale Slope` : m in the equation specified by Rescale Intercept"},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_image_related = ['Rows', 'Columns', 'Pixel Spacing',\n                      'Bits Allocated', 'High Bit', 'Pixel Representation',\n                      'Rescale Intercept', 'Rescale Slope']\nall_patient_ids = train_df.Patient.unique()\nall_dicoms = [DICOMImages(id) for id in all_patient_ids]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 1. Does a patient can have different format of DICOM metadata?"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_unique_dict(df):\n    dict_unique = {'Patient': df['Patient ID'].unique()[0]}\n    for col in df.columns:\n        try:\n            dict_unique.update( {f\"{col}\": df[col].nunique()} )\n        except TypeError:\n            dict_unique.update( {f\"{col}\": df[col].astype(str).nunique()} )\n    return dict_unique ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"uniqued_dicom_df = pd.DataFrame([get_unique_dict(dicom.df) for dicom in tqdm(all_dicoms, leave=False)])\nuniqued_dicom_df.to_csv(\"uniqued_dicom_df.csv\", header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print((uniqued_dicom_df[cols_image_related] == 1).sum(axis=0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"uniqued_dicom_df[uniqued_dicom_df['Pixel Spacing'] != 1][['Patient'] + cols_image_related]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We found a patient with 5 unique `Pixel Spacing` values. Let's see in more details for this particular patient."},{"metadata":{"trusted":true},"cell_type":"code","source":"investigate_id = 'ID00099637202206203080121'\ninvestigate_df = DICOMImages(investigate_id).df\ninvestigate_df['Pixel Spacing'] = investigate_df['Pixel Spacing'].astype(str)\ninvestigate_df[cols_image_related].drop_duplicates(subset=['Pixel Spacing'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the largest difference between largest-pixel-spacing slice vs lowest-pixel-spacing slice is around `(0.37, 0.37)` milimeters. We may consider removing this patient from our training model."},{"metadata":{"trusted":true},"cell_type":"code","source":"investigate_slices_data = DICOMImages(investigate_id).slices_data\nplot_slices_data(investigate_slices_data)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* We have much closer lookup view slice on first few images and a rather distant view on the last images.\n* Some slices seem to look 'brighter' then the other slice"},{"metadata":{},"cell_type":"markdown","source":"### 2. How many slices does a patient have?"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2, 1, figsize=(16, 6))\nsns.distplot(uniqued_dicom_df['SOP Instance UID'], kde=False, bins=100, ax=ax[0])\nsns.distplot(uniqued_dicom_df[uniqued_dicom_df['SOP Instance UID'] < 200]['SOP Instance UID'], kde=False, bins=100, ax=ax[1])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most of patients has only around less than 100 slices. When taking slices from patient with a lot of slices we may need to actually downsample the slices."},{"metadata":{},"cell_type":"markdown","source":"### 3. Does all patient have the same dicom format?\n\nWe will now see the distribution for those image related metadata."},{"metadata":{"trusted":true},"cell_type":"code","source":"all_dicoms_df = pd.DataFrame([\n    # We can take the first slice only since the value is unique\n    # for the columns that we'd like to investigate\n    DICOMImages.get_dicom_metadata(dicom.get_slice(0)) for dicom in tqdm(all_dicoms, leave=False)\n])\nprint(all_dicoms_df.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_pixel_spacing_extra = ['Pixel Spacing (row)', 'Pixel Spacing (col)']\nall_dicoms_df[cols_pixel_spacing_extra] = pd.DataFrame(\n    all_dicoms_df['Pixel Spacing'].tolist(), columns=cols_pixel_spacing_extra\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Update our `cols_image_related` to include new columns\ncols_image_related = list(set(cols_image_related + cols_pixel_spacing_extra))\ncols_image_related.remove('Pixel Spacing')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_image_related_num = cols_pixel_spacing_extra\ncols_image_related_cat = [c for c in cols_image_related if c not in cols_image_related_num]\ndistplot_numerical(all_dicoms_df, cols_image_related_num, grid_c=2)\ndistplot_categorical(all_dicoms_df, cols_image_related_cat, grid_c=4)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From here we can take nots:\n* Resizing to either (512, 512) or (768, 768) is reasonable\n* We may need to apply window that is less than 0 (-400 HU if from the paper) before perfoming lung segmentation, depending on the segmentation algorithm\n\nTake two Patent ID with different `Rescale Intercept` to see how our segmentation performs."},{"metadata":{"trusted":true},"cell_type":"code","source":"r_intercept_0_mask = all_dicoms_df['Rescale Intercept'] == 0.0","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4. Are all image square?"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_xy_scatter_sized(df, col_x, col_y):    \n    dff = df.groupby(col_x)[col_y].value_counts()\n    x, y = zip(*dff.index.values)\n    s = dff.to_numpy()\n    \n    plt.figure(figsize=(8, 8))\n    plt.scatter(x, y, s=10*s)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_xy_scatter_sized(all_dicoms_df, 'Rows', 'Columns')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"not_square_mask = all_dicoms_df['Rows'] != all_dicoms_df['Columns']\nnot_square_patient = all_dicoms_df[not_square_mask]['Patient ID'].values[6]\nplot_slices_data(DICOMImages(not_square_patient).sampled_slices_data(24))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After checking all images, a center crop should do find for those images."},{"metadata":{},"cell_type":"markdown","source":"## Image Data EDA\n\nAfter watching the video in the DICOM introductory notebook, I have this simple plotting idea which maybe similar to how radiologist use their software. Lung presets were drawn from this [source](https://arxiv.org/pdf/1811.02651.pdf) which is -700 HU to -600 HU. We will sample a patient CT-scan."},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_id_0_intercept = all_dicoms_df[r_intercept_0_mask]['Patient ID'].values[0]\nsample_id_not0_intercept = all_dicoms_df[~r_intercept_0_mask]['Patient ID'].values[0]\n\nprint(f\"Patient ID (0 Intercept): {sample_id_0_intercept}\")\nprint(f\"Patient ID (Not-0 Intercept): {sample_id_not0_intercept}\")\n\nsampled_slices_data_0_intercept = DICOMImages(sample_id_0_intercept).sampled_slices_data(n_samples=30)\nsampled_slices_data_not0_intercept = DICOMImages(sample_id_not0_intercept).sampled_slices_data(n_samples=30)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Lung Segmentation\n\nFor each of segmentation method we will see:\n* For 0 rescale intercept slices\n* For not 0 rescale intercept slices\n* For 0 rescale intercept slices (windowed/thresholded)\n* For not 0 rescale intercept slices (windowed/thresholded)\n\nFor every slices we will only sample 30 slices with equal intervals.\n\nLet's first try simple image plot with thresholded image to between range of -1000 HU up to -400 HU."},{"metadata":{"trusted":true},"cell_type":"code","source":"import cv2\nfrom skimage import measure, morphology, segmentation\n\n\ndef threshold_slices_data(slices_data, low=-1000, high=-400):\n    copy = slices_data.copy()\n    copy[copy < low] = low\n    copy[copy > high] = high\n    return copy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_slices_data(sampled_slices_data_0_intercept)\nplot_slices_data(sampled_slices_data_not0_intercept)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_slices_data(threshold_slices_data(sampled_slices_data_0_intercept, low=-1000, high=-400))\nplot_slices_data(threshold_slices_data(sampled_slices_data_not0_intercept, low=-1000, high=-400))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Watershed Lung Segmentation"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def watershed_separate_lungs(image, threshold_low=-1000, output_shape=(512, 512), **kwargs):\n    \"\"\"\n    Segment lung image using watershed algorithm\n    \n    Parameters\n    ----------\n    image : numpy.ndarray\n        hounsfield units (HU) transformed image\n    threshold_low : int, default=-1000\n        lower HU threshold for image\n    output_shape : tuple, default=(512, 512)\n        desired output masked segmented lung image\n    kwargs\n        kwarg for generate_markers()\n        \n    Returns\n    -------\n    segmented : numpy.ndarray(shape=output_shape)\n        segmented lung image\n    \"\"\"\n    #Creation of the markers as shown above:\n    marker_internal, marker_external, marker_watershed = generate_markers(\n        image,\n        output_shape=output_shape,\n        **kwargs\n    )\n    \n    #Creation of the Sobel-Gradient\n    sobel_filtered_dx = scipy.ndimage.sobel(image, 1)\n    sobel_filtered_dy = scipy.ndimage.sobel(image, 0)\n    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n    sobel_gradient *= 255.0 / np.max(sobel_gradient)\n    \n    #Watershed algorithm\n    watershed = morphology.watershed(sobel_gradient, marker_watershed)\n    \n    #Reducing the image created by the Watershed algorithm to its outline\n    outline = scipy.ndimage.morphological_gradient(watershed, size=(3,3))\n    outline = outline.astype(bool)\n    \n    #Performing Black-Tophat Morphology for reinclusion\n    #Creation of the disk-kernel and increasing its size a bit\n    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [1, 1, 1, 1, 1, 1, 1],\n                       [0, 1, 1, 1, 1, 1, 0],\n                       [0, 0, 1, 1, 1, 0, 0]]\n    blackhat_struct = scipy.ndimage.iterate_structure(blackhat_struct, 8)\n    #Perform the Black-Hat\n    outline += scipy.ndimage.black_tophat(outline, structure=blackhat_struct)\n    \n    #Use the internal marker and the Outline that was just created to generate the lungfilter\n    lungfilter = np.bitwise_or(marker_internal, outline)\n    #Close holes in the lungfilter\n    #fill_holes is not used here, since in some slices the heart would be reincluded by accident\n    lungfilter = scipy.ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n    \n    #Apply the lungfilter (note the filtered areas being assigned to specified threshold HU)\n    segmented = np.where(lungfilter == 1,\n                         image,\n                         threshold_low*np.ones(output_shape))\n    \n    return segmented\n\ndef generate_markers(image, threshold=-600, output_shape=(512, 512)):\n    \"\"\"\n    Create watershed marker matrix.\n    \n    Parameters\n    ----------\n    image : numpy.ndarray\n        hounsfield units (HU) transformed image\n    threshold : int, default=-600\n        threshold of internal marker, defaulting -600 for lung segmentation\n    output_shape : tuple, default=(512, 512)\n        desired output shape of marker_watershed\n        \n    Returns\n    -------\n    marker_internal : numpy.ndarray\n    marker_external : numpy.ndarray\n    watershed_marker : numpy.ndarray\n    \"\"\"\n    marker_internal = image < threshold\n    marker_internal = segmentation.clear_border(marker_internal)\n    marker_internal_labels = measure.label(marker_internal)\n    areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n    areas.sort()\n    if len(areas) > 2:\n        for region in measure.regionprops(marker_internal_labels):\n            if region.area < areas[-2]:\n                for coordinates in region.coords:                \n                       marker_internal_labels[coordinates[0], coordinates[1]] = 0\n    marker_internal = marker_internal_labels > 0\n    \n    external_a = scipy.ndimage.binary_dilation(marker_internal, iterations=10)\n    external_b = scipy.ndimage.binary_dilation(marker_internal, iterations=55)\n    marker_external = external_b ^ external_a\n    \n    marker_watershed = np.zeros(output_shape, dtype=np.int)\n    marker_watershed += marker_internal * 255\n    marker_watershed += marker_external * 128\n    \n    return marker_internal, marker_external, marker_watershed\n\n\ndef plot_watershed_segmentation(slices_data, cmap='Blues_r'):\n    cnt = slices_data.shape[0]\n    rows = cnt // 10\n    fig, axes = plt.subplots(rows, 10, figsize=(14, rows*2))\n    pbar = tqdm(zip(slices_data, axes.reshape(-1)), leave=False, total=cnt)\n    for img, ax in pbar:\n        segmented_img = watershed_separate_lungs(img,\n                                                 threshold_low=-2000,\n                                                 output_shape=(512, 512),\n                                                 threshold=-400)\n        ax.imshow(segmented_img, cmap=cmap)\n        ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nplot_watershed_segmentation(sampled_slices_data_0_intercept)\nplot_watershed_segmentation(sampled_slices_data_not0_intercept)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Morphological Closing Segmentation"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"from skimage.filters import threshold_otsu, median\nfrom skimage.segmentation import clear_border\nfrom skimage import morphology\nfrom scipy.ndimage import binary_fill_holes\n\n\ndef lung_segment(img, display=False):\n    thresh = threshold_otsu(img)\n    binary = img <= thresh\n\n    lungs = median(clear_border(binary))\n    lungs = morphology.binary_closing(lungs, selem=morphology.disk(7))\n    lungs = binary_fill_holes(lungs)\n\n    final = lungs*img\n    final[final == 0] = np.min(img)\n    \n    if display:\n        fig, ax = plt.subplots(1, 4, figsize=(15, 15))\n\n        ax[0].set_title('HU Image')\n        ax[0].imshow(img, cmap='gray')\n        ax[0].axis('off')\n\n        ax[1].set_title('Thresholded Image')\n        ax[1].imshow(binary, cmap='gray')\n        ax[1].axis('off')\n\n        ax[2].set_title('Lungs Mask')\n        ax[2].imshow(lungs, cmap='gray')\n        ax[2].axis('off')\n\n        ax[3].set_title('Final Image')\n        ax[3].imshow(final, cmap='gray')\n        ax[3].axis('off')\n    \n    return final, lungs\n\n\ndef plot_morphological_closing_segmentation(slices_data):\n    cnt = slices_data.shape[0]\n    rows = cnt // 10\n    fig, axes = plt.subplots(rows, 10, figsize=(14, rows*2))\n    pbar = tqdm(zip(slices_data, axes.reshape(-1)), leave=False, total=cnt)\n    for img, ax in pbar:\n        segmented_img, segmented_lungs = lung_segment(img)\n        ax.imshow(segmented_img, cmap='Blues_r')\n        ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nplot_morphological_closing_segmentation(sampled_slices_data_0_intercept)\nplot_morphological_closing_segmentation(sampled_slices_data_not0_intercept)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets try thresholding the image first then performing the morphological segmentation."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nplot_morphological_closing_segmentation(threshold_slices_data(sampled_slices_data_0_intercept))\nplot_morphological_closing_segmentation(threshold_slices_data(sampled_slices_data_not0_intercept))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Allunia Final Segmentation"},{"metadata":{"trusted":true},"cell_type":"code","source":"def allunia_final_segment(slice, hu_max=-320):\n    binary_image = np.array(slice > hu_max, dtype=np.int8)+1\n    labels = measure.label(binary_image)\n\n    background_label_1 = labels[0,0]\n    background_label_2 = labels[0,-1]\n    background_label_3 = labels[-1,0]\n    background_label_4 = labels[-1,-1]\n\n    #Fill the air around the person\n    binary_image[background_label_1 == labels] = 2\n    binary_image[background_label_2 == labels] = 2\n    binary_image[background_label_3 == labels] = 2\n    binary_image[background_label_4 == labels] = 2\n\n    #We have a lot of remaining small signals outside of the lungs that need to be removed. \n    #In our competition closing is superior to fill_lungs \n    selem = morphology.disk(4)\n    binary_image = morphology.closing(binary_image, selem)\n\n    binary_image -= 1 #Make the image actual binary\n    binary_image = 1-binary_image # Invert it, lungs are now 1\n    return binary_image\n\n\ndef plot_allunia_segmentation(slices_data, cmap='Blues_r'):\n    cnt = slices_data.shape[0]\n    rows = cnt // 10\n    fig, axes = plt.subplots(rows, 10, figsize=(14, rows*2))\n    pbar = tqdm(zip(slices_data, axes.reshape(-1)), leave=False, total=cnt)\n    for img, ax in pbar:\n        segmented_img = allunia_final_segment(img)\n        ax.imshow(segmented_img, cmap=cmap)\n        ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nplot_allunia_segmentation(sampled_slices_data_0_intercept)\nplot_allunia_segmentation(sampled_slices_data_not0_intercept)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Howard Chen Segmentation (raddq.com)"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.cluster import KMeans\n\n\ndef raddq_segmentation(img, display=False):\n    \"\"\"\n    What changes:\n    * Last dilation from (10, 10) to (15, 15)\n    \"\"\"\n    row_size = img.shape[0]\n    col_size = img.shape[1]\n    \n    mean = np.mean(img)\n    std = np.std(img)\n    img = img-mean\n    img = img/std\n    # Find the average pixel value near the lungs\n    # to renormalize washed out images\n    middle = img[int(col_size/5):int(col_size/5*4),int(row_size/5):int(row_size/5*4)] \n    mean = np.mean(middle)  \n    max = np.max(img)\n    min = np.min(img)\n    # To improve threshold finding, I'm moving the \n    # underflow and overflow on the pixel spectrum\n    img[img==max]=mean\n    img[img==min]=mean\n    #\n    # Using Kmeans to separate foreground (soft tissue / bone) and background (lung/air)\n    #\n    kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))\n    centers = sorted(kmeans.cluster_centers_.flatten())\n    threshold = np.mean(centers)\n    thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image\n\n    # First erode away the finer elements, then dilate to include some of the pixels surrounding the lung.  \n    # We don't want to accidentally clip the lung.\n\n    eroded = morphology.erosion(thresh_img,np.ones([3,3]))\n    dilation = morphology.dilation(eroded,np.ones([8,8]))\n\n    labels = measure.label(dilation) # Different labels are displayed in different colors\n    label_vals = np.unique(labels)\n    regions = measure.regionprops(labels)\n    good_labels = []\n    for prop in regions:\n        B = prop.bbox\n        if B[2]-B[0]<row_size/10*9 and B[3]-B[1]<col_size/10*9 and B[0]>row_size/5 and B[2]<col_size/5*4:\n            good_labels.append(prop.label)\n    mask = np.ndarray([row_size,col_size],dtype=np.int8)\n    mask[:] = 0\n\n    #  After just the lungs are left, we do another large dilation\n    #  in order to fill in and out the lung mask \n    for N in good_labels:\n        mask = mask + np.where(labels==N,1,0)\n    mask = morphology.dilation(mask, np.ones([15, 15])) # 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\n\n\ndef plot_raddq_segmentation(slices_data, cmap='Blues_r'):\n    cnt = slices_data.shape[0]\n    rows = cnt // 10\n    fig, axes = plt.subplots(rows, 10, figsize=(14, rows*2))\n    pbar = tqdm(zip(slices_data, axes.reshape(-1)), leave=False, total=cnt)\n    for img, ax in pbar:\n        segmented_img = raddq_segmentation(img)\n        ax.imshow(segmented_img, cmap=cmap)\n        ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nplot_raddq_segmentation(sampled_slices_data_0_intercept)\nplot_raddq_segmentation(sampled_slices_data_not0_intercept)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nplot_raddq_segmentation(threshold_slices_data(sampled_slices_data_0_intercept))\nplot_raddq_segmentation(threshold_slices_data(sampled_slices_data_not0_intercept))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Healthy vs Non Healthy Scan\n\nHaving an understanding on how important is windowing in DICOM images, we will try to find how does FVC relate to a patient CT-scan results. The first question that came to my mind was can we 'naively' see how the progression of the lung images correlate to the FVC? As we know that a patient only take CT Scan once, what we can do is take images from high FVC & low FVC patient and compare the slices."},{"metadata":{"trusted":true},"cell_type":"code","source":"worst_patient, best_patient = train_df\\\n    .loc[train_df['Weeks'] == 0, ['Patient', 'FVC']]\\\n    .sort_values(by='FVC')\\\n    .iloc[[0, -1]].values\nworst_patient, best_patient","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"worst_patient_slices_data = DICOMImages(worst_patient[0]).sampled_slices_data(30)\nbest_patient_slices_data = DICOMImages(best_patient[0]).sampled_slices_data(30)\nplot_slices_data(threshold_slices_data(worst_patient_slices_data, low=-1000, high=-400))\nplot_slices_data(threshold_slices_data(best_patient_slices_data, low=-1000, high=-400))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sadly, not much can be inferred from these image comparison as we lack the domain knowledge needed for the analysis. Let's move on and try just a simple one middle image modelling!"},{"metadata":{},"cell_type":"markdown","source":"# Dataset Creation\n\nPreprocessing steps:\n* If not square, crop to lowest dimension\n* Threshold image\n* Morphological Segmentation\n* Manual unique another center crop\n* Center crop image 5% from all sides\n* Recenter image using center of maxx\n\nAs the dataset is quite small, we will exhaustively see all patient image preprocessed."},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage.transform import resize\nfrom matplotlib.patches import Circle\n\n\ndef crop_to_square(image):\n    width, height = image.shape\n    if width != height:\n        min_ = min(width, height)\n        if min_ == width:\n            top = (height - width) // 2\n            squared = image[ :, top:top+width].copy()\n        else:\n            left = (width - height) // 2\n            squared = image[ left:left+height, :].copy()\n    else:\n        squared = image\n    assert squared.shape[0] == squared.shape[1]\n    return squared\n\n\ndef morphological_segmentation(img):\n    segmented_img, _ = lung_segment(img)\n    return segmented_img\n\n\ndef center_crop_pad(image, pct=0.08):\n    \"\"\"\n    Parameters\n    ----------\n    image : numpy.ndarray\n        slice of image\n    pct : float or int, default=0.05\n        if float, crop pct% of from all image sides, output \n            shape of (1-pct)% of original image\n        if int, crop pct pixels from all image sides\n    \"\"\"\n    original_width, original_height = image.shape\n    if type(pct) == float:\n        left = right = int(pct * original_width)\n        top = bottom = int(pct * original_height)\n    elif type(pct) == int:\n        left = right = pct\n        top = bottom = pct\n\n    cropped_image = image[ left:original_width-right, top:original_height-bottom ].copy()\n    padded_image = np.pad(cropped_image, [(left, right), (top, bottom)], mode='minimum')\n    return padded_image\n\n\ndef recenter_image(slice_data):\n    copy_ = slice_data.copy()\n    width, height = copy_.shape\n    min_val = copy_.min()\n    copy_[copy_ != min_val] = 1\n    copy_[copy_ == min_val] = 0\n    cx, cy = scipy.ndimage.measurements.center_of_mass(copy_)\n    return cx, cy\n\ndef plot_preprocess_steps(slice_data):\n    squared = crop_to_square(slice_data)\n    segmented = morphological_segmentation(squared)\n    crop_pad = center_crop_pad(segmented)\n    resized = resize(crop_pad, output_shape=(512, 512))\n    cx, cy = recenter_image(resized)\n    \n    fig, axes = plt.subplots(3, 2, figsize=(14, 18))\n    axes = axes.reshape(-1)\n    axes[0].imshow(slice_data, cmap='Blues_r')\n    axes[0].set_title(\"Original Image\")\n    axes[1].imshow(squared, cmap='Blues_r')\n    axes[1].set_title(\"Squared Image\")\n    axes[2].imshow(segmented, cmap='Blues_r')\n    axes[2].set_title(\"Segmented Image\")\n    axes[3].imshow(crop_pad, cmap='Blues_r')\n    axes[3].set_title(\"Center Crop & Padded Image\")\n    axes[4].imshow(resized, cmap='Blues_r')\n    axes[4].add_patch(Circle((cx, cy), 24, color='y'))\n    axes[4].set_title(\"Resized Image\")\n    axes[5].imshow(resized, cmap='Blues_r')\n    axes[5].set_title(\"Centered Image\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_id_0_intercept = DICOMImages(all_dicoms_df[r_intercept_0_mask]['Patient ID'].sample(1).values[0])\nsample_id_not0_intercept = DICOMImages(all_dicoms_df[~r_intercept_0_mask]['Patient ID'].sample(1).values[0])\nis_not_square = all_dicoms_df['Rows'] != all_dicoms_df['Columns']\nsample_id_not_square = DICOMImages(all_dicoms_df[is_not_square]['Patient ID'].sample(1).values[0])\n\nprint(f\"Patient ID (0 Intercept): {sample_id_0_intercept.id}\")\nprint(f\"Patient ID (Not-0 Intercept): {sample_id_not0_intercept.id}\")\nprint(f\"Patient ID (Not Square): {sample_id_not_square.id}\")\n\nsampled_slices_data_0_intercept = sample_id_0_intercept.sampled_slices_data(n_samples=30)\nsampled_slices_data_not0_intercept = sample_id_not0_intercept.sampled_slices_data(n_samples=30)\nsampled_slices_data_not_square = sample_id_not_square.sampled_slices_data(n_samples=30)\n\nmiddle_slice_data_0_intercept = sample_id_0_intercept.middle_slice_data\nmiddle_slice_data_not0_intercept = sample_id_not0_intercept.middle_slice_data\nmiddle_slice_data_not_square = sample_id_not_square.middle_slice_data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"binaried_image = middle_slice_data_0_intercept.copy()\nbinaried_image[binaried_image == 0] = -1000\nplot_preprocess_steps(threshold_slices_data(binaried_image, low=-1000, high=-400))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_preprocess_steps(threshold_slices_data(middle_slice_data_not0_intercept, low=-1000, high=-400))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"binaried_image = middle_slice_data_not_square.copy()\nbinaried_image[binaried_image == 0] = -1000\nplot_preprocess_steps(threshold_slices_data(binaried_image, low=-1000, high=-400))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Semi-Automatic Bounding Box Alignment\n\nThe idea is to manually draw the bounding box to remove 'patient bed', 'small blobs', etc. after segmentation. Preset bbox will be drawn based on not -1000 value pixel. We will tune for all the patient dicoms.\n\n[**NOT IMPLEMENTED**]\n* Rotate the image based on patient orientation during CT-scan (or just perform rotation augmentation)\n* Contrast adjustment (or contrast augmentation)\n* Center image (or shift augmentation"},{"metadata":{"trusted":true},"cell_type":"code","source":"import IPython.display as ipd\nfrom matplotlib.patches import Rectangle\n\n\ndef segment_lung(slice_data, image_type, segment_func):\n    if image_type == 'zero':\n        slice_data[slice_data == 0] = -1000\n    segmented_image = segment_func(threshold_slices_data(slice_data, low=-1000, high=-400))\n    return segmented_image\n\n\ndef infer_bounding_box(segmented_image):\n    y_match, x_match = np.where(segmented_image != -1000)\n    y_min, x_min = y_match.min(), x_match.min()\n    y_max, x_max = y_match.max(), x_match.max()\n    width = abs(x_max - x_min)\n    height = abs(y_max - y_min)\n    return BoundingBox((x_min, y_min), width, height)\n\n\ndef draw_with_bounding_box(segmented_image, bbox):\n    fig, axes = plt.subplots(figsize=(14, 6))\n    axes.imshow(segmented_image, cmap='Blues_r')\n    bbox_patch = Rectangle(*bbox.attribute_list,\n                           fill=False,\n                           color='yellow')\n    axes.add_patch(bbox_patch)\n    plt.show()\n    \n    \nclass BoundingBox:\n    \"\"\"Initiation of bbox follows matplotlib Rectangle patch\"\"\"\n    def __init__(self, xy, width, height):\n        self.x, self.y = xy\n        self.width = width\n        self.height = height\n        \n    @property\n    def attribute_list(self):\n        return [(self.x, self.y), self.width, self.height]\n    \n    def __repr__(self):\n        return f\"Bbox (bottom left width height): {self.x} {self.y} {self.width} {self.height}\"\n\n    \ndef crop_recenter(image, bbox, pad_value=-1000):\n    x, y, width, height = bbox.x, bbox.y, bbox.width, bbox.height\n    cropped_image = image[ y:y+height, x:x+width ]\n    out_height, out_width = image.shape\n    \n    padded_image = np.ones(image.shape, dtype=np.int16) * pad_value\n    x_start = (out_width - width) // 2\n    y_start = (out_height - height) // 2\n    padded_image[ y_start:y_start+height, x_start:x_start+width ] = cropped_image\n    return padded_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from ipywidgets import IntSlider, interact, fixed, interact_manual, Text\n\n\ndef manual_bbox(slice_data, image_type, x=None, y=None, x_max=None, y_max=None):\n    width = x_max - x\n    height = y_max - y\n    segment_func = morphological_segmentation\n    segmented_image = segment_lung(slice_data, image_type, segment_func=segment_func)\n    bbox = infer_bounding_box(segmented_image)\n    print(\"Inferred bbox (bottom-left), width, height:\", bbox.x, bbox.y, bbox.width, bbox.height)\n    bbox.x = x or bbox.x\n    bbox.y = y or bbox.y\n    bbox.width = width or bbox.width\n    bbox.height = height or bbox.height\n\n    fig, axes = plt.subplots(1, 3, figsize=(14, 6))\n    axes[0].imshow(slice_data, cmap='Blues_r')\n    axes[0].set_title('Original Image')\n\n    # Segmented image with bbox\n    bbox_patch = Rectangle(*bbox.attribute_list,\n                           fill=False,\n                           color='yellow')\n    axes[1].imshow(segmented_image, cmap='Blues_r')\n    axes[1].add_patch(bbox_patch)\n    axes[1].set_title(f'Segmented ({segment_func.__name__})')\n\n    # Crop centered image according to bbox\n    crop_recentered_image = crop_recenter(segmented_image, bbox)\n    axes[2].imshow(crop_recentered_image, cmap='Blues_r')\n    axes[2].set_title('Bbox Crop & Centered')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Filter out those ids we won't use in training."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define all the 'unique' ids\nvarying_pixel_spacing_ids = ['ID00099637202206203080121']\nfailed_segmentation_ids = ['ID00026637202179561894768']\nbad_ids = ['ID00011637202177653955184', 'ID00052637202186188008618']\n\n# Define the filter mask\nnot_bad_ids = ~all_dicoms_df['Patient ID'].isin(bad_ids)\nzero_intercept = all_dicoms_df['Rescale Intercept'] == 0\n\n# Apply the filter mask\nall_zero_intercept_ids = all_dicoms_df[not_bad_ids & zero_intercept]['Patient ID']\nall_zero_dicoms = [DICOMImages(id) for id in all_zero_intercept_ids]\nall_not_zero_intercept_ids = all_dicoms_df[not_bad_ids & ~zero_intercept]['Patient ID']\nall_not_zero_dicoms = [DICOMImages(id) for id in all_not_zero_intercept_ids]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Zero Intercept Dicoms\n\nFailed segmentation is coded (0, 0, 1, 1)\n\nThis is actually interactive, you just need to edit the `is_final` variable. For commit run we only show the first patient middle slice."},{"metadata":{"trusted":true},"cell_type":"code","source":"threshold_map_zero_intercept = []\nmore_broken_ids = []\nmore_exceptions = []\nis_final = 'q'  # change to 'n' to enable interactive mode\ni = 1\ntotal = len(all_zero_dicoms)\nsegment_func = morphological_segmentation\nfor dicoms in all_zero_dicoms:\n    try:\n        ipd.clear_output(wait=True)\n        middle_slice_data = dicoms.middle_slice_data\n        image_type = dicoms.image_type['name']\n        segmented_image = segment_lung(middle_slice_data, image_type, segment_func=segment_func)\n        bbox = infer_bounding_box(segmented_image)\n        manual_bbox(middle_slice_data, image_type, bbox.x, bbox.y, bbox.x + bbox.width, bbox.y + bbox.height)\n\n        while is_final not in ['y', 'q']:\n            print(f\"{i}/{total}\")\n            is_final = input()\n            if is_final in ['y', 'q']:\n                break\n\n            bbox.x, bbox.y, bbox.width, bbox.height = map(int, input().split())\n            print(bbox)\n    #         print(f\"New bbox (bottom-left), width, height: {bbox.attribute_list}\")\n            manual_bbox(middle_slice_data, image_type, bbox.x, bbox.y, bbox.x + bbox.width, bbox.y + bbox.height)\n\n        if is_final == 'q':\n            break\n\n        threshold_map_zero_intercept.append({\n            \"patient\": dicoms.id,\n            \"x\": bbox.x,\n            \"y\": bbox.y,\n            \"width\": bbox.width,\n            \"height\": bbox.height\n        })\n        i += 1\n        is_final = 'n'\n    except Exception as e:\n        more_exceptions.append(e)\n        more_broken_ids.append(dicoms.id)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Don't forget to save all the hardwork! :)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# df_map_zero = pd.DataFrame(threshold_map_zero_intercept)\n# df_map_zero.to_csv('threshold_map_zero.csv', header=True, index=False)\n\n# df_map_zero_exceptions = pd.DataFrame(list(zip(more_broken_ids, more_exceptions)), columns=['id', 'exception'])\n# df_map_zero_exceptions.to_csv('threshold_map_zero_exceptions.csv', header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Not Zero Intercept Dicoms\n\nFailed segmentation is coded (0, 0, 1, 1)\n\nAgain, this is actually interactive, you just need to edit the `is_final` variable. For commit run we only show the first patient middle slice."},{"metadata":{"trusted":true},"cell_type":"code","source":"threshold_map_not_zero_intercept = []\nmore_broken_ids = []\nmore_exceptions = []\nis_final = 'q'  # # change to 'n' to enable interactive mode\ni = 1\ntotal = len(all_not_zero_dicoms)\nsegment_func = morphological_segmentation\nfor dicoms in all_not_zero_dicoms:\n    try:\n        ipd.clear_output(wait=True)\n        middle_slice_data = dicoms.middle_slice_data\n        image_type = dicoms.image_type['name']\n        segmented_image = segment_lung(middle_slice_data, image_type, segment_func=segment_func)\n        bbox = infer_bounding_box(segmented_image)\n        manual_bbox(middle_slice_data, image_type, bbox.x, bbox.y, bbox.x + bbox.width, bbox.y + bbox.height)\n\n        while is_final not in ['y', 'q']:\n            print(f\"{i}/{total}\")\n            is_final = input()\n            if is_final in ['y', 'q']:\n                break\n\n            bbox.x, bbox.y, bbox.width, bbox.height = map(int, input().split())\n            print(bbox)\n    #         print(f\"New bbox (bottom-left), width, height: {bbox.attribute_list}\")\n            manual_bbox(middle_slice_data, image_type, bbox.x, bbox.y, bbox.x + bbox.width, bbox.y + bbox.height)\n\n        if is_final == 'q':\n            break\n\n        threshold_map_not_zero_intercept.append({\n            \"patient\": dicoms.id,\n            \"x\": bbox.x,\n            \"y\": bbox.y,\n            \"width\": bbox.width,\n            \"height\": bbox.height\n        })\n        i += 1\n        is_final = 'n'\n    except Exception as e:\n        more_exceptions.append(e)\n        more_broken_ids.append(dicoms.id)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Don't forget to save all the hardwork! :)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# df_map_not_zero = pd.DataFrame(threshold_map_not_zero_intercept)\n# df_map_not_zero.to_csv('threshold_map_not_zero.csv', header=True, index=False)\n\n# df_map_not_zero_exceptions = pd.DataFrame(list(zip(more_broken_ids, more_exceptions)), columns=['id', 'exception'])\n# df_map_not_zero_exceptions.to_csv('threshold_map_not_zero_exceptions.csv', header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Dumpster\n\n\nCode that I don't want to throw away. This code lets you explore different threshold for an image slice by using slider."},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"Interactive Part\"\"\"\nsample_dicoms = DICOMImages('ID00376637202297677828573')\nmiddle_slice_data = sample_dicoms.middle_slice_data\nimage_type = sample_dicoms.image_type['name']\nsegmented_image = segment_lung(middle_slice_data, image_type, segment_func=segment_func)\nbbox = infer_bounding_box(segmented_image)\nx_ = IntSlider(min=0, max=middle_slice_data.shape[1], value=bbox.x)\ny_ = IntSlider(min=0, max=middle_slice_data.shape[1], value=bbox.y)\nx_max_ = IntSlider(min=0, max=middle_slice_data.shape[0], value=bbox.x + bbox.width)\ny_max_ = IntSlider(min=0, max=middle_slice_data.shape[0], value=bbox.y + bbox.height)\nvalues = interact_manual(\n    manual_bbox,\n    slice_data=fixed(middle_slice_data),\n    image_type=fixed(image_type),\n    x=x_,\n    y=y_,\n    x_max=x_max_,\n    y_max=y_max_,\n)","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}