{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook aggregates the shape of the image and checks the actual image.  \nThis shows that the margins at the edges of the image are not consistent and need to be considered for pre-processing.\n\nupdate:  \nver4: Data that takes a long time to create has been changed to be read from kaggle [Datasets](https://www.kaggle.com/currypurin/osic-image-eda).  \nver7: added crop preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport pickle\nfrom pathlib import Path\nimport gc\nimport pandas as pd\nimport numpy as np\nfrom PIL import Image\nimport matplotlib.pyplot as plt\nfrom tqdm.auto import tqdm\nimport pydicom\nfrom pydicom.tag import Tag\nimport gc\n\npd.options.display.max_rows=200\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"INPUT = Path(\"../input/osic-pulmonary-fibrosis-progression/\")\ntrain = pd.read_csv(INPUT / 'train.csv')\ndataset_dir = Path(\"../input/osic-image-eda/\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ls ../input/osic-image-eda/ ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# name and number of dicoms"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_n_dicom_df(train):\n    df_list = []\n    for patient_id in train['Patient'].unique():\n        patient_dir = INPUT / 'train' / patient_id\n        path_list = list(patient_dir.glob(\"*\"))\n        n_dicom = len(path_list)\n        n_list = [int(str(i).split('/')[-1].split('.')[0]) for i in path_list]\n        sort_n_list = sorted(n_list)\n        tmp_df = pd.DataFrame({'Patient': [patient_id],\n                               'n_dicom': [n_dicom],\n                               'n_list': [sort_n_list]})\n        df_list.append(tmp_df)\n        \n    n_dicom_df = pd.concat(df_list, sort=False)\n    return n_dicom_df.reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_dicom_df = get_n_dicom_df(train)\nn_dicom_df.to_csv('n_dicom_df.csv', index=False)\nn_dicom_df.head(15)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The number of dicoms is different for each patient. Also, the file may not start with one."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(n_dicom_df['n_dicom'], bins=20)\nplt.title('Number of dicom per patient');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_dicom_df['n_dicom'].value_counts().head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# image shape"},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df_path = dataset_dir / \"shape_df.csv\"\nif shape_df_path.is_file():\n    shape_df = pd.read_csv(shape_df_path)\nelse:\n    !conda install -c conda-forge gdcm -y\n    gr = train.groupby('Patient')\n    df_list = []\n    for patient_id, group_df in tqdm(gr):\n        height_list = []\n        width_list = []\n        shape_list = []\n        tmp_df_list = []\n        for dcm_path in (INPUT / 'train' / patient_id).glob(\"*\"):\n            try:\n                dicom = pydicom.dcmread(dcm_path)\n                tmp_df = pd.DataFrame({'Patient': [patient_id],\n                                       'height': [dicom[Tag(\"Rows\")].value],\n                                       'width': [dicom[Tag(\"Columns\")].value],\n                                       'shape': [str(dicom.pixel_array.shape)]})\n                tmp_df_list.append(tmp_df)\n            except:\n                print(dcm_path)\n        if len(tmp_df_list) >= 1:\n            df = pd.concat(tmp_df_list)\n            df.drop_duplicates(inplace=True)\n            df_list.append(df.reset_index(drop=True))\n    shape_df = pd.concat(df_list)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df.to_csv('shape_df.csv', index=False)\nshape_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df.groupby('Patient').count().max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_list = shape_df['shape'].value_counts().index\nshape_df['shape'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"shapes with different height and width  \n(752, 888) , (734, 888) , (843, 888) , (733, 888) , (1100, 888) ,(788, 888)  "},{"metadata":{},"cell_type":"markdown","source":"# img"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ndef imshow_dcm(height, width):\n    _shape_df = shape_df[shape_df['shape'] == str((height, width))]\n    dcm_list = []\n    path_list = []\n    for i in range(4):\n        patient_id = np.random.choice(_shape_df['Patient'])\n        dcm_dir = INPUT / f'train/{patient_id}'\n        dcm_path = np.random.choice(list(dcm_dir.glob(\"*\")))\n        path_list.append(str(dcm_path).split('/')[-2:])\n        dicom = pydicom.dcmread(dcm_path)\n        dcm_list.append(dicom)\n    for i in range(4):\n        plt.subplot(2, 2, i+1)\n        plt.imshow(dcm_list[i].pixel_array, cmap=plt.cm.bone)\n        plt.title(path_list[i])\n        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 512 x 512"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(512, 512)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (768, 768)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(768, 768)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (752, 888)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(752, 888)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (632, 632)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(632, 632)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (734, 888)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(734, 888)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (843, 888)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(843, 888)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (733, 888)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(733, 888)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (1100, 888)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(1100, 888)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (1302, 1302)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(1302, 1302)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## (788, 888)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nimshow_dcm(788, 888)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The ones with different height and width seem to have margins and need to be cropped. Let's try cropping them."},{"metadata":{},"cell_type":"markdown","source":"# Crop"},{"metadata":{"trusted":true},"cell_type":"code","source":"dcm_path = Path(\"../input/osic-pulmonary-fibrosis-progression/train/ID00094637202205333947361/8.dcm\")\ndicom = pydicom.dcmread(dcm_path)\nimg = dicom.pixel_array\nplt.imshow(img, cmap=plt.cm.bone)\nplt.title(f\"shape: {img.shape}\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Areas with the same number of pixels on the edges are not required. Crop it.\n\ndef crop_image(img: np.ndarray):\n    edge_pixel_value = img[0, 0]\n    mask = img != edge_pixel_value\n    return img[np.ix_(mask.any(1),mask.any(0))]\n\nplt.figure(figsize=(16, 8))\nplt.subplot(121)\nplt.imshow(img, cmap=plt.cm.bone)\nplt.title(img.shape)\n\nplt.subplot(122)\nplt.imshow(crop_image(img), cmap=plt.cm.bone)\nplt.title(crop_image(img).shape);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This process is based on codes from this great notebook https://www.kaggle.com/ratthachat/aptos-eye-preprocessing-in-diabetic-retinopathy.  \nLet's crop the other shapes"},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df = shape_df.merge(n_dicom_df, on='Patient', how='left')\ncrop_df = shape_df[shape_df[\"shape\"].isin([\"(752, 888)\" , \"(734, 888)\" , \"(843, 888)\", \"(733, 888)\" , \"(1100, 888)\", \"(788, 888)\"])]\ncrop_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Display the 1.dicom of each Patient.\n\ndef get_image_array_from_dicom(patient_id, n):\n    dcm_file_path = INPUT / f'train/{patient_id}/{n}.dcm'\n    dicom = pydicom.dcmread(dcm_file_path)\n    return dicom.pixel_array\n\nfor patient_id in crop_df['Patient']:\n    n = 1\n    image = get_image_array_from_dicom(patient_id, n)\n    plt.figure(figsize=(16, 8))\n    plt.subplot(121)\n    plt.imshow(image, cmap=plt.cm.bone)\n    plt.title(image.shape)\n    \n    plt.subplot(122)\n    if image.shape[0] != image.shape[1]:\n        image = crop_image(image)\n    plt.imshow(image, cmap=plt.cm.bone)\n    plt.title(image.shape)    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# preprocess"},{"metadata":{},"cell_type":"markdown","source":"First, check the thickness between the images."},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_image_position_diff(patient_id, n_list):\n    dicom_path_list = [INPUT / 'train' / patient_id / f'{n}.dcm' for n in n_list]\n    dicoms = [pydicom.read_file(path_) for path_ in dicom_path_list]\n    diff_list_ = []\n    for i in range(len(dicoms)-1):\n        try:\n            diff = np.abs(dicoms[i].ImagePositionPatient[2] - dicoms[i + 1].ImagePositionPatient[2])\n        except AttributeError:\n            diff = np.nan\n        diff_list_.append(diff)\n    return diff_list_\n\ndef get_slicethickness_slope_intercept(patient_id, n_list):\n    dicom_path_list = [INPUT / 'train' / patient_id / f'{n}.dcm' for n in n_list]\n    dicoms = [pydicom.read_file(path_) for path_ in dicom_path_list]\n    thickness_list_ = []\n    pixelspacing_list_ = []\n    slope_list_ = []\n    intercept_list_ = []\n    for i in range(len(dicoms)):\n        try:\n            slice_thickness_ = dicoms[i].SliceThickness\n        except AttributeError:\n            slice_thickness_ = np.nan\n        try:\n            pixelspacing_ = dicoms[i].PixelSpacing\n        except AttributeError:\n            pixelspacing_ = np.nan\n        try:\n            slope_ = dicoms[i].RescaleSlope\n        except AttributeError:\n            slope_ = np.nan\n        try:\n            intercept_ = dicoms[i].RescaleIntercept\n        except AttributeError:\n            intercept_ = np.nan\n\n        thickness_list_.append(slice_thickness_)\n        pixelspacing_list_.append(pixelspacing_)\n        slope_list_.append(slope_)\n        intercept_list_.append(intercept_)\n        \n    return thickness_list_, pixelspacing_list_, slope_list_, intercept_list_","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"diff_list = []\nfor i in tqdm(range(len(shape_df))):\n    diff_list.append(get_image_position_diff(shape_df.loc[i, 'Patient'], shape_df.loc[i, 'n_list']))\n\nthickness_list = []\npixelspacing_list = []\nslope_list = []\nintercept_list = []\nfor i in tqdm(range(len(shape_df))):\n    t, p, s, i = get_slicethickness_slope_intercept(shape_df.loc[i, 'Patient'], shape_df.loc[i, 'n_list'])\n    thickness_list.append(t)\n    pixelspacing_list.append(p)\n    slope_list.append(s)\n    intercept_list.append(i)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df['diff_list'] = diff_list\nshape_df['thickness_list'] = thickness_list\nshape_df['pixelspacing_list'] = pixelspacing_list\nshape_df['slope_list'] = slope_list\nshape_df['intercept_list'] = intercept_list","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## diff"},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df['diff_list_std'] = shape_df['diff_list'].apply(lambda x: np.array(x).std())\nshape_df['diff_list_n_nan'] = shape_df['diff_list'].apply(lambda x: np.sum(pd.Series(x).isna()))\n\nshape_df[['Patient', 'diff_list', 'diff_list_std', 'diff_list_n_nan', 'thickness_list', 'pixelspacing_list']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## slope and intercept"},{"metadata":{"trusted":true},"cell_type":"code","source":"shape_df[['Patient', 'slope_list', 'intercept_list']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.all(shape_df['slope_list'].apply(lambda x:np.all(x)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with open('shape_df_ver2.pickle', 'wb') as f:\n    pickle.dump(shape_df, f)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Todo\n\n* Create a feature from a dicom image."},{"metadata":{},"cell_type":"markdown","source":"# References:\n\n1. https://www.kaggle.com/ratthachat/aptos-eye-preprocessing-in-diabetic-retinopathy\n2. https://www.kaggle.com/jameschapman19/pytorch-tabular-qr-histogram"}],"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}