{"cells":[{"metadata":{},"cell_type":"markdown","source":"## RSNA: STR Plumonary Embolism Detection\n\n\n### What is plumonary Embolism ?\n\n\n\nPlumonary Embolism is the blood clot that is formed in the artery of the lungs.
\nIt is a very serious disease which causes death of one-third people who go undiagnosed or undetected\nwith the disease.
\nBut immediate treatment could help prevent permenant lung damage.\n\n### What causes Plumonary Embolism ?\n\nPlumonary Embolism is caused due to the blood clot formed deep in the vein of our body.
\nMain reason for the blood clot is [deep vein thrombosis](https://www.healthline.com/health/deep-venous-thrombosis)
The blood clot which causes plumonary embolism are formed in legs or pelvis.\n\n### Reason for blood clots.\n\n**Injury Damage**: bone fracture , tissue damage etc.
\n**Inactivity**: sitting for long period of time or lyinig on bed because of illness.
\n**Medical Condition**: There are some medical condition which can cause blood clot to form in our body.
\n\n### symptoms of plumonary embolism.\n\nThe symptoms of plumonary embolism depends on the size of blood clot.
\nMain symptoms is shortness of breath which can be gradual or sudden.\nother symptoms.\n\n* anxiety\n* clammy or bluish skin\n* chest pain that may extend into your arm, jaw, neck, and shoulder\n* fainting\n* irregular heartbeat\n* lightheadedness\n* rapid breathing\n* rapid heartbeat\n* restlessness\n* spitting up blood\n* weak pulse\n\nRead more about Plumonary Embolism [here](https://www.healthline.com/health/pulmonary-embolus)."},{"metadata":{},"cell_type":"markdown","source":"## 1. What is the competition about ?π‘\n\nHere we are provided with the CT Plumonary angiography(CTPA) which is common type of medical imaging.
\nData Contains hundereds of such images and we have to use our model to predict wether the image has PE or Not.\n\n**Let me know if any information or code is incorrect i will correct it and
\n if you find this notebook usefull please UPVOTE π**\n"},{"metadata":{},"cell_type":"markdown","source":"## 2.Metrics: Weighted Log-loss π\n\nThe metrics used by the competion is weighte log-loss which is weighted over some labels.\n\nWe have to predict total 10 labels 9 for exam/study level and 1 for image level.\n\nso submission file should have number of rows equal to
\n(number of images) + (number of exam/study label * number of exam/study)\n\nlabels for exam/study level.\n\n* Label: Weight\n* Negative for PE : 0.0736196319\n* Indeterminate \t:0.09202453988\n* Chronic \t:0.1042944785\n* Acute & Chronic \t:0.1042944785\n* Central PE \t:0.1877300613\n* Left PE \t:0.06257668712\n* Right PE \t:0.06257668712\n* RV/LV Ratio >= 1 :0.2346625767\n* RV/LV Ratio < 1 :0.0782208589\n\n### 2.1 Exam/Study- level weighted log-loss\n\nyij be the label for the exam i and label j. yij = 1 if present else 0 and pij
\nis the predicted probability. weight of the label j be wj. so the weighted log-loss is given by.\n\nexam_log_loss = -wj * [yij * log(pij) + (1 - yij)* log(1-pij)]\n\nThen mean is taken over the log_loss for all such labels j.\n\n### 2.2 Image level weighted log-loss\n\nyik =1 if PE is present in the image else 0. where i is exam number and k is image number.
\nNow the weigtage of the label PE present or not is w = 0.0736196319(\"Negative for PE\").
\nqi be the ratio of positive images to total images.\n\nimage_log_loss = -w * qi [yik * log(pik) + (1-yik)*log(1-pik)]\n\nThe total loss is the average of all image and exam loss, divided by the average of all row weights.
\nTo get the average of all rows weights, sum the weights of all images and all exam-level labels and divide by number of rows.\n\nNow if that was not confusing enough for you there is another catch.\n\n[evaluation](https://www.kaggle.com/c/rsna-str-pulmonary-embolism-detection/overview/evaluation) page.\n\n### 2.3 All the labels must be logically consistent.\n\nAll the labels in the submission file must be logically consistent or your submission will be disqualified.\n\nWhat is meant by logically consistent.\n\nAt the image level, any image with predicted probability > 0.5 is considered as being positive for PE will count as a positive image\n\nAt the exam level, we have\n\n1. Negative, Indeterminate, (Positive) and it can only be one of these. If any image is predicted positive, there cannot also be a predicted probability of Negative > 0.5 nor can there be a predicted probability of Indeterminate > 0.5.\n\nSimilarly, if no image is positive (p > 0.5), then there must be one and only one negative or indeterminate with p > 0.5\n\n1. Right, left, central -- if any image is predicted positive (p > 0.5) then at least one of these labels must be assigned p > 0.5; more than one of these labels may be assigned p > 0.5. When no images are predicted positive, then none of these labels may be assigned p > 0.5\n\n2. RV/LV ratio. It can be only one of these and it must be present if at least one image is positive.\n\n * if any image on the exam is positive, one of these must have p > 0.5 both cannot have p > 0.5\n\n3. Acute, Chronic, Acute & Chronic -- it cannot be both chronic & acute and chronic so\n\n * only one can have p > 0.5\n * it is also possible that neither has p > 0.5\n * in other words, it is inconsistent to say chronic has p > 0.5 and acute & chronic has p > 0.5.\n\nCode for checking consitency of the submission is provided [here](https://www.kaggle.com/anthracene/host-confirmed-label-consistency-check)
\nwhich will be used in this notebook"},{"metadata":{},"cell_type":"markdown","source":"### Importing Libraries π"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import os\nimport sys\nimport glob\nimport tqdm\nfrom typing import Dict\nimport cv2\nfrom collections import Counter\nimport random\n\nimport pydicom as dicom\nfrom joblib import Parallel, delayed\n\n\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport matplotlib.animation as animation\nfrom IPython.display import HTML\n\n\nplt.style.use(\"fivethirtyeight\")\n\nimport plotly.express as px\nimport plotly.graph_objs as go\nfrom plotly.offline import iplot\nimport plotly.figure_factory as ff\n\n#supress warnings\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfrom colorama import Fore, Back, Style\ny_ = Fore.YELLOW\nr_ = Fore.RED\ng_ = Fore.GREEN\nb_ = Fore.BLUE\nm_ = Fore.MAGENTA\nsr_ = Style.RESET_ALL","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Getting dataπ½"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"folder_path = \"../input/rsna-str-pulmonary-embolism-detection\"\ntrain_path = folder_path + \"/train/\"\ntest_path = folder_path + \"/test/\"\n \ntrain_data = pd.read_csv(folder_path + \"/train.csv\")\ntest_data = pd.read_csv(folder_path + \"/test.csv\")\nsample = pd.read_csv(folder_path + \"/sample_submission.csv\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(f\"{y_}Number of rows in train data: {r_}{train_data.shape[0]}\\n{y_}Number of columns in train data: {r_}{train_data.shape[1]}\")\nprint(f\"{g_}Number of rows in test data: {r_}{test_data.shape[0]}\\n{g_}Number of columns in test data: {r_}{test_data.shape[1]}\")\nprint(f\"{b_}Number of rows in submission data: {r_}{sample.shape[0]}\\n{b_}Number of columns in submission data:{r_}{sample.shape[1]}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**There are four labels in training set which are just for information and requires no predictions OA Constrast, QA Motion, True filling defect not PE and Flow artifact**."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"columns_only_for_info = [\"qa_motion\",\"qa_contrast\",\"flow_artifact\",\"true_filling_defect_not_pe\"]\ncols_ID = [\"StudyInstanceUID\",\"SeriesInstanceUID\",\"SOPInstanceUID\"]\n\ndef highlight_cols(x):\n df = x.copy()\n df.loc[:,:] = 'background-color: lightgreen'\n df[columns_only_for_info] = 'background-color: red'\n return df \n\ntrain_data.head().style.apply(highlight_cols, axis=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Location of each image is as \"StudyInstanceUID\"/SeriesInstanceUID/SOPInstanceUID.dcm
\nlet us perform sanity check wether StudyInstanceUID , SeriesInstanceUID and SOP instance UID in data match with actual folders"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def sanity_check():\n #get all StudyInstanceUID\n train_StudyInstanceUID = os.listdir(train_path)\n test_StudyInstanceUID = os.listdir(test_path)\n \n train_SeriesInstanceUID = list()\n test_SeriesInstanceUID = list()\n \n #get all SeriesInstanceUID \n for x in train_StudyInstanceUID:\n train_SeriesInstanceUID.extend(os.listdir(train_path+x))\n \n \n for x in test_StudyInstanceUID:\n test_SeriesInstanceUID.extend(os.listdir(test_path+x))\n \n train_SOPInstanceUID = list()\n test_SOPInstanceUID = list()\n \n #get all SOPInstanceUID\n for x,y in zip(train_StudyInstanceUID,train_SeriesInstanceUID):\n train_SOPInstanceUID.extend(os.listdir(train_path+x+\"/\"+y))\n\n for x,y in zip(test_StudyInstanceUID,test_SeriesInstanceUID):\n test_SOPInstanceUID.extend(os.listdir(test_path+x+\"/\"+y))\n \n #removing the extention\n train_SOPInstanceUID = [x.split(\".\")[0] for x in train_SOPInstanceUID]\n test_SOPInstanceUID = [x.split(\".\")[0] for x in test_SOPInstanceUID]\n \n ## Note: there might be better way to do this\n train_UIDs = [train_StudyInstanceUID,train_SeriesInstanceUID,train_SOPInstanceUID]\n test_UIDs = [test_StudyInstanceUID,test_SeriesInstanceUID,test_SOPInstanceUID]\n cols = [\"StudyInstanceUID\",\"SeriesInstanceUID\",\"SOPInstanceUID\"]\n \n\n for y,x in zip(train_UIDs,cols):\n if train_data[x].nunique() == len(y):\n print(f\"{g_}Number of {x} in train data and folder are same\")\n\n if Counter(train_data[x].unique()) == Counter(y):\n print(f\"{g_}{x} in train csv and folder are same\")\n\n else:\n print(f\"{r_}{x} in train csv and folder are not same\")\n not_in_train_folder = y[Counter(train_data[x].unique()) != Counter(y)]\n print(\"Not in train folder: \",not_in_train_folder)\n\n else:\n print(f\"{r_}Number of {x} in train data and folder are same\")\n\n for y,x in zip(test_UIDs,cols):\n if test_data[x].nunique() == len(y):\n print(f\"{g_}Number of {x} in test data and folder are same\")\n\n if Counter(test_data[x].unique()) == Counter(y):\n print(f\"{g_}{x} in test csv and folder are same\")\n\n else:\n print(f\"{r_}{x} in test csv and folder are not same\")\n not_in_test_folder = y[Counter(test_data[x].unique()) != Counter(y)]\n print(\"Not in test folder\",not_in_test_folder)\n else:\n print(f\"{r_}Number of {x} in test data and folder are same\")\n ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sanity_check()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data[\"ImagePath\"] =train_path+ train_data[cols_ID[0]]+\"/\"+train_data[cols_ID[1]]+\"/\"+train_data[cols_ID[2]]+\".dcm\"\ntest_data[\"ImagePath\"] = test_path+ test_data[cols_ID[0]]+\"/\"+test_data[cols_ID[1]]+\"/\"+test_data[cols_ID[2]]+\".dcm\"\n\ntrain_data[\"ImageDict\"] =train_path+ train_data[cols_ID[0]]+\"/\"+train_data[cols_ID[1]]\ntest_data[\"ImageDict\"] = test_path+ test_data[cols_ID[0]]+\"/\"+test_data[cols_ID[1]]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"dicom.dcmread(train_data.loc[0,\"ImagePath\"])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Basic Physics Behind CT scan.\n\n* So Idea behind CT(computed tomography) scan is that when a x-ray beam is passed through a tissue\n [attenuation factor](https://radiopaedia.org/articles/attenuation-coefficient?lang=us) of tissue is calculated.\n\n* attenuation factor is a measure of by how much strength of the beam is reduced when it is passed through\n certain material. \n \n* As different body parts have different attenuation factor it helps to seprate those parts easily.\n\n### Hounsfield Unit.\n\n* In 1979 **[Godfrey Hounsfield](https://en.wikipedia.org/wiki/Godfrey_Hounsfield)** got Nobel Price for his part in making CT scanners\n\n* Task of CT scanner is to find density of different tissues in our body now to
\n convert this density to gray scale image Hounsfield came up with a way of linear transforming
\n the density to gray scale image.
\n\n* Hounsfield Unit is calculated by considering density of water as 0 unit and air as -1000 unit
\n and all other density are transformed accordingly. Normally body parts lies in the range
\n -1000(ex air in lungs) to +1000.\n \n* This -1000 to 1000 range is adjueste for 256 gray values to come up with image.\n\n* As the values spectral composition of the x-ray depands on various parameters and voltage Hounsfield unit
\n makes it easier to compare CT scans from different machines.\n\n### window and level\n* Gray scale has 256 different units to represent colours but our eyes are not able to see small changes in colour
\n in gray images. \n\n* so instead of spreading whole Hounsfield Range to a gray image we select window of Hounsfield
\n and spread 256 values between these Hounsfield values \n \n* everything below this range is black and above the range is white.
\n\n* Level is center of this window. \n\n* We can adjust this window size according to our need to get better view of particular tissue we
\n want to observer.\n\nwatch [this](https://www.youtube.com/watch?v=KZld-5W99cI&feature=youtu.be) video for more info\n\n### voxel size\n\n* voxel is basically 3d version of pixel as pixel is representation of image in 2d voxel is a 3d representation.\n here it means pixelspacing in x and y direction and slice thickness in z direction\n \n* There are two type of voxel isotropic in which step in all sides are same and non isotropic steps are different.\n\n* smaller voxel sizes give better image.\n\n\n\n\n\n"},{"metadata":{},"cell_type":"markdown","source":"### Metadata of Image"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"N = 10000\n\ndef get_value(data):\n if type(data) == dicom.multival.MultiValue:\n return np.int(data[0])\n else:\n return np.int(data)\n \ndef get_meta_features(path):\n data = dicom.dcmread(path)\n \n slicethickness = data.SliceThickness\n windowwidth = get_value(data.WindowWidth)\n rows = data.Rows\n columns = data.Columns\n windowcenter = get_value(data.WindowCenter)\n intercepts = data.RescaleIntercept\n slopes = data.RescaleSlope\n pixelspacingcolumn = data.PixelSpacing[1]\n pixelspacingrows = data.PixelSpacing[0]\n kvp = data.KVP\n tableheight = data.TableHeight\n xray = data.XRayTubeCurrent\n exposure = data.Exposure\n modality = data.Modality\n rotationdirection = data.RotationDirection\n instancenumber = data.InstanceNumber\n \n final_data = [slicethickness,windowwidth,rows,columns,windowcenter,intercepts,\n slopes, pixelspacingcolumn,pixelspacingrows,kvp,tableheight,\n xray,exposure,modality,rotationdirection,instancenumber]\n return final_data\n\nmeta_data = Parallel(n_jobs = -1, verbose = 1)(map(delayed(get_meta_features),train_data[\"ImagePath\"].sample(n=N)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_data = pd.DataFrame(meta_data,\n columns = [\"SliceThickness\",\n \"WindowWidth\",\n \"Rows\",\n \"Columns\",\n \"WindowCenter\",\n \"Intercept\",\n \"Slope\",\n \"PixelSpacingRows\",\n \"PixelSpacingColumns\",\n \"KVP\",\n \"TableHeight\",\n \"XRay\",\n \"Exposure\",\n \"Modality\",\n \"RotationDirection\",\n \"InstanceNumber\"])\n\nmeta_data[\"Area\"] = meta_data[\"Rows\"] * meta_data[\"Columns\"]\nmeta_data[\"PixelArea\"] = meta_data[\"PixelSpacingRows\"] * meta_data[\"PixelSpacingColumns\"]\nmeta_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3 EDA"},{"metadata":{},"cell_type":"markdown","source":"### 3.1 Check Row and Columns range"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def dist(column,color):\n sns.distplot(meta_data[column],label=column,color=color)\n plt.legend()\n\nplt.figure(figsize=(15,7))\nplt.subplot(121)\ndist(\"Rows\",\"blue\")\nplt.subplot(122)\ndist(\"Columns\",\"green\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are 512 rows and columns in all 10000 images so it might be safe to assume that all the images follows same."},{"metadata":{},"cell_type":"markdown","source":"### 3.2 Distribution of Pixel Spacing Columns and Rows"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,7))\nplt.subplot(121)\ndist(\"PixelSpacingRows\",\"purple\")\nplt.subplot(122)\ndist(\"PixelSpacingColumns\",\"red\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.3 Distribution of Pixel Area"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(dpi=100)\ndist(\"PixelArea\",\"yellow\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.4 Distribution of Window width"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(dpi=100)\ndist(\"WindowWidth\",\"orange\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.5 Distribution of KVP"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(dpi=100)\ndist(\"KVP\",\"blue\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.6 Distribution of TableHeight"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(dpi=100)\ndist(\"TableHeight\",\"brown\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.7 Distribution of XRay "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(dpi=100)\ndist(\"XRay\",\"pink\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3.8 countplot of modality and rotation direction\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.subplot(121)\nsns.countplot(meta_data['Modality'])\nplt.subplot(122)\nsns.countplot(meta_data[\"RotationDirection\"])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It means all the Modality are CT and rotation direction are CW"},{"metadata":{},"cell_type":"markdown","source":"### 3.9 Distribution of Image Raw and Rescaled values"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def distribution_of_image_values(n,train=True):\n samples = train_data.sample(n=n) if train else test_data.sample(n=n)\n image_paths = samples[\"ImagePath\"].values\n \n plt.figure(figsize=(15,7))\n \n for i,image_path in enumerate(image_paths):\n image_data = dicom.dcmread(image_path)\n try:\n image = image_data.pixel_array.flatten()\n rescaled_image = image * image_data.RescaleSlope + image_data.RescaleIntercept\n \n plt.subplot(121)\n sns.distplot(image.flatten())\n plt.title(\"Raw Image\")\n \n plt.subplot(122)\n sns.distplot(rescaled_image.flatten())\n plt.title(\"Rescaled Image\")\n except:\n pass\n plt.show() ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"distribution_of_image_values(100)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"import gc\ndel meta_data\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4. Visulizing Images πΌοΈ\n\n\n#### What is DICOM image ? \n\n**DICOM(Digital Image and Communication in Medecine)** is a standard developed and
\nmaintained by **National Electrical Manufacturers Association (NEMA)** for storing
\nand transfering the medical images like CT(computerised Tomography),
\nMagnetic resonanse image(MRI) and other types of medical images.\n\nDICOM is very good protocol and intresting, to read further click [here](https://en.wikipedia.org/wiki/DICOM)\n"},{"metadata":{},"cell_type":"markdown","source":"### 4.1 Single Image πΌοΈ"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def show_image(train=True):\n image_path = train_data[\"ImagePath\"].sample(n=1).values[0] if train\\\n else test_data[\"ImagePath\"].sample(n=1).values[0]\n print(f\"{y_} Image {r_}{image_path}\")\n image = dicom.dcmread(image_path)\n image = image.pixel_array\n plt.figure(figsize=(7,7))\n plt.imshow(image,cmap='gray')\n plt.axis('off')\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_image()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4.2 Grid of sorted images of some random patient "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def show_grid(cmap='gray',train=True):\n single_sample = train_data.sample(n=1) if train else test_data.sample(n=1)\n image_dict = single_sample[\"ImageDict\"].values[0]\n \n images = [dicom.read_file(image_dict+\"/\"+filename) for filename in os.listdir(image_dict)]\n images.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n plt.figure(figsize=(10,10))\n \n for i,image in enumerate(images[:100]):\n plt.subplot(10,10,i+1)\n plt.imshow(image.pixel_array,cmap=cmap)\n plt.axis('off')\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_grid()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_grid(cmap='jet',train=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_grid(cmap='RdYlBu')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4.3 Animation"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"def show_animation(train=True):\n single_sample = train_data.sample(n=1) if train else test_data.sample(n=1)\n image_dict = single_sample[\"ImageDict\"].values[0]\n\n images = [dicom.read_file(image_dict+\"/\"+filename) for filename in os.listdir(image_dict)]\n images.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n fig = plt.figure()\n ims = list()\n for image in images:\n img = plt.imshow(image.pixel_array,cmap='gray',animated=True)\n plt.axis('off')\n ims.append([img])\n ani = animation.ArtistAnimation(fig,ims,interval=100,blit=False,repeat_delay=1000)\n return ani\n\nani = show_animation() ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"HTML(ani.to_jshtml())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4.4 3d- Reconstruction"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_x_y_cordinate(image,n1,n2):\n arr = np.argwhere((image >= n1) & (image <= n2))\n return arr\n \ndef reconstruct(train=True,n1=-1000,n2=2000,s=2,color='b',alpha=0.01,number=1000):\n single_sample = train_data.sample(n=1) if train else test_data.sample(n=1)\n image_dict = single_sample[\"ImageDict\"].values[0]\n \n images = [dicom.read_file(image_dict+\"/\"+filename) for filename in os.listdir(image_dict)]\n images.sort(key = lambda x: float(x.ImagePositionPatient[2]))\n \n fig = plt.figure(figsize=(10,10))\n ax = fig.add_subplot(111,projection='3d')\n \n for i,image in enumerate(images[:number]):\n img = image.pixel_array\n arr = get_x_y_cordinate(img,n1,n2)\n x = arr[:,0]\n y = arr[:,1]\n z = np.full(shape=len(x),fill_value = i+images[i].SliceThickness)\n ax.scatter(x,y,z,s=s,c=color,alpha=alpha)\n \n ax.set_xlabel('X Label')\n ax.set_ylabel('Y Label')\n ax.set_zlabel('Z Label')\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reconstruct(n1=1500,n2=2000,alpha=0.1,s=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"well I though it would work"},{"metadata":{},"cell_type":"markdown","source":"## Pytorch Baseline Model π₯"},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.functional as F\nimport torch.optim as optim\nimport torchvision\nfrom torchvision import models\nfrom torch.utils.data import Dataset,DataLoader\nimport cv2\nimport albumentations as A\nfrom albumentations.pytorch.transforms import ToTensorV2\n\nfrom sklearn.model_selection import KFold\n\nimport vtk\nfrom vtk.util import numpy_support\nfrom tqdm.auto import tqdm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"SEED = 42\n\ndef seed_everything(seed):\n random.seed(seed)\n os.environ['PYTHONHASSEED'] = str(seed)\n np.random.seed(seed)\n torch.manual_seed(seed)\n torch.cuda.manual_seed(seed)\n torch.backends.cudnn.deterministic = True\n torch.backends.cudnn.benchmark = True\n\nseed_everything(SEED)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"target_columns = ['pe_present_on_image', 'negative_exam_for_pe', 'rv_lv_ratio_gte_1', \n 'rv_lv_ratio_lt_1','leftsided_pe', 'chronic_pe','rightsided_pe', \n 'acute_and_chronic_pe', 'central_pe', 'indeterminate']\n\n#vtk is used because dicom is giving some error\n\nreader = vtk.vtkDICOMImageReader()\ndef get_img(path):\n reader.SetFileName(path)\n reader.Update()\n _extent = reader.GetDataExtent()\n ConstPixelDims = [_extent[1]-_extent[0]+1, _extent[3]-_extent[2]+1, _extent[5]-_extent[4]+1]\n\n ConstPixelSpacing = reader.GetPixelSpacing()\n imageData = reader.GetOutput()\n pointData = imageData.GetPointData()\n arrayData = pointData.GetArray(0)\n ArrayDicom = numpy_support.vtk_to_numpy(arrayData)\n ArrayDicom = ArrayDicom.reshape(ConstPixelDims, order='F')\n ArrayDicom = cv2.resize(ArrayDicom,(512,512))\n return ArrayDicom\n\n\ndef convert_to_rgb(array):\n array = array.reshape((512, 512, 1))\n return np.stack([array, array, array], axis=2).reshape((3,512, 512))","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"class RsnaDataset(Dataset):\n \n def __init__(self,df,transforms=None):\n super().__init__()\n self.image_paths = df['ImagePath'].unique()\n self.df = df\n self.transforms = transforms\n \n def __getitem__(self,index):\n \n image_path = self.image_paths[index]\n data = self.df[self.df['ImagePath']==image_path]\n labels = data[target_columns].values.reshape(-1)\n image = get_img(image_path)\n image = convert_to_rgb(image)\n \n if self.transforms:\n image = self.transforms(image=image)['image']\n \n image = torch.tensor(image,dtype=torch.float) \n labels = torch.tensor(labels,dtype=torch.float)\n \n return image,labels\n \n def __len__(self):\n return self.image_paths.shape[0] ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"classes = len(target_columns)\nmodel = models.resnet18(pretrained=True)\nin_features = model.fc.in_features\nmodel.fc = nn.Linear(in_features,classes)\n\nconfig={\n \"learning_rate\":0.001,\n \"train_batch_size\":32,\n \"valid_batch_size\":32,\n \"epochs\":10,\n \"nfolds\":3,\n \"number_of_samples\":7000\n }\n\ntrain_data = train_data.sample(n=config[\"number_of_samples\"]).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def run(plot_losses=True):\n \n def train_loop(train_loader,model,loss_fn,device,optimizer,lr_scheduler=None):\n model.train()\n total_loss = 0\n tqdm_loader = tqdm(train_loader)\n for i, (images, targets) in enumerate(tqdm_loader):\n images,targets = images.to(device), targets.to(device)\n \n optimizer.zero_grad()\n outputs = model(images)\n\n loss = loss_fn(outputs,targets)\n loss.backward()\n \n total_loss += loss.item()\n\n optimizer.step()\n if lr_scheduler != None:\n lr_scheduler.step()\n \n total_loss /= len(train_loader)\n return total_loss\n \n def valid_loop(valid_loader,model,loss_fn,device):\n model.eval()\n total_loss = 0\n predictions = list()\n tqdm_loader = tqdm(valid_loader)\n\n for i, (images, targets) in enumerate(tqdm_loader):\n images, targets = images.to(device),targets.to(device)\n \n outputs = model(images) \n\n loss = loss_fn(outputs,targets)\n predictions.extend(outputs.detach().cpu().numpy())\n \n total_loss += loss.item()\n total_loss /= len(valid_loader)\n \n return total_loss,np.array(predictions) \n \n kfold = KFold(n_splits=config[\"nfolds\"])\n\n fold_train_losses = list()\n fold_valid_losses = list()\n \n device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n print(f\"{device} is used\")\n \n loss_fn = nn.BCEWithLogitsLoss()\n \n for k , (train_idx,valid_idx) in enumerate(kfold.split(train_data)):\n \n x_train,x_valid = train_data.loc[train_idx,:],train_data.loc[valid_idx,:]\n \n model.to(device)\n\n train_ds = RsnaDataset(x_train)\n train_dl = DataLoader(train_ds,\n batch_size = config[\"train_batch_size\"],\n shuffle=True\n )\n\n valid_ds = RsnaDataset(x_valid)\n valid_dl = DataLoader(valid_ds,\n batch_size = config[\"valid_batch_size\"],\n shuffle=False\n )\n \n optimizer = optim.Adam(model.parameters(),lr=config[\"learning_rate\"])\n lr_scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer,T_max= 300,eta_min= 0.000001)\n\n print(f\"Fold {k}\")\n best_loss = 999\n train_losses = list()\n valid_losses = list()\n \n for i in range(config['epochs']):\n train_loss = train_loop(train_dl,model,loss_fn,device,optimizer,lr_scheduler)\n valid_loss,predictions = valid_loop(valid_dl,model,loss_fn,device)\n \n train_losses.append(train_loss)\n valid_losses.append(valid_loss)\n \n print(f\"epoch:{i} Training | loss:{train_loss} Validation | loss:{valid_loss} \")\n \n if valid_loss <= best_loss:\n print(f\"{g_}Validation loss Decreased from {best_loss} to {valid_loss}{sr_}\")\n best_loss = valid_loss\n torch.save(model.state_dict(),f'model{k}.bin')\n \n fold_train_losses.append(train_losses)\n fold_valid_losses.append(valid_losses)\n \n if plot_losses == True:\n plt.figure(figsize=(20,14))\n for i, (t,v) in enumerate(zip(fold_train_losses,fold_valid_losses)):\n plt.subplot(2,5,i+1)\n plt.title(f\"Fold {i}\")\n plt.plot(t,label=\"train_loss\")\n plt.plot(v,label=\"valid_loss\")\n plt.legend()\n plt.show() ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"run()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}