{"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":"code","source":"'''\nThis notebook can run in local as long as PATH is set properly \npointing to the data folder.\nThis notebook saves intermidiate results in PATH folder, \nso write permission to PATH folder is needed to run this notebook.\n\nIn kaggle, we don't have write permission. So this notebook can't run in kaggle.\nIt uses RandomForestClassifier from scikit-learn to train a model from nii files.\nThen the model is used to predict the spine location for each dcm file. \nUsing the spine location for each dcm file with train.csv, each dcm file can be labelled\nwith C1, C2, ...C7. \nOnce each dcm file is labelled, Conv2D can be used to train a model. \nWhich architecture to use is not decided yet. Should be 1 input 1 outputs, \nor 1 input 7 outputs, or 1 input 8 outputs? \nOnce the work is done, it will be posted here for comments. \nPlease comment! It is great opportunity to learn from experts!\n'''\n\n'''\nload train.csv\n'''\n%matplotlib inline\n\nimport os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n#set up PATH to the downloaded files in your environment\n#this is in my local\n# PATH = os.path.expanduser('~')\n# PATH = os.path.join(PATH, 'Downloads','rsna-2022-cervical-spine-fracture-detection')\n\nPATH = '../input/rsna-2022-cervical-spine-fracture-detection'\ndf_train = pd.read_csv(os.path.join(PATH, \"train.csv\"))\ndisplay(df_train.head())","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-08-19T14:28:54.078411Z","iopub.execute_input":"2022-08-19T14:28:54.078931Z","iopub.status.idle":"2022-08-19T14:28:54.115538Z","shell.execute_reply.started":"2022-08-19T14:28:54.078894Z","shell.execute_reply":"2022-08-19T14:28:54.114691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ndraw pie chart \n'''\nax = df_train[\"patient_overall\"].value_counts().plot.pie(ylabel=\"patient_overall\", autopct=\"%.1f%%\")","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:33:24.707068Z","iopub.execute_input":"2022-08-19T14:33:24.707471Z","iopub.status.idle":"2022-08-19T14:33:24.874183Z","shell.execute_reply.started":"2022-08-19T14:33:24.707422Z","shell.execute_reply":"2022-08-19T14:33:24.873024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ndraw bar chart for C1, C2, ...C7\n'''\ncx = [f\"C{i + 1}\" for i in range(7)]\nax = df_train[cx].sum().plot.bar(xlabel=\"spine\", ylabel=\"counts\", grid=True)\nprint(cx)","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:33:47.943165Z","iopub.execute_input":"2022-08-19T14:33:47.94366Z","iopub.status.idle":"2022-08-19T14:33:48.247404Z","shell.execute_reply.started":"2022-08-19T14:33:47.94362Z","shell.execute_reply":"2022-08-19T14:33:48.246226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ndraw bar chart for occurances\n'''\nax = df_train[cx].sum(axis=1).value_counts().plot.bar(xlabel=\"occurances\", ylabel=\"counts\", grid=True)","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:34:19.903593Z","iopub.execute_input":"2022-08-19T14:34:19.904628Z","iopub.status.idle":"2022-08-19T14:34:20.096099Z","shell.execute_reply.started":"2022-08-19T14:34:19.904573Z","shell.execute_reply":"2022-08-19T14:34:20.094737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ndraw heat map for C1, C2...C7\n'''\nimport seaborn as sn\ncorr = df_train[cx].corr()\nax = sn.heatmap(corr, cmap=\"Blues\", annot=True)","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:34:42.127924Z","iopub.execute_input":"2022-08-19T14:34:42.128337Z","iopub.status.idle":"2022-08-19T14:34:43.047725Z","shell.execute_reply.started":"2022-08-19T14:34:42.128292Z","shell.execute_reply":"2022-08-19T14:34:43.046382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nload train_bounding_boxes.csv\n'''\ndf_bounding = pd.read_csv(os.path.join(PATH, \"train_bounding_boxes.csv\"))\ndisplay(df_bounding.head())","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:35:02.200585Z","iopub.execute_input":"2022-08-19T14:35:02.200976Z","iopub.status.idle":"2022-08-19T14:35:02.236309Z","shell.execute_reply.started":"2022-08-19T14:35:02.20094Z","shell.execute_reply":"2022-08-19T14:35:02.235101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nload 1 dcm file to check its attributes\n'''\nimport pydicom\nfrom pydicom.pixel_data_handlers import apply_voi_lut\nspl = df_bounding.sample(1).iloc[0]\ndicom_path = os.path.join(PATH, \"train_images\", spl[\"StudyInstanceUID\"], f'{spl[\"slice_number\"]}.dcm')\ndicom = pydicom.dcmread(dicom_path)\ndisplay(dicom)\ndisplay(dicom.dir())\nimg = apply_voi_lut(dicom.pixel_array, dicom)","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:36:53.094535Z","iopub.execute_input":"2022-08-19T14:36:53.095847Z","iopub.status.idle":"2022-08-19T14:36:53.460114Z","shell.execute_reply.started":"2022-08-19T14:36:53.095786Z","shell.execute_reply":"2022-08-19T14:36:53.458927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ndisplay dcm with patch from train_bounding_boxes.csv\n'''\nfrom matplotlib.patches import Rectangle\n\nfig, ax = plt.subplots()\nax_im = ax.imshow(img, cmap=\"gray\")\nrect = Rectangle((spl['x'], spl['y']), spl['width'], spl['height'], linewidth=1, edgecolor='r', facecolor='none')\n# Add the patch to the Axes\nax.add_patch(rect)\n\n_= plt.colorbar(ax_im)","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:37:32.865518Z","iopub.execute_input":"2022-08-19T14:37:32.865977Z","iopub.status.idle":"2022-08-19T14:37:33.186363Z","shell.execute_reply.started":"2022-08-19T14:37:32.865941Z","shell.execute_reply":"2022-08-19T14:37:33.184909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nload 1 nii file \n'''\nimport nibabel as nib\nniiImg = nib.load(os.path.join(PATH, \"segmentations\", \"1.2.826.0.1.3680043.9926.nii\"))\nseg = niiImg.get_fdata()\nseg = seg[:,::-1,::-1].transpose(2,1,0)\ndisplay(seg.shape)\ndisplay(seg[169])\nnp.unique(seg[169])","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:37:57.024228Z","iopub.execute_input":"2022-08-19T14:37:57.025577Z","iopub.status.idle":"2022-08-19T14:37:57.906995Z","shell.execute_reply.started":"2022-08-19T14:37:57.025532Z","shell.execute_reply":"2022-08-19T14:37:57.905612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ndisplay nii file\n'''\nfig, axes = plt.subplots(4,4, figsize=(12,12))\nfor i, ax in enumerate(axes.reshape(-1)):\n    ax.imshow(seg[59+i])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:38:17.470812Z","iopub.execute_input":"2022-08-19T14:38:17.47132Z","iopub.status.idle":"2022-08-19T14:38:19.528246Z","shell.execute_reply.started":"2022-08-19T14:38:17.471283Z","shell.execute_reply":"2022-08-19T14:38:19.526865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nmake a dataframe for all nii files, consitting of StudyInstanceUID and path to nii file\n'''\nseg_paths = glob.glob(os.path.join(PATH, \"segmentations\", f\"*.nii\"))\nseg_df = pd.DataFrame({'path': seg_paths})\nseg_df['StudyInstanceUID'] = seg_df['path'].apply(lambda x:x.split('/')[-1][:-4])\nseg_df = seg_df[['StudyInstanceUID','path']]\nprint('seg_df shape:', seg_df.shape)\nseg_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:38:46.70217Z","iopub.execute_input":"2022-08-19T14:38:46.703064Z","iopub.status.idle":"2022-08-19T14:38:46.72477Z","shell.execute_reply.started":"2022-08-19T14:38:46.703001Z","shell.execute_reply":"2022-08-19T14:38:46.723503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nGet information from the .dcm files\n'''\n\n# From https://www.kaggle.com/code/andradaolteanu/rsna-fracture-detection-dicom-images-explore\ndef get_observation_data(path):\n    dataset = pydicom.read_file(path)\n    # Dictionary to store the information from the image\n    observation_data = {\n        \"Rows\" : dataset.get(\"Rows\"),\n        \"Columns\" : dataset.get(\"Columns\"),\n        \"SOPInstanceUID\" : dataset.get(\"SOPInstanceUID\"),\n        \"ContentDate\" : dataset.get(\"ContentDate\"),\n        \"SliceThickness\" : dataset.get(\"SliceThickness\"),\n        \"InstanceNumber\" : dataset.get(\"InstanceNumber\"),\n        \"ImagePositionPatient\" : dataset.get(\"ImagePositionPatient\"),\n        \"ImageOrientationPatient\" : dataset.get(\"ImageOrientationPatient\"),\n    }\n    # String columns\n    str_columns = [\"SOPInstanceUID\", \"ContentDate\", \n                   \"SliceThickness\", \"InstanceNumber\"]\n    for k in str_columns:\n        observation_data[k] = str(dataset.get(k)) if k in dataset else None\n\n    return observation_data","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:39:10.777326Z","iopub.execute_input":"2022-08-19T14:39:10.77816Z","iopub.status.idle":"2022-08-19T14:39:10.787679Z","shell.execute_reply.started":"2022-08-19T14:39:10.778102Z","shell.execute_reply":"2022-08-19T14:39:10.786275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ntest get_observation_data() with 1 dcm file\n'''\npath = os.path.join(PATH, \"train_images/1.2.826.0.1.3680043.10001/109.dcm\")\nexample = get_observation_data(path)\nprint(example)\ndisplay(example.keys())","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:39:28.56186Z","iopub.execute_input":"2022-08-19T14:39:28.562263Z","iopub.status.idle":"2022-08-19T14:39:28.584627Z","shell.execute_reply.started":"2022-08-19T14:39:28.56222Z","shell.execute_reply":"2022-08-19T14:39:28.583521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nRetrieves the desired metadata from the .dcm files and saves it into dataframe.\n'''\n\nfrom tqdm import tqdm\ndef get_metadata():\n    \n    exceptions = 0\n    dicts = []\n\n    for k in tqdm(range(len(df_train))):\n        if (k % 100)==0:\n            print(f'Iteration: {k}')\n            \n        dt = df_train.iloc[k, :]\n\n        # Get all .dcm paths for this Instance\n        dcm_paths = glob.glob(os.path.join(PATH, \"train_images\", f\"{dt.StudyInstanceUID}/*\"))\n        for path in dcm_paths:\n            try:\n                # Get datasets\n                dataset = get_observation_data(path)\n                dicts.append(dataset)\n            except Exception as e:\n                exceptions += 1\n                continue\n\n    # Convert into df\n    meta_train_data = pd.DataFrame(data=dicts, columns=example.keys())\n    \n    # Export information\n    meta_train_data.to_csv(os.path.join(PATH, \"meta_train.csv\"), index=False)\n    \n    print(f\"Metadata created. Number of total fails: {exceptions}.\")","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:57:47.395361Z","iopub.execute_input":"2022-08-19T14:57:47.395828Z","iopub.status.idle":"2022-08-19T14:57:47.405345Z","shell.execute_reply.started":"2022-08-19T14:57:47.395794Z","shell.execute_reply":"2022-08-19T14:57:47.404493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ngenerate meta_train.csv file. it takes time\n'''\nget_metadata()","metadata":{"execution":{"iopub.status.busy":"2022-08-19T14:57:53.653805Z","iopub.execute_input":"2022-08-19T14:57:53.654351Z","iopub.status.idle":"2022-08-19T17:18:45.189812Z","shell.execute_reply.started":"2022-08-19T14:57:53.65431Z","shell.execute_reply":"2022-08-19T17:18:45.186987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nclean up meta_train.csv and save into meta_train_clean.csv\n'''\nmeta_train = pd.read_csv(os.path.join(PATH,\"meta_train.csv\"))\nmeta_train[\"StudyInstanceUID\"] = meta_train[\"SOPInstanceUID\"].apply(lambda x: \".\".join(x.split(\".\")[:-2]))\n# # Create column with full image size (height x width)\nmeta_train[\"ImageSize\"] = meta_train[\"Rows\"].astype(str) + \" x \" + meta_train[\"Columns\"].astype(str)\nmeta_train.drop(columns=\"ContentDate\", axis=1, inplace=True)\n\nmeta_train.to_csv(os.path.join(PATH,\"meta_train_clean.csv\"), index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nload meta_train_clean.csv for those in seg_df,which is loaded from segmentation folder \n'''\n# Metadata was extracted previously (check out my RSNA dataset)\nmeta_train = pd.read_csv(os.path.join(PATH,\"meta_train_clean.csv\"))\n\n# Only select patients with segmentations\nmeta_seg = meta_train[meta_train['StudyInstanceUID'].isin(seg_df['StudyInstanceUID'])].reset_index(drop=True)\nprint('meta_seg shape:', meta_seg.shape)\nmeta_seg.head(2)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\npopulate C1,C2...C7 and save into meta_segmentation.csv\n'''\nk=0\n# Loop over 87 patients with segmentations\nfor path, UID in zip(seg_df['path'], seg_df['StudyInstanceUID']):\n    # Get segmentations for patient\n    seg_nib = nib.load(path)\n    seg = seg_nib.get_fdata()\n    seg = seg[:, ::-1, ::-1].transpose(2, 1, 0) # Align orientation with train images\n    num_slices, _, _ = seg.shape\n    \n    # Loop over slices\n    for i in range(num_slices):\n        mask = seg[i]\n        unique_vals = np.unique(mask)\n        \n        # Loop over unique values (except 0)\n        for j in unique_vals[1:]:\n            \n            # Ignore thoratic spine etc\n            if j <= 7:   \n#                meta_seg.loc[(meta_seg['StudyInstanceUID']==UID)&(meta_seg['Slice']==i),f'C{int(j)}'] = 1\n                meta_seg.loc[(meta_seg['StudyInstanceUID']==UID)&(meta_seg['InstanceNumber']==i),f'C{int(j)}'] = 1\n                \n    # Iteration tracker\n    if (k%10)==0:\n        print(f'Iteration:{k}')\n    k+=1\n\n# Save extracted targets\nmeta_seg.to_csv(os.path.join(PATH,\"meta_segmentation.csv\"), index=False)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nload meta_segmentation.csv and load ImagePositionPatient as list\n'''\nmeta_seg = pd.read_csv(os.path.join(PATH,'meta_segmentation.csv'),converters={'ImagePositionPatient': pd.eval})\nmeta_seg.head(1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\ncheck content of meta_seg\n'''\npd.set_option('display.max_rows', None)\nmeta_seg[meta_seg['StudyInstanceUID']==meta_seg['StudyInstanceUID'].unique()[0]].loc[110:340]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nclean up data for meta_seg and assign it to df4\n'''\nimport pandas as pd\ndf2 = pd.DataFrame(meta_seg['ImagePositionPatient'].tolist(),columns=['ImagePositionPatient_x', 'ImagePositionPatient_y', 'ImagePositionPatient_z'])\ndf3 = pd.concat([df2, meta_seg], axis=1)\ndf4 = df3.drop('ImagePositionPatient',axis=1)\ndf4['C1'].fillna(0,inplace=True)\ndf4['C2'].fillna(0,inplace=True)\ndf4['C3'].fillna(0,inplace=True)\ndf4['C4'].fillna(0,inplace=True)\ndf4['C5'].fillna(0,inplace=True)\ndf4['C6'].fillna(0,inplace=True)\ndf4['C7'].fillna(0,inplace=True)\ndf4.shape\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nsplit data into train set and validation set\n'''\nfrom sklearn.model_selection import train_test_split\n\nfeatures = ['InstanceNumber','SliceThickness','ImagePositionPatient_x','ImagePositionPatient_y','ImagePositionPatient_z']\ntargets = ['C1','C2','C3','C4','C5','C6','C7']\n\ndf4 = df4.reindex(columns=[*features,*targets])\ndisplay(df4.head(2))\n# Features and targets\nX = df4[features]\ny = df4[targets]\n\n# Train-validation split\nX_train, X_valid, y_train, y_valid = train_test_split(X,y,train_size=0.8,test_size=0.2,random_state=0)\nX_train.head(1)\n\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nuse RandomForestClassifier to train a model to predict location of dcm in the spine\n'''\nfrom sklearn.ensemble import RandomForestClassifier\n\n# Classifier\nclf = RandomForestClassifier()\nclf.fit(X_train, y_train)\n\n# Score model\nprint('Classifier average accuracy:', clf.score(X_valid,y_valid))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\nFeature importances\n'''\npd.DataFrame({'Feature':features, 'Importance':clf.feature_importances_}).sort_values(by='Importance', ascending=False)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\npredict a few to test\n'''\nnp.set_printoptions(threshold=np.inf)\nclf.predict(X)[110:120,:]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\npredict targets for train set\n'''\n# Read in metadata for entire train set\nmeta_train2 = pd.read_csv(os.path.join(PATH,'meta_train_clean.csv'),converters={'ImagePositionPatient': pd.eval})\nmeta_train3 = pd.DataFrame(meta_train2[\"ImagePositionPatient\"].tolist(),columns=['ImagePositionPatient_x', 'ImagePositionPatient_y', 'ImagePositionPatient_z'])\n\nmeta_train4 = pd.concat([meta_train3, meta_train2], axis=1)\nmeta_train5 = meta_train4[['InstanceNumber','SliceThickness','ImagePositionPatient_x','ImagePositionPatient_y','ImagePositionPatient_z']]\n\n# Initialise targets\nmeta_train5[targets]=0\n\n# Predict targets for entire train set\nmeta_train5[targets] = clf.predict(meta_train5[features])\ndisplay(meta_train5.head(10))\n\n# Save to csv\nmeta_train5.to_csv(os.path.join(PATH, 'meta_train_with_vertebrae.csv'), index=False)\n\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"meta_train_with_vertebrae.csv has column: C1,C2...C7. From the file, we can see some dcm files, all C1,C2...C7 are 0. Some has only one 1 for C1,C2...C7. Some has two 1 for C1,C2...C7, but the two are adjacent. Looks like the prediction makes sense, at least. ","metadata":{}}]}