{"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":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        #print(os.path.join(dirname, filename))\n        pass\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-05-09T16:05:19.943548Z","iopub.execute_input":"2022-05-09T16:05:19.944338Z","iopub.status.idle":"2022-05-09T16:05:22.376160Z","shell.execute_reply.started":"2022-05-09T16:05:19.944236Z","shell.execute_reply":"2022-05-09T16:05:22.375439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading and preprocessing data","metadata":{}},{"cell_type":"markdown","source":"First we need to convert dicom files to objects that we can use; CSV files, jpg or png images e.g. To do this, we can use pydicom library for reading dicom files and then converting them. For more information you can visit these links:\n\n[https://asvcode.github.io/MedicalImaging/medical_imaging/dicom/fastai/2020/04/28/Medical-Imaging-Using-Fastai.html](http://)\n\n[https://pydicom.github.io/pydicom/stable/tutorials/dataset_basics.html](http://)","metadata":{}},{"cell_type":"code","source":"from fastai.basics import *\nfrom fastai.callback.all import *\nfrom fastai.vision.all import *\nfrom fastai.medical.imaging import *\n\nimport pydicom, kornia, cv2\n\ntrain_dcm_path = get_dicom_files(\"../input/unifesp-x-ray-body-part-classifier/train\")\ntest_dcm_path = get_dicom_files(\"../input/unifesp-x-ray-body-part-classifier/test\")\n\nprint(\"train len:\", len(train_dcm_path), \"\\ntest len:\", len(test_dcm_path))","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:22.377750Z","iopub.execute_input":"2022-05-09T16:05:22.378001Z","iopub.status.idle":"2022-05-09T16:05:28.141009Z","shell.execute_reply.started":"2022-05-09T16:05:22.377968Z","shell.execute_reply":"2022-05-09T16:05:28.140198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can convert some parts of dicom file to pandas dataset. Thanks to Mr.Bombonato, he convert these dicom files to csv files you can find here:\n\n[https://www.kaggle.com/datasets/ibombonato/unifesp-xray-body-part-dicom-metadata-csv](http://)","metadata":{}},{"cell_type":"code","source":"train_dcm2csv_df = pd.read_csv('../input/unifesp-xray-body-part-dicom-metadata-csv/dicom_metadata_train.csv')\ntrain_dcm2csv_df","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:28.142200Z","iopub.execute_input":"2022-05-09T16:05:28.142451Z","iopub.status.idle":"2022-05-09T16:05:28.223408Z","shell.execute_reply.started":"2022-05-09T16:05:28.142415Z","shell.execute_reply":"2022-05-09T16:05:28.222703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Furthermore, there's another csv file in the main dataset that includes \"SOP IDs\" and \"targets\".","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv('../input/unifesp-x-ray-body-part-classifier/train.csv')\ntrain_df","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:28.224718Z","iopub.execute_input":"2022-05-09T16:05:28.225640Z","iopub.status.idle":"2022-05-09T16:05:28.251522Z","shell.execute_reply.started":"2022-05-09T16:05:28.225598Z","shell.execute_reply":"2022-05-09T16:05:28.250768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.info()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:28.253772Z","iopub.execute_input":"2022-05-09T16:05:28.254048Z","iopub.status.idle":"2022-05-09T16:05:28.274031Z","shell.execute_reply.started":"2022-05-09T16:05:28.254010Z","shell.execute_reply":"2022-05-09T16:05:28.273194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df['Target'].unique()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:28.275161Z","iopub.execute_input":"2022-05-09T16:05:28.275876Z","iopub.status.idle":"2022-05-09T16:05:28.284225Z","shell.execute_reply.started":"2022-05-09T16:05:28.275837Z","shell.execute_reply":"2022-05-09T16:05:28.283362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The labels are represented as integers that map to the following:\n\n\n* Abdomen = 0\n* Ankle = 1\n* Cervical Spine = 2\n* Chest = 3\n* Clavicles = 4\n* Elbow = 5\n* Feet = 6\n* Finger = 7\n* Forearm = 8\n* Hand = 9\n* Hip = 10\n* Knee = 11\n* Lower Leg = 12\n* Lumbar Spine = 13\n* Others = 14\n* Pelvis = 15\n* Shoulder = 16\n* Sinus = 17\n* Skull = 18\n* Thigh = 19\n* Thoracic Spine = 20\n* Wrist = 21","metadata":{}},{"cell_type":"markdown","source":"As you see, some of dicom files have more that one class (target)! So, we have a multilabel classification problem. This is my strategy:\n\nFirst we add all labels to the \"train_df\" as new columns by their orders and set all to zero. Then by considering the numbers in the target column of each record, we replace 1 in the columns related that number. \n\nFor example, suppose that target have two numbers: \"1 , 6\". These numbers represent \"ankle, feet\". We replace 1 in \"ankle\" and \"feet\" columns of that record.","metadata":{}},{"cell_type":"code","source":"train_df['Abdomen'] = 0\ntrain_df['Ankle'] = 0\ntrain_df['Cervical Spine'] = 0\ntrain_df['Chest'] = 0\ntrain_df['Clavicles'] = 0\ntrain_df['Elbow'] = 0\ntrain_df['Feet'] = 0\ntrain_df['Finger'] = 0\ntrain_df['Forearm'] = 0\ntrain_df['Hand'] = 0\ntrain_df['Hip'] = 0\ntrain_df['Knee'] = 0\ntrain_df['Lower Leg'] = 0\ntrain_df['Lumbar Spine'] = 0\ntrain_df['Others'] = 0\ntrain_df['Pelvis'] = 0\ntrain_df['Shoulder'] = 0\ntrain_df['Sinus'] = 0\ntrain_df['Skull'] = 0\ntrain_df['Thigh'] = 0\ntrain_df['Thoracic Spine'] = 0\ntrain_df['Wrist'] = 0\ntrain_df","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:28.285683Z","iopub.execute_input":"2022-05-09T16:05:28.285984Z","iopub.status.idle":"2022-05-09T16:05:28.353365Z","shell.execute_reply.started":"2022-05-09T16:05:28.285950Z","shell.execute_reply":"2022-05-09T16:05:28.352229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labels = ['Abdomen', 'Ankle', 'Cervical Spine', 'Chest', 'Clavicles', 'Elbow', 'Feet', 'Finger', 'Forearm', 'Hand', 'Hip', 'Knee', 'Lower Leg', 'Lumbar Spine', 'Others', 'Pelvis', 'Shoulder', 'Sinus', 'Skull', 'Thigh', 'Thoracic Spine', 'Wrist']\n\nfor i in range(len(train_df)):\n    lbl_list = train_df.Target[i].split()\n    \n    for j in lbl_list:\n        train_df.loc[i, labels[int(j)]] = 1\n        \n\ntrain_df[:10]","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:28.354891Z","iopub.execute_input":"2022-05-09T16:05:28.355182Z","iopub.status.idle":"2022-05-09T16:05:28.989141Z","shell.execute_reply.started":"2022-05-09T16:05:28.355146Z","shell.execute_reply":"2022-05-09T16:05:28.988460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we merge two dataframs based on \"SOPInstanceUID\" columns. So we have an integrated dataframe which help us to create \"X\" and \"y\" of model much easier.","metadata":{}},{"cell_type":"code","source":"train_merged_df = pd.merge(train_dcm2csv_df, train_df, on='SOPInstanceUID')\ntrain_merged_df","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:28.990318Z","iopub.execute_input":"2022-05-09T16:05:28.991070Z","iopub.status.idle":"2022-05-09T16:05:29.033242Z","shell.execute_reply.started":"2022-05-09T16:05:28.991029Z","shell.execute_reply":"2022-05-09T16:05:29.032467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"An example of converting dicom file to jpg image:","metadata":{}},{"cell_type":"code","source":"tmp = train_dcm_path[np.random.randint(0, 1737)].dcmread()\ndcm_img = tmp.pixel_array.astype(float)\nimg = np.uint8((np.maximum(dcm_img,0)/dcm_img.max())*255)\nplt.imshow(img, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:29.035193Z","iopub.execute_input":"2022-05-09T16:05:29.036592Z","iopub.status.idle":"2022-05-09T16:05:31.534644Z","shell.execute_reply.started":"2022-05-09T16:05:29.036564Z","shell.execute_reply":"2022-05-09T16:05:31.533937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This step is the most time consuming part of the code; \n\nFirst we should reading dicom files and convert them to jpg images but there is a big problem; The size of dicom images are too big; for example 3000 * 4000 pixels! There isn't enough memory so we should resize all of them to 128 * 128 pixels but unfortunately it will have consequences like information loss. If you have enough memory, you can resize them to 512 * 512 pixels for example. Anyway, the images create X array for training.\n\nThe \"y\" of model, is an array of 1738 * 22. In another words, the label of each image (each SOPInstanceUID) is a vector contains 22 binary numbers that represent 22 labels (body parts). ","metadata":{}},{"cell_type":"code","source":"import cv2\nfrom matplotlib import pyplot as plt\n\nX = []\ny = []\nfor i in range(len(train_dcm_path)):\n    tmp = train_dcm_path[i].dcmread()\n    \n    dcm_img = tmp.pixel_array.astype(float)\n    rescaled_img = np.uint8((np.maximum(dcm_img,0)/dcm_img.max())*255)\n    image = cv2.resize(rescaled_img, dsize=(128,128), interpolation=cv2.INTER_AREA)/255 #resizing and normalizing\n    X.append(np.expand_dims(image,axis=-1))\n    \n    uid = str(tmp['SOPInstanceUID'].value)\n    y.append(np.ndarray.flatten(np.array(train_merged_df.loc[train_merged_df['SOPInstanceUID'] == uid][train_merged_df.columns[64:]])))\n    \nX = np.array(X)\ny = np.array(y)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T16:05:31.536026Z","iopub.execute_input":"2022-05-09T16:05:31.536317Z","iopub.status.idle":"2022-05-09T18:17:22.388007Z","shell.execute_reply.started":"2022-05-09T16:05:31.536281Z","shell.execute_reply":"2022-05-09T18:17:22.387235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training and testing custom CNN model","metadata":{}},{"cell_type":"markdown","source":"First we should split X and y data to train and test data.","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.15, random_state=42)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:40:57.592735Z","iopub.execute_input":"2022-05-09T18:40:57.592999Z","iopub.status.idle":"2022-05-09T18:40:57.686242Z","shell.execute_reply.started":"2022-05-09T18:40:57.592969Z","shell.execute_reply":"2022-05-09T18:40:57.685453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To avoid of overfitting, we can use \"EarlyStopping\" callback; So we need to split train data to \"train\" and \"validation\" parts.","metadata":{}},{"cell_type":"code","source":"X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size=0.1, random_state=42)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:41:01.042559Z","iopub.execute_input":"2022-05-09T18:41:01.042882Z","iopub.status.idle":"2022-05-09T18:41:01.105951Z","shell.execute_reply.started":"2022-05-09T18:41:01.042846Z","shell.execute_reply":"2022-05-09T18:41:01.105201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### A sample of X_train and related y_train ###\n\nplt.imshow(X_train[7], cmap='gray')\nprint(y_train[7])","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:41:02.319597Z","iopub.execute_input":"2022-05-09T18:41:02.319855Z","iopub.status.idle":"2022-05-09T18:41:02.538054Z","shell.execute_reply.started":"2022-05-09T18:41:02.319827Z","shell.execute_reply":"2022-05-09T18:41:02.537210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we can build our CNN model.","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras.models import Model, Sequential\nfrom tensorflow.keras.layers import InputLayer, Conv2D, Dropout, MaxPooling2D, Flatten, Dense\n\ncnn_model = Sequential()\ncnn_model.add(InputLayer(input_shape=(128,128,1)))\n              \ncnn_model.add(Conv2D(64, (3,3), activation='relu', kernel_initializer='he_normal', padding='same'))\n#cnn_model.add(Conv2D(64, (3, 3), activation='relu', kernel_initializer='he_normal', padding='same'))\ncnn_model.add(MaxPooling2D((2, 2)))\n\ncnn_model.add(Conv2D(128, (3, 3), activation='relu', kernel_initializer='he_normal', padding='same'))\n#cnn_model.add(Conv2D(128, (3, 3), activation='relu', kernel_initializer='he_normal', padding='same'))\ncnn_model.add(MaxPooling2D((2, 2)))\n\ncnn_model.add(Conv2D(256, (3, 3), activation='relu', kernel_initializer='he_normal', padding='same'))\ncnn_model.add(Conv2D(256, (3, 3), activation='relu', kernel_initializer='he_normal', padding='same'))\ncnn_model.add(MaxPooling2D((2, 2)))\n\ncnn_model.add(Conv2D(512, (3, 3), activation='relu', kernel_initializer='he_normal', padding='same'))\ncnn_model.add(MaxPooling2D((2, 2)))\n\ncnn_model.add(Flatten())\ncnn_model.add(Dense(4096, activation=\"relu\"))\ncnn_model.add(Dense(2048, activation=\"relu\"))\ncnn_model.add(Dense(22, activation=\"sigmoid\"))\n\ncnn_model.compile(optimizer='adam', loss='BinaryCrossentropy', metrics=['accuracy'])\ncnn_model.summary()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:44:05.894213Z","iopub.execute_input":"2022-05-09T18:44:05.894620Z","iopub.status.idle":"2022-05-09T18:44:05.994129Z","shell.execute_reply.started":"2022-05-09T18:44:05.894576Z","shell.execute_reply":"2022-05-09T18:44:05.993444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = cnn_model.fit(x = X_train,\n                        y = y_train,\n                        batch_size = 8,\n                        callbacks = [tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=5)],\n                        validation_data = (X_val, y_val),\n                        verbose = 1,\n                        epochs = 50)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:44:11.813264Z","iopub.execute_input":"2022-05-09T18:44:11.813549Z","iopub.status.idle":"2022-05-09T18:45:01.327114Z","shell.execute_reply.started":"2022-05-09T18:44:11.813512Z","shell.execute_reply":"2022-05-09T18:45:01.326437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting results of trained model involvs accuracy, validation accuracy, loss and validation loss","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\nfig, ax = plt.subplots(nrows=1, ncols=2, figsize=(13,5))\n\nax[0].plot(history.history['accuracy'])\nax[0].plot(history.history['val_accuracy'])\nax[0].set_title('model accuracy')\nax[0].set_ylabel('accuracy')\nax[0].set_xlabel('epoch')\nax[0].legend(['train', 'val'], loc='upper left')\n\nax[1].plot(history.history['loss'][1:])\nax[1].plot(history.history['val_loss'])\nax[1].set_title('model loss')\nax[1].set_ylabel('val_loss')\nax[1].set_xlabel('epoch')\nax[1].legend(['train', 'val'], loc='upper right')\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:45:10.115394Z","iopub.execute_input":"2022-05-09T18:45:10.116139Z","iopub.status.idle":"2022-05-09T18:45:34.789012Z","shell.execute_reply.started":"2022-05-09T18:45:10.116095Z","shell.execute_reply":"2022-05-09T18:45:34.788336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The accuracy of the model on test data:","metadata":{}},{"cell_type":"code","source":"cnn_model.evaluate(x = X_test,\n                   y = y_test,\n                   verbose=1)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:45:42.376855Z","iopub.execute_input":"2022-05-09T18:45:42.377910Z","iopub.status.idle":"2022-05-09T18:45:42.621231Z","shell.execute_reply.started":"2022-05-09T18:45:42.377855Z","shell.execute_reply":"2022-05-09T18:45:42.620255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction\n\nFirst select a dicom file from test data and convert it to jpg image same as training phase. Then pass it to \"predict\" function and get output.\n\nAfter some tests, I set threshold = 0.1 which means if score is greater than 0.1, the related label is in that test image. ","metadata":{}},{"cell_type":"code","source":"test_candid = test_dcm_path[np.random.randint(0, 742)].dcmread()\ndcm_img = test_candid.pixel_array.astype(float)\nconverted_img = np.uint8((np.maximum(dcm_img,0)/dcm_img.max())*255)\nresized_img = cv2.resize(converted_img, dsize=(128,128), interpolation=cv2.INTER_AREA)/255\ntest_img = np.expand_dims(resized_img, axis=-1)\nplt.imshow(test_img, cmap='gray')\n\nimage = np.expand_dims(test_img, axis=0)\nprdct_arr = cnn_model.predict(image)\n\nbody_parts = []\nprdct_arr = prdct_arr.flatten()\nfor i in range(len(prdct_arr)):\n    if (prdct_arr[i] > 0.1):\n        body_parts.append(labels[i])\n\nprint(body_parts)","metadata":{"execution":{"iopub.status.busy":"2022-05-09T18:45:58.161975Z","iopub.execute_input":"2022-05-09T18:45:58.162237Z","iopub.status.idle":"2022-05-09T18:46:05.141217Z","shell.execute_reply.started":"2022-05-09T18:45:58.162207Z","shell.execute_reply":"2022-05-09T18:46:05.140539Z"},"trusted":true},"execution_count":null,"outputs":[]}]}