{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Import","metadata":{}},{"cell_type":"code","source":"import os\nfrom collections import defaultdict\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly import figure_factory as FF\n\nimport scipy.ndimage\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\n\nimport random\nimport pydicom\n\n\n\n# matplotlib.rcParams['figure.figsize'] = (20.0, 10.0)","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Read files","metadata":{}},{"cell_type":"code","source":"INPUT_DIR = '/kaggle/input/osic-pulmonary-fibrosis-progression'\n\ntrainset = pd.read_csv(f'{INPUT_DIR}/train.csv')\ntestset = pd.read_csv(f'{INPUT_DIR}/test.csv')\nsample_sub = pd.read_csv(f'{INPUT_DIR}/sample_submission.csv')","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# General Stats","metadata":{}},{"cell_type":"code","source":"print('train_set contained {} rows with unique {} patients'.format(len(trainset), trainset['Patient'].nunique()))\ndisplay(trainset.head())\n\nprint('test_set contained {} rows with unique {} patients'.format(len(testset), trainset['Patient'].nunique()))\ndisplay(testset.head())\n\nprint('sample_submission contained {} rows'.format(len(sample_sub)))\ndisplay(sample_sub.head())","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Summary of given data","metadata":{}},{"cell_type":"code","source":"summ = pd.DataFrame({\n    'data': ['train.csv', 'test.csv', 'sample_submission.csv'],\n    'rows': [len(trainset), len(testset), len(sample_sub)],\n    'patient': [trainset['Patient'].nunique(), testset['Patient'].nunique(), sample_sub['Patient_Week'].nunique()]\n})\nsumm.set_index('data', inplace=True)\ndisplay(summ)","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Randomly picks from each Smoking Status","metadata":{}},{"cell_type":"code","source":"def individual_patient(patient_id):\n    patient_df = trainset[trainset['Patient'] == patient_id]\n#     display(patient_df)\n    \n    fig = make_subplots(rows=2, cols=1, specs=[[{'type':'table'}], \n                                               [{'secondary_y': True}]])\n    \n    fig.add_trace(go.Table(header=dict(values=list(patient_df.columns[1:]),\n                                       align='center',\n                                       fill_color='#3c446a',\n                                       font=dict(color='white')),\n                           cells=dict(values=[patient_df[i] for i in patient_df.columns[1:]],\n                                      align='center'),\n                           columnwidth = [50, 50, 100, 50, 50, 100]), row=1, col=1)\n    \n    fig.add_trace(go.Scatter(x=patient_df['Weeks'], \n                             y=patient_df['FVC'], \n                             mode='lines+markers+text', \n                             text=patient_df['FVC'], \n                             name='FVC'), row=2, col=1, secondary_y=False)\n    \n    fig.add_trace(go.Scatter(x=patient_df['Weeks'],\n                             y=patient_df['Percent'],\n                             mode='markers', \n                             text=round(patient_df['Percent'], 2),\n                             name='Percent'), row=2, col=1, secondary_y=True)\n    \n#     fig.update_traces(textposition='top center')\n    fig.update_layout(title_text=f'<b>FVC</b> (line) and <b>Percent</b> (marker) of patient : <b>{patient_id}</b>',\n                      xaxis_title=\"Weeks\",\n                      width=800,\n                      height=700)\n    fig.update_yaxes(title_text=\"Forced vital capacity\", secondary_y=False)\n    fig.update_yaxes(title_text=\"Percent\", secondary_y=True)\n    \n    fig.show()\n    ","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for ss in trainset['SmokingStatus'].unique():\n    for sample in random.sample(trainset[trainset['SmokingStatus'] == ss]['Patient'].tolist(), 2):\n        individual_patient(sample)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Distribution of Sex and Age","metadata":{}},{"cell_type":"code","source":"fig = px.histogram(trainset, x='Age', color='Sex', marginal='box', \n                   histnorm='probability density', opacity=0.7)\nfig.update_layout(title='Distribution of Age between Male and Female',\n                  width=800, height=500)\nfig.show()","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We are going to explore further with Smoking Status feature involved. ","metadata":{}},{"cell_type":"code","source":"parti_patient = trainset.drop_duplicates(subset='Patient')\n\n\nfig = px.histogram(parti_patient,\n                  x='Age',\n                  facet_row='SmokingStatus',\n                  facet_col='Sex',\n                  )\nfig.for_each_annotation(lambda a: a.update(text=a.text.replace(\"SmokingStatus=\", \"\")))\nfig.update_layout(title='Distribution of Age sperated by Sex (col) and Smoking Status (row)',\n                  autosize=True, width=800, height=600,\n                  font_size=14)\nfig.show()","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note that at Female & Currently smokes distribution contains only two patient with 49 and 68 years old","metadata":{}},{"cell_type":"code","source":"# age of each smoking status categorized by sex\nm_exsmk_age = trainset.query('Sex == \"Male\" and SmokingStatus == \"Ex-smoker\"').drop_duplicates(subset='Patient')['Age']\nm_cursmk_age = trainset.query('Sex == \"Male\" and SmokingStatus == \"Currently smokes\"').drop_duplicates(subset='Patient')['Age']\nm_nevsmk_age = trainset.query('Sex == \"Male\" and SmokingStatus == \"Never smoked\"').drop_duplicates(subset='Patient')['Age']\n\n\nf_exsmk_age = trainset.query('Sex == \"Female\" and SmokingStatus == \"Ex-smoker\"').drop_duplicates(subset='Patient')['Age']\nf_cursmk_age = trainset.query('Sex == \"Female\" and SmokingStatus == \"Currently smokes\"').drop_duplicates(subset='Patient')['Age']\nf_nevsmk_age = trainset.query('Sex == \"Female\" and SmokingStatus == \"Never smoked\"').drop_duplicates(subset='Patient')['Age']\n\n\n# for pie chart\npie_labels = ['Male & Ex-smoker', 'Male & Currently smokes', 'Male & Never smoked','Female & Ex-smoker', 'Female & Currently smokes', 'Female & Never smoked']\nss_values = [m_exsmk_age, m_cursmk_age, m_nevsmk_age,\n             f_exsmk_age, f_cursmk_age, f_nevsmk_age]\npie_values = [*map(lambda x : len(x), ss_values)]","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = make_subplots(rows=1, cols=2, specs=[[{'type':'domain'}, {'type':'domain'}]], subplot_titles=['Male', 'Female'])\nfig.add_trace(go.Pie(labels=pie_labels[:3], values=pie_values[:3], ), row=1, col=1)\nfig.add_trace(go.Pie(labels=pie_labels[3:], values=pie_values[3:]), row=1, col=2)\nfig.update_layout(title='Rate of Ex-smoker, Never smoked, Currently smokes between Male & Female',\n                  autosize=True, width=800, height=400)\nfig.update_traces(hoverinfo='label', textinfo='percent+value')\nfig.show()","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Number of male patient : 139 (78.98%)\n- Number of female patient : 37 (21.02%)","metadata":{}},{"cell_type":"markdown","source":"## Pivot table\n Pivot table below has shown each general its statistics such as maximum value, minimum vale, mean, standard deviation\nEach value was categorized by Sex followed by SmokingStatus","metadata":{}},{"cell_type":"code","source":"pivot_smkstat_sex = pd.pivot_table(trainset, index=['Sex', 'SmokingStatus'], aggfunc={'Age': ['max', 'min', np.mean, np.std],\n                                                                                      'FVC': ['max', 'min', np.mean, np.std],\n                                                                                      'Percent':['max', 'min', np.mean, np.std]})\ndisplay(pivot_smkstat_sex)","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Included age range\nWe are going to be more specific at age range of each category.\nThe table below is same as previous one except included age range following:\n    * from 40 ~ 50\n    * from 50 ~ 60\n    * from 60 ~ 70\n    * from 70 ~ 80\n    * from 80 ~ 90\nNote that some categories might not show certain age range because it does not exist in training set","metadata":{}},{"cell_type":"code","source":"age_range = pd.cut(trainset['Age'], np.arange(40, 100, 10))\npivot_smkstat_sex = pd.pivot_table(trainset, index=['Sex', 'SmokingStatus', age_range], aggfunc={'FVC': ['max', 'min', np.mean, np.std],\n                                                                                                 'Percent':['max', 'min', np.mean, np.std]})\nmean_of_fvc = pivot_smkstat_sex['FVC'][['mean','std']].round(2).reset_index()\nmean_of_fvc = mean_of_fvc.sort_values(['Sex','SmokingStatus', 'Age'])\nmean_of_fvc['Age'] = mean_of_fvc['Age'].astype(str).map({'(40, 50]':'40_50', '(50, 60]':'50_60', '(60, 70]':'60_70',\n                                                         '(70, 80]':'70_80','(80, 90]':'80_90'})\nmean_of_fvc = mean_of_fvc.rename(columns={\"mean\":\"mean_FVC\", 'std':'std_FVC', 'Age':'AgeRange'})\ndisplay(mean_of_fvc)","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.bar(mean_of_fvc,\n             x='SmokingStatus',\n             y='mean_FVC', color='Sex',\n             barmode='group',\n             facet_col='AgeRange',\n             error_y='std_FVC',\n             category_orders={'Age':['40_50', '50_60', '60_70', '70_80', '80_90']}\n            )\nfig.update_layout(title='Mean of FVC catagorized by Sex and Age range',\n                  yaxis_title=\"Mean value of Forced vital capacity \",\n                  width=800,\n                  height=500)\nfig.show()","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note that the plot is not in order 40_50, <u><b>60_70, 50_60</b></u>, 70_80, 80_90","metadata":{}},{"cell_type":"code","source":"age_range_df = trainset.melt(id_vars=['Sex','Age','Weeks','SmokingStatus'], value_vars=['FVC'])\nage_range_df['AgeRange'] = pd.cut(age_range_df['Age'], np.arange(40, 100, 10))\nage_range_df = age_range_df.loc[:,['Sex','Age','Weeks','SmokingStatus','AgeRange','value']]\nage_range_df['AgeRange'] = age_range_df['AgeRange'].astype(str).map({'(40, 50]':'40_50', '(50, 60]':'50_60', '(60, 70]':'60_70',\n                                                                     '(70, 80]':'70_80','(80, 90]':'80_90'})","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.scatter(age_range_df,\n                 x='Weeks',\n                 y='value',\n                 facet_col='AgeRange',\n                 facet_row='SmokingStatus',\n                 category_orders={'AgeRange':['40_50', '50_60', '60_70', '70_80', '80_90']},\n                 opacity=0.5,\n                 )\nfig.for_each_annotation(lambda a: a.update(text=a.text.replace(\"SmokingStatus=\", \"\")))\nfig.update_layout(title='Scatter plot of Age range and Smoking Status')\nfig.show()","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Relationship between FVC and Percent","metadata":{}},{"cell_type":"code","source":"fig = px.density_contour(trainset,\n                         x ='Percent',\n                         y ='FVC',\n                         marginal_x=\"histogram\",\n                         marginal_y=\"histogram\",\n                         color='SmokingStatus',\n                         \n)\nfig.update_layout(title='Relationship between Percent and FVC',\n                  width=800,\n                  height=400)\nfig.show()","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is pretty obvious that the higher FVC, the higher Percent.\nAs data desciption said:\n* FVC(Forced Vital Capacity) - the recorded lung capacity in ml\n* Percent- a computed field which approximates the patient's FVC as a percent of the typical FVC for a person of similar characteristics","metadata":{}},{"cell_type":"code","source":"below_100 = trainset.query('Percent < 100')\nmore_100 = trainset.query('Percent > 100')\nbetween_5 = trainset.query('97.5 <= Percent <= 102.5')\n\nx_bar = below_100.groupby('SmokingStatus').size().index.to_list()","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = make_subplots(rows=1,\n                    cols=3,\n                    subplot_titles=['Percent < 100%', '97.5% <= Percent <= 102.5%', 'Percent > 100%'])\nfig.add_trace(go.Bar(x=x_bar,\n                     y=below_100.groupby('SmokingStatus').size(),\n                     ), row=1, col=1)\nfig.add_trace(go.Bar(x=x_bar,\n                     y=between_5.groupby('SmokingStatus').size(),\n                     ), row=1, col=2)\nfig.add_trace(go.Bar(x=x_bar,\n                     y=more_100.groupby('SmokingStatus').size(),\n                     ), row=1, col=3)\n\n\nfig.update_layout(title='Count plot of percent region',\n                  width=800,\n                  height=500,\n                  showlegend=False)","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainset.head()","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df = px.data.gapminder()\n# px.scatter(df, x=\"gdpPercap\", y=\"lifeExp\", animation_frame=\"year\", animation_group=\"country\",\n#            size=\"pop\", color=\"continent\", hover_name=\"country\",\n#            log_x=True, size_max=55, range_x=[100,100000], range_y=[25,90])\n\n","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# DICOM","metadata":{}},{"cell_type":"markdown","source":"## What is DICOM ?\n\nDICOM (Digital Imaging and Communications in Medicine) is a standard protocol for the management and transmission of medical images and related data and is used in many healthcare facilities.\n\n## What is DICOM used for ?\n\nDICOM is the international standard to communicate and manage medical images and data. Its mission is to ensure the interoperability of systems used to produce, store, share, display, send, query, process, retrieve and print medical images, as well as to manage related workflows.\n\n## Why DICOM is important ?\n\nWith the introduction of advanced imaging technologies, such as CT scans, and the growing use of computing in clinical work, ACR and NEMA saw a need for a standard method to transfer images and associated information between different vendor devices, according to the International Organization for Standardization. These devices produce a variety of digital image formats.\n\n\nsource : https://searchhealthit.techtarget.com/definition/DICOM-Digital-Imaging-and-Communications-in-Medicine","metadata":{}},{"cell_type":"markdown","source":"I am going to use DefalutDict to store all train dicom file path \n    * dicom dict = {\"{Patient_id}\": ['/kaggle/input/osic-pulmonary-fibrosis-progression/train/{Patient_id}']}\nIt helps us to easily access all different patient id dicom filepath\n\nNote: You might be familiar with some part of code down below\nSince I am new DICOM file format, I've searched for any related of similar appoarch.\nI found this extremely helpful notebook from Data Science bowl 2017 (https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial) written by Guido Zuidhof\nIf you are new to dicom like me and looking for an informative and demonstrative notebook, \nI encourage you to have a look at mentioned notebook and give him an vote as well.","metadata":{}},{"cell_type":"code","source":"DICOM_DIR = '/kaggle/input/osic-pulmonary-fibrosis-progression/train'\n\ndicom_dict = defaultdict(list)\n\n\nfor dirname in os.listdir(DICOM_DIR):\n    path = os.path.join(DICOM_DIR, dirname)\n    dicom_dict[dirname].append(path)\n    \np_id = sorted(trainset['Patient'].unique())","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# What are Hounsfield Units (HU)?\n<img src=\"https://pbrainmd.files.wordpress.com/2015/10/hounsfield-2.jpg\" width=\"500\" height=\"600\">\n\nThe Hounsfield Units (HU) make up the grayscale in medical CT imaging. It is a scale from black to white of 4096 values (12 bit) and ranges from -1024 HU to 3071 HU (zero is also a value). It is defined by the following:\n\n-1024 HU is black and represents air (in the lungs). 0 HU represents water (since we consist mostly out of water, there is a large peak here). 3071 HU is white and represents the densest tissue in a human body, tooth enamel. All other tissues are somewhere within this scale; fat is around -100 HU, muscle around 100 HU and bone spans from 200 HU (trabecular/spongeous bone) to about 2000 HU (cortical bone).\n\nMetal implants typically have very high Hounsfield units. Therefore, they are attributed the maximum value in typical 12-bit CT scans (3071).\n\nsource :\n*  https://www.materialise.com/en/faq/what-are-hounsfield-units-hu\n*  https://pbrainmd.wordpress.com/2015/10/02/hounsfield-units/","metadata":{"trusted":true}},{"cell_type":"code","source":"# ### load_scan:\n# 1. take a string path where patient dicom files were stored.\n# 2. store every slice into list and sort in ImagePositionPatient order.\n\n# ### dicom_file:\n# 1. take index number of patient which stored in dict_dicom earlier\n# 2. this might be useful when you need to pick some random patient\n# 3. It also takes specific patient Id in case you need.\n# 4. Note that this function is going to read all file in taken path.\n\n# ### get_pixels_hu\n# 1. take dicom file which had called through dicom_file function\n# 2. It stacks up all the load slices of certain patient\n# 3. stacked slices will be calculated into Hounsfield Units\n\n\n\ndef load_scan(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n        \n    return slices\n\ndef dicom_file(idx_num, patient_id=None):\n    if patient_id:\n        return load_scan(dicom_dict[patient_id][0])\n    return load_scan(dicom_dict[p_id[idx_num]][0])\n\ndef get_pixels_hu(slices):\n    image = np.stack([s.pixel_array for s in slices])\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 0\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    for slice_number in range(len(slices)):\n        \n        intercept = slices[slice_number].RescaleIntercept\n        slope = slices[slice_number].RescaleSlope\n        \n        if slope != 1:\n            image[slice_number] = slope * image[slice_number].astype(np.float64)\n            image[slbice_number] = image[slice_number].astype(np.int16)\n            \n        image[slice_number] += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = dicom_file(0)\ntest_hu = get_pixels_hu(test)\nprint('Patient {}'.format(test[0].PatientName))\nprint('Slices : {}\\nPixels : ({} x {})'.format(test_hu.shape[0], test_hu.shape[1], test_hu.shape[2]))\n","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# interactive plot is too heavy\n# fig = go.Figure()\n# fig.add_trace(go.Histogram(x=test_hu.flatten(), nbinsx=80, histnorm='percent'))\n# fig.update_layout(title='Hounsfield Units(HU) of Patient ID00007637202177411956430',\n#                   width=800,\n#                   height=600)\n\nplt.figure(figsize=(12, 8))\nax = sns.distplot(test_hu.flatten(), bins=80, norm_hist=True)\nax.set_title('Hounsfield Units of patient ID00007637202177411956430', fontsize=25)\nplt.show()","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def show_dicom_pic(dicom_pixel, p_id=None):\n    fig = plt.figure(figsize=(18, len(dicom_pixel)//2))\n    for idx, pic in enumerate(dicom_pixel):\n        fig.add_subplot(len(dicom_pixel)//5, 5, idx+1)\n        plt.imshow(pic, cmap='gray')\n        plt.axis('off')\n        plt.subplots_adjust(wspace=0.1, hspace=0.3)\n        plt.title(str(idx + 1))\n        plt.colorbar()\n    if p_id:\n        fig.suptitle('Patient {}'.format(p_id), fontsize=25)\n    fig.show()","metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_dicom_pic(test_hu)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3D Visualization","metadata":{"trusted":true}},{"cell_type":"markdown","source":"Useful links: \n\n* https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/\n* https://medium.com/@hengloose/a-comprehensive-starter-guide-to-visualizing-and-analyzing-dicom-images-in-python-7a8430fcb7ed","metadata":{}},{"cell_type":"code","source":"test2 = dicom_file(40)\ntest2_hu = get_pixels_hu(test2)\nprint('Patient {}'.format(test2[0].PatientName))\nprint('Slices : {}\\nPixels : ({} x {})'.format(test2_hu.shape[0], test2_hu.shape[1], test2_hu.shape[2]))\n\nplt.figure(figsize=(12, 8))\nax = sns.distplot(test2_hu.flatten(), bins=80, norm_hist=True)\nax.set_title('Hounsfield Units of patient {}'.format(test2[0].PatientName), fontsize=25)\nplt.show()","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = np.array([scan[0].SliceThickness, scan[0].PixelSpacing[0], scan[0].PixelSpacing[1]], dtype=np.float32)\n\n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n    \n    return image, new_spacing","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_mesh(image, threshold):\n    p = image.transpose(2, 1, 0)\n    \n    verts, faces, normals, values = measure.marching_cubes_lewiner(p, threshold)\n    return verts, faces\n\ndef static_3d(image, threshold=-300):\n    \n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    \n    verts, faces = make_mesh(image, threshold)\n    x, y, z = zip(*verts)\n    \n    mesh = Poly3DCollection(verts[faces], alpha=0.1)\n    face_color = [0.5, 0.5, 1]\n    mesh.set_facecolor(face_color)\n    \n    ax.add_collection3d(mesh)\n    ax.set_xlim(0, max(x))\n    ax.set_ylim(0, max(y))\n    ax.set_zlim(0, max(z))\n    plt.show()\n    \ndef interactive_3d(image, threshold=-300):\n    verts, faces = make_mesh(image, threshold)\n    x, y, z = zip(*verts)\n    fig = FF.create_trisurf(x=x,\n                            y=y,\n                            z=z,\n                            plot_edges=False,\n                            simplices=faces)\n    iplot(fig)","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"resampled_test2_hu, spacing = resample(test2_hu, test2)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"static_3d(resampled_test2_hu)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def largest_label_volume(im, bg=-1):\n    vals, counts = np.unique(im, return_counts=True)\n\n    counts = counts[vals != bg]\n    vals = vals[vals != bg]\n\n    if len(counts) > 0:\n        return vals[np.argmax(counts)]\n    else:\n        return None\n\ndef segment_lung_mask(image, fill_lung_structures=True):\n    \n    # not actually binary, but 1 and 2. \n    # 0 is treated as background, which we do not want\n    binary_image = np.array(image > -320, dtype=np.int8)+1\n    labels = measure.label(binary_image)\n    \n    # Pick the pixel in the very corner to determine which label is air.\n    #   Improvement: Pick multiple background labels from around the patient\n    #   More resistant to \"trays\" on which the patient lays cutting the air \n    #   around the person in half\n    background_label = labels[0,0,0]\n    \n    #Fill the air around the person\n    binary_image[background_label == labels] = 2\n    \n    \n    # Method of filling the lung structures (that is superior to something like \n    # morphological closing)\n    if fill_lung_structures:\n        # For every slice we determine the largest solid structure\n        for i, axial_slice in enumerate(binary_image):\n            axial_slice = axial_slice - 1\n            labeling = measure.label(axial_slice)\n            l_max = largest_label_volume(labeling, bg=0)\n            \n            if l_max is not None: #This slice contains some lung\n                binary_image[i][labeling != l_max] = 1\n    return binary_image","metadata":{"_kg_hide-input":true,"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"segmented_lungs = segment_lung_mask(resampled_test2_hu, False)\nsegmented_lungs_fill = segment_lung_mask(resampled_test2_hu, True)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Segmented lungs","metadata":{}},{"cell_type":"code","source":"static_3d(segmented_lungs, 1.5)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Segmented lungs filled","metadata":{}},{"cell_type":"code","source":"static_3d(segmented_lungs_fill, 1.5)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Difference between both","metadata":{}},{"cell_type":"code","source":"static_3d(segmented_lungs_fill - segmented_lungs, -0.5)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"interactive_3d(segmented_lungs_fill - segmented_lungs, -0.5)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# What's next?\n* extract valuable feature from dicom file\n","metadata":{}},{"cell_type":"code","source":"","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]}]}