{"cells":[{"metadata":{},"cell_type":"markdown","source":"<font color=\"red\" size=5><center>RSNA-STR Pulmonary Embolism Detection - 肺塞栓症検知 -</center></font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Fork元\n注釈： 本記事は https://www.kaggle.com/nitindatta/pulmonary-embolism-dicom-preprocessing-eda を和訳した物です。","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Introduction \n\nこのノートでは、DICOMファイルとCTスキャンについて学びます。また、seabornとmatplotlibを使って表形式のデータを可視化します。最後に何が追加でできるかを説明します。\n\n### **肺塞栓症(Pulmonary Embolism; PE)とは**\n*　肺塞栓症は、血液のかたまり（血栓）や、まれに他の固形物が血液の流れに乗って肺の動脈（肺動脈）に運ばれ、そこをふさいでしまう（塞栓）病気です。\n* 肺塞栓症の症状には、息切れ、特に息を吸うときの胸の痛み、血痰などがあります。\n* ほとんどの場合、肺塞栓症は、足から移動した血栓によって引き起こされます。（エコノミークラス症候群によっても肺塞栓症は引き起こされます。\n\n\n<font color=\"red\" size=3>このカーネルが役に立てたなら、Upvoteしていただければ幸いです。</font>"},{"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":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nimport matplotlib.pyplot as plt\nimport seaborn as sns\n%matplotlib inline\nfrom IPython.display import HTML\n\nsns.set_style('darkgrid')\nimport pydicom\nimport scipy.ndimage\nimport gdcm\nimport imageio\nfrom IPython import display\n\n\nfrom skimage import measure \nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nfrom skimage.morphology import disk, opening, closing\nfrom tqdm import tqdm\n\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nimport plotly.figure_factory as ff\nfrom plotly.graph_objs import *\ninit_notebook_mode(connected=True) \nfrom PIL import Image\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n\nfrom os import listdir, mkdir","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"basepath = \"../input/rsna-str-pulmonary-embolism-detection/\"\nlistdir(basepath)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(basepath + \"train.csv\")\ntest = pd.read_csv(basepath + \"test.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.shape, test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"train.head().T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# データ概要\n\n\n* `StudyInstanceUID` - データ内の各検査の一意のID。\n* `SeriesInstanceUID` - 検査中の各シリーズに固有のIDです。\n* `SOPInstanceUID` - 検査内の各画像に固有のIDです。\n* `pe_present_on_image` - 画像レベルで，画像上にPEが存在するかどうかを示します．\n* `negative_exam_for_pe` - 検査レベルで、PEが存在する画像があるかどうかを示します。\n* `qa_motion` - 検査で放射線技師がモーションアーチファクトの問題を指摘したかどうかを示します。\n* `qa_contrast` - 放射線技師が検査で造影に問題があると指摘したかどうかを示します。\n* `flow_artifact` - 参考値\n* `rv_lv_ratio_gte_1` - 検査レベル, 検査に含まれる RV/LV 比が >= 1 であるかどうかを示します.\n* `rv_lv_ratio_lt_1` - 検査レベル, 検査に含まれるRV/LV比が1未満であるかどうかを示します.\n* `leftsided_pe` - 検査中の画像の左側に PE が存在することを示します。\n* `chronic_pe` - 検査レベルの PE が慢性的なものであることを示します。\n* `true_filling_defect_not_pe` - PE ではない疾患を示します。\n* `rightsided_pe` - 試験レベルで、試験中の画像の右側に PE が存在することを示します。\n* `acute_and_chronic_pe` - 試験に含まれるPEが急性および慢性の両方であることを示します。\n* `central_pe` - 試験の画像の中心部に PE が存在することを示します。\n* `indeterminate` - 検査は PE に対して陰性ではないが、QA の問題により試験レベルの最終的なラベルセットを作成できなかったことを示しています。"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(\"Number of unique Study instances are\", train['StudyInstanceUID'].nunique())\nprint(\"Number of unique Series instances are\", train['SeriesInstanceUID'].nunique())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"📌 研究と系列の両方が同じ数であるため、各研究には1つの系列しかないと推論できます。\n\nNULL値、各列のタイプ、メモリ使用量などのデータについて、いくつかのサニティチェックを行います\n\n### 欠損値は存在するか"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('Null values in train data:',train.isnull().sum().sum())\nprint('Null values in test data:',test.isnull().sum().sum())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 欠損値なし"},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"train.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"test.info()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"約240MBのテーブルデータを用いる。"},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scans(dcm_path):\n    files = listdir(dcm_path)\n    f = [pydicom.dcmread(dcm_path + \"/\" + str(file)) for file in files]\n    return f","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example = basepath + \"train/\" + train.StudyInstanceUID.values[0] +'/'+ train.SeriesInstanceUID.values[0]\nfile_names = listdir(example)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scans = load_scans(example)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### dicomデータの例"},{"metadata":{},"cell_type":"markdown","source":"CTスキャンについて\n\n* CTスキャンは、X線を照射された物体や組織の放射線密度に関する情報を取得します。\n* 横方向のスライスは、いくつかの異なる方向から測定を行った後、スキャンを再構成します。\n* CTスキャンはすでにHUフォーマットになっています。\n* CTスキャンでは約4000個のグレー値が得られるが、それは我々の目では捉えられない。そこで、我々は\"windowing\"を実行します。\n* 水はHU 0、空気は-1000\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"scans[0]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(12,6))\nfor n in range(5):\n    image = scans[n].pixel_array.flatten()\n    rescaled_image = image * scans[n].RescaleSlope + scans[n].RescaleIntercept\n    sns.distplot(image.flatten());\nplt.title(\"HU unit distributions for 5 examples\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"📌 上のグラフでは、5つの例の画素分布をプロットしています。\n\n次に、CTスキャン画像と合わせて画素配列分布を見ていきます。"},{"metadata":{},"cell_type":"markdown","source":"## Utility Functions"},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"# dicom画像のロード\ndef load_slice(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in listdir(path)]\n    slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\n# HU配列に変換\ndef transform_to_hu(slices):\n    images = np.stack([file.pixel_array for file in slices])\n    images = images.astype(np.int16)\n\n    # convert ouside pixel-values to air:\n    # I'm using <= -1000 to be sure that other defaults are captured as well\n    images[images <= -1000] = 0\n    \n    # convert to HU\n    for n in range(len(slices)):\n        \n        intercept = slices[n].RescaleIntercept\n        slope = slices[n].RescaleSlope\n        \n        if slope != 1:\n            images[n] = slope * images[n].astype(np.float64)\n            images[n] = images[n].astype(np.int16)\n            \n        images[n] += np.int16(intercept)\n    \n    return np.array(images, dtype=np.int16)\n\ndef resample(image, scan, new_spacing=[1,1,1]):\n    spacing = np.array([float(scans_0[0].SliceThickness), \n                        float(scans_0[0].PixelSpacing[0]), \n                        float(scans_0[0].PixelSpacing[0])])\n\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)\n    \n    return image, new_spacing\n\ndef make_mesh(image, threshold=-300, step_size=1):\n    p = image.transpose(2,1,0)\n    verts, faces, norm, val = measure.marching_cubes_lewiner(p, threshold, step_size=step_size, allow_degenerate=True)\n    return verts, faces\n\n\ndef plt_3d(verts, faces):\n    print(\"Drawing\")\n    x,y,z = zip(*verts) \n    fig = plt.figure(figsize=(10, 10))\n    ax = fig.add_subplot(111, projection='3d')\n\n    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n    mesh = Poly3DCollection(verts[faces], linewidths=0.05, alpha=1)\n    face_color = [1, 1, 0.9]\n    mesh.set_facecolor(face_color)\n    ax.add_collection3d(mesh)\n\n    ax.set_xlim(0, max(x))\n    ax.set_ylim(0, max(y))\n    ax.set_zlim(0, max(z))\n#     ax.set_axis_bgcolor((0.7, 0.7, 0.7))\n    ax.set_facecolor((0.7,0.7,0.7))\n    plt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"sns.set_style('white')\nhu_scans = transform_to_hu(scans)\n\nfig, ax = plt.subplots(1,2,figsize=(15,4))\n\n\nax[0].set_title(\"CT-scan in HU\")\nax[0].imshow(hu_scans[0], cmap=\"plasma\")\nax[1].set_title(\"HU values distribution\");\nsns.distplot(hu_scans[0].flatten(), ax=ax[1],color='red', kde_kws=dict(lw=2, ls=\"--\",color='blue'));\nax[1].grid(False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"📌 グラフから、空気は約-1000のHUを持ち、次に高いのは水であることから、面積の大部分が空気で満たされていることが推察できます。HUは約0"},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient = load_slice('../input/rsna-str-pulmonary-embolism-detection/train/0003b3d648eb/d2b2960c2bbf')\nfirst_patient_pixels = transform_to_hu(first_patient)\n\ndef sample_stack(stack, rows=6, cols=6, start_with=10, show_every=5):\n    fig,ax = plt.subplots(rows,cols,figsize=[18,20])\n    for i in range(rows*cols):\n        ind = start_with + i*show_every\n        ax[int(i/rows),int(i % rows)].set_title(f'slice {ind}')\n        ax[int(i/rows),int(i % rows)].imshow(stack[ind],cmap='bone')\n        ax[int(i/rows),int(i % rows)].axis('off')\n    plt.show()\n\nsample_stack(first_patient_pixels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"順番に並べて5枚ずつ飛ばして、より多くの種類のスライスを見られるようにしています。"},{"metadata":{"trusted":true},"cell_type":"code","source":"imageio.mimsave(\"/tmp/gif.gif\", first_patient_pixels, duration=0.1)\ndisplay.Image(filename=\"/tmp/gif.gif\", format='png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first_patient_scan = '../input/rsna-str-pulmonary-embolism-detection/train/0003b3d648eb/d2b2960c2bbf'\nscans_0 = load_scans(first_patient_scan)\nimgs_after_resamp, spacing = resample(first_patient_pixels, scans_0, [1,1,1])\nv, f = make_mesh(imgs_after_resamp, threshold = 350)\nplt_3d(v, f)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"im_path = []\ntrain_path = '../input/rsna-str-pulmonary-embolism-detection/train/'\nfor i in listdir(train_path): \n    for j in listdir(train_path + i):\n        x = i+'/'+j\n        im_path.append(x)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def get_window_value(feature):\n    if type(feature) == pydicom.multival.MultiValue:\n        return np.int(feature[0])\n    else:\n        return np.int(feature)\n\npixelspacing_r = []\npixelspacing_c = []\nslice_thicknesses = []\nids = []\nid_pth = []\nrow_values = []\ncolumn_values = []\nwindow_widths = []\nwindow_levels = []\n\nfor i in im_path:\n    ids.append(i.split('/')[0]+'_'+i.split('/')[1])\n    example_dcm = listdir(train_path  + i + \"/\")[0]\n    id_pth.append(train_path + i)\n    dataset = pydicom.dcmread(train_path + i + \"/\" + example_dcm)\n    \n    window_widths.append(get_window_value(dataset.WindowWidth))\n    window_levels.append(get_window_value(dataset.WindowCenter))\n    \n    spacing = dataset.PixelSpacing\n    slice_thicknesses.append(dataset.SliceThickness)\n    \n    row_values.append(dataset.Rows)\n    column_values.append(dataset.Columns)\n    pixelspacing_r.append(spacing[0])\n    pixelspacing_c.append(spacing[1])\n    \nscan_properties = pd.DataFrame(data=ids, columns=[\"ID\"])\nscan_properties.loc[:, \"rows\"] = row_values\nscan_properties.loc[:, \"columns\"] = column_values\nscan_properties.loc[:, \"area\"] = scan_properties[\"rows\"] * scan_properties[\"columns\"]\nscan_properties.loc[:, \"pixelspacing_r\"] = pixelspacing_r\nscan_properties.loc[:, \"pixelspacing_c\"] = pixelspacing_c\nscan_properties.loc[:, \"pixelspacing_area\"] = scan_properties.pixelspacing_r * scan_properties.pixelspacing_c\nscan_properties.loc[:, \"slice_thickness\"] = slice_thicknesses\nscan_properties.loc[:, \"id_pth\"] = id_pth\nscan_properties.loc[:, \"window_width\"] = window_widths\nscan_properties.loc[:, \"window_level\"] = window_levels\nscan_properties.head().T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### pixelspacing(ピクセル間隔)\n* dicom ファイルにある pixelspacing 属性は重要な属性です。これは、1つのピクセルがどのくらいの物理的な距離をカバーしているかを教えてくれます。横断スライスの平面内のX方向とY方向を記述する2つの値しかないことがわかります。\n* 1人の患者の場合、このピクセル間隔は通常、すべてのスライスで同じです。\n* しかし、患者の間では医師やクリニックの個人的または組織的な好みによってピクセル間隔が異なる場合があり、またスキャナーのタイプにも依存します。その結果、肺のサイズで2つの画像を比較した場合、自動的に大きい方が臓器の物理的なサイズが大きいことを意味するわけではありません。"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sns.set_style('darkgrid')\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(pixelspacing_r, ax=ax[0], color='green', kde_kws=dict(lw=3, ls=\"--\",color='red'))\nax[0].set_title(\"Pixel spacing distribution \\n in row direction \")\nax[0].set_ylabel(\"Counts in train\")\nax[0].set_xlabel(\"mm\")\nsns.distplot(pixelspacing_c, ax=ax[1], color=\"Blue\",kde_kws=dict(lw=3, ls=\"--\",color='red'))\nax[1].set_title(\"Pixel spacing distribution \\n in column direction\");\nax[1].set_ylabel(\"Counts in train\");\nax[1].set_xlabel(\"mm\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"我々は、値が本当に患者から患者に多くの違いがあることを見ることができます! 彼らはmmで与えられているように、CTスキャンは通常、512行と列の値をカバーしています。"},{"metadata":{},"cell_type":"markdown","source":"### 1回のCTスキャンでカバーされる物理的な領域とスライス量\n\nさて、ct-scanでカバーされる物理的な距離を計算するためのいくつかの重要な量がわかりました!"},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties[\"r_distance\"] = scan_properties.pixelspacing_r * scan_properties.rows\nscan_properties[\"c_distance\"] = scan_properties.pixelspacing_c * scan_properties[\"columns\"]\nscan_properties[\"area_cm2\"] = 0.1* scan_properties[\"r_distance\"] * 0.1*scan_properties[\"c_distance\"]\nscan_properties[\"slice_volume_cm3\"] = 0.1*scan_properties.slice_thickness * scan_properties.area_cm2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(scan_properties.area_cm2, ax=ax[0], color=\"Limegreen\",kde_kws=dict(lw=3, ls=\"--\",color='red'))\nsns.distplot(scan_properties.slice_volume_cm3, ax=ax[1], color=\"Mediumseagreen\",kde_kws=dict(lw=3, ls=\"--\",color='red'))\nax[0].set_title(\"CT-slice area in $cm^{2}$\")\nax[1].set_title(\"CT-slice volume in $cm^{3}$\")\nax[0].set_xlabel(\"$cm^{2}$\")\nax[1].set_xlabel(\"$cm^{3}$\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties.head(3).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scan_properties.describe().T","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"scan_properties.to_csv('Pulmonary_Embolism_CT_scans_data.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"scan_cols = scan_properties.copy()\nscan_cols.drop(['rows','columns','area'],axis=1,inplace=True)\n\ncorr = scan_cols.corr()\nmask = np.zeros_like(corr)\nmask[np.triu_indices_from(mask)] = True\nwith sns.axes_style(\"white\"):\n    f, ax = plt.subplots(figsize=(10, 10))\n    ax = sns.heatmap(corr,mask=mask,square=True,linewidths=.8,cmap=\"viridis\",annot=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"異常なヒートマップが見られますが、これは値のほとんどが、いくつかの特徴を高度に相関させるようにして、得られた値だからです。"},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = train.copy()\ncols.drop(['StudyInstanceUID','SeriesInstanceUID','SOPInstanceUID'],axis=1,inplace=True)\ncolumns = cols.columns","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(7,2,figsize=(16,28))\nfor i,col in enumerate(columns): \n    plt.subplot(7,2,i+1)\n    sns.countplot(cols[col],palette='hot')   ","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"corr = cols.corr()\nmask = np.zeros_like(corr)\nmask[np.triu_indices_from(mask)] = True\nwith sns.axes_style(\"white\"):\n    f, ax = plt.subplots(figsize=(12, 12))\n    ax = sns.heatmap(corr,mask=mask,square=True,linewidths=.8,cmap=\"summer\",annot=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Acknowledgements\n1. [Excellent work of Laura Fink](https://www.kaggle.com/allunia/pulmonary-fibrosis-dicom-preprocessing)\n2. [Insights used by prk007](https://www.kaggle.com/prk007/insights-from-tabular-and-image-data)\n3. [3D reconstruction by Md. Redwan Karim Sony](https://www.kaggle.com/redwankarimsony/rsna-str-3d-stacking-3d-plot-segmentation/comments)"},{"metadata":{},"cell_type":"markdown","source":""}],"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}