{"cells":[{"metadata":{},"cell_type":"markdown","source":"> **Diagnosing Diabetic Retinopathy from Fundus Photos **\n\nCatherine Ko; Computational Modeling and Simulation; March 6, 2020\n\n**About Diabetic Retinopathy**\n\nDiabetic Retinopathy is a diabetic eye disease that can lead to vision loss and blindness in extreme cases. It occurs when high blood sugar damages blood vessels in the retina (back of the eye). Diabetic retinopathy is typically diagnosed by ophthalmologists via fundus photography. They take a picture of your retina and look for the telltale signs as shown below:\n\n![](https://www.researchgate.net/profile/Asiri_Wijesinghe/publication/303481072/figure/fig1/AS:394097530556416@1470971581841/Retinal-lesions-in-DR-such-as-microaneurysms-exudates-and-hemorrhages-regions-13.png)\n\nMicroaneyrsm - a tiny swelling in the side of a blood vessel characterized by small red dots on the retina\n\nHemmorages - bleeding from a ruptured blood vessel\n\nSoft exudates (cotton-wool spots) - internal superficial leakage from retinal arteries\n\nHard exudates - extracellular lipid leaked from retinal capillaries characterized by yellow grains\n\n\n"},{"metadata":{},"cell_type":"markdown","source":"From a Kaggle competition, I was able to access ~3,600 retina images. \n\nEach image in the training dataset was labeled with one of the following diagnoses: \n\n0 - No DR\n\n1 - Mild\n\n2 - Moderate\n\n3 - Severe\n\n4 - Proliferative DR\n\nMy goal for this project was to diagnosis whether a patient has diabetic retinopathy (DR) (and to what degree) from running these fundus photos through my machine learning algorithm.\n\nFirst, I imported the data from the Kaggle competition by reading in the training csv."},{"metadata":{"_cell_guid":"","_uuid":"","trusted":true},"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\n#from keras.preprocessing import image\nimport cv2\nfrom tqdm import tqdm \nfrom PIL import Image, ImageEnhance","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"N=3662 # number of images \n\ndf_train = pd.read_csv(\"../input/aptos2019-blindness-detection/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at a portion of the training data csv:"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As shown above, the csv gives the image's id_code in the first column and the corresponding diagnosis for severity of diabetic retinopathy in the second column."},{"metadata":{},"cell_type":"markdown","source":"Now let's take a look at the distribution of diabetic retinopathy diagnosis by plotting a histogram of the counts and printing the raw counts:"},{"metadata":{"trusted":true},"cell_type":"code","source":"(counts,bin_edges,_)=plt.hist(df_train['diagnosis'],align='left',bins=5,edgecolor='purple',linewidth=1.5,color='lavender')\nplt.title('Histogram of diabetic retinopathy diagnosis', fontsize=14,weight='bold')\nplt.xticks(bin_edges[:-1], np.arange(0,5,1))\n\nax=plt.gca()\nax.set_xlabel('Diabetic Retinopathy Diagnosis',fontsize=12);\nax.set_ylabel('Frequency',fontsize=12);\n\ndf_train['diagnosis'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at the actual retina images we got from the Kaggle competition:"},{"metadata":{"trusted":true},"cell_type":"code","source":"paths='../input/aptos2019-blindness-detection/train_images/'+df_train['id_code'][:5]+'.png'\nplt.figure(figsize=(20,4))\nfor index, (path,label) in enumerate(zip(paths,df_train['diagnosis'][:5])):\n    plt.subplot(1,5,index+1)\n    plt.imshow(Image.open(path))\n    plt.title('Diagnosis: %i\\n' % label, fontsize = 20)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Image Preprocssing **\n\nBefore running the images through my machine learning algorithm, I wanted to standardize all the images and improve lighting conditions to make the telltale signs of diabetic retinopathy that we discussed earlier more obvious. \n\nSince these images come from a variety of clinics and imaging conditions, the sizes of the image are different and so I created a function called preprocess_image to resize the images to 224x224 and increase the image contrast by 1.5x. I chose to increase the contrast to make the telltale signs of diabetic retinopathy more obvious. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess_image(path, desired_size=224):\n    '''\n    resize image to desired size x desired size\n    and also increase contrast by 1.5x\n    '''\n    im = Image.open(path)\n    im = im.resize((desired_size, )*2,resample=Image.LANCZOS)\n    \n    # increase contrast of the images\n    enhancer = ImageEnhance.Brightness(im)\n    factor = 1.5 #factor > 1 increases the contrast\n    im_output = enhancer.enhance(factor)\n    return im_output","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After creating the function in charge of preprocessing the images, I ran all 3,662 images through my preprocessing function and stored the image data in a 4D array called pics_data. "},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time \n\n# reading in images into the array pics_data\n# https://www.kaggle.com/xhlulu/aptos-2019-densenet-keras-starter\n\n# create an empty 4d array to store the images\npics_data = np.empty((N, 224, 224,3), dtype=np.uint8)\n\nfor i, image_name in enumerate(tqdm(df_train['id_code'][:N])):\n    pics_data[i, :, :,:] = preprocess_image(f'../input/aptos2019-blindness-detection/train_images/{image_name}.png')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After processing my images, here's a glimpse of the images after preprocessing. Notice how the contrast is increased and the images are all the same size now. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,4))\nfor index, (path,label) in enumerate(zip(paths,df_train['diagnosis'][:5])):\n    plt.subplot(1,5,index+1)\n    plt.imshow(preprocess_image(path))\n    plt.title('Diagnosis: %i\\n' % label, fontsize = 20)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here I'm just reformatting the shape of my arrays to make sure I can run it through the next steps."},{"metadata":{"trusted":true},"cell_type":"code","source":"# store corresponding diagnoses in an array named y_data\ny_data=np.array(df_train['diagnosis'][:N])\n\n# print the shape of y_data array\ny_data.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# print the shape of pics_data\npics_data.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# reshape y_data from 1D --> 2D \ny_data2D=y_data.reshape(N,1)\n\n# reshape pics_data from 1D --> 2D \npics_data2D=pics_data.reshape(N,150528)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the next part, I split up the training images into a training set of images and a testing set of images. I also split up the corresponding diagnoses into a training set and testing set. We will use the training sets to train the model and use the testing sets to test the model and see how the model performs on data it hasn't seen before after training. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# split dataset into training and test sets in a way that is blind to the programmer\nfrom sklearn.model_selection import train_test_split\nx_train, x_test, y_train, y_test = train_test_split(pics_data2D, y_data2D, test_size=0.25, random_state=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, it is time to make the machine learning model! Below, I created a logistic regressor and trained it using my training datasets."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nfrom sklearn.linear_model import LogisticRegression\n# create logistic regressor\nlogisticRegr = LogisticRegression(random_state=0)\n# train logistic regressor using training sets\nlogisticRegr.fit(x_train, y_train.ravel()) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now the logistic regressor has been trained! \nTime to use the newly-trained model to make diagnoses predictions on our testing set of images."},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = logisticRegr.predict(x_test) # predict entire test set\ny_test_reshaped=y_test.reshape(916,)\ndf = pd.DataFrame({'Actual':y_test_reshaped,'Predicted':predictions})\ndf.tail()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at how accurate our model was. Below I printed some stats that will help us evaluate the error of this model."},{"metadata":{"trusted":true},"cell_type":"code","source":"# evaluating the algorithm\n\nfrom sklearn.metrics import classification_report, confusion_matrix\nprint(confusion_matrix(y_test,predictions))\nprint(classification_report(y_test,predictions))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Accuracy of the logistic regression model above is ~72% as shown in the number in the f1-score column and accuracy row. The f1-score for diagnosis 0 (no DR) was the highest at 93%, then diagnosis of 2 (moderate DR) at 61%, diagnosis of 1 (mild DR) at 45%, and diagnosis of 4 (proliferative DR) at 23%, and finally 3 (severe DR) at 14%. So in general, the model was better at correctly predicting you if did not have DR. Since the model's f1-score was not that high for people with severe and proliferative DR, the model would be likely to predict that you have less severe DR than you might actually have (more false negatives in a sense).\n\nThis would be less ideal than having a model with more false positives because this means that people with severe/proliferative DR are more likely to not be concerned enough with the actual state of their eye health."},{"metadata":{},"cell_type":"markdown","source":"**How exactly is the model being trained?**\n\nIn logistic regression (a form of supervised learning), when we are training the dataset, the program takes a whole bunch of x vectors from x_train set and runs them through the model $\\theta$. Then the program compares the calculated hypothesis (predicted) values h(x) and compares them to the actual y values stored in y_train. Mathematically, the computer calculates the cost function which is the sum of all the squared differences between actual - predicted y. \n\nThen the algorithm uses gradient descent and partial derivative calculus stuff behind the scenes to change the numbers in the model vector $\\theta$ until the cost function (error) is minimized! \n\nI'm using logistic regression as opposed to linear regression (another form of supervised learning) since I want the output to be categorical (either 0,1,2,3,4 depending on the severity of DR). \n"},{"metadata":{},"cell_type":"markdown","source":"**Future Directions**\n\nSome future directions for this project could be to improve the accuracy of my logistic regression and possibly create a neural network. I tried some other logistic regressors (with different number of training cycles and solver methods) and neural nets below but didn't have the time to bring those models to fruition."},{"metadata":{},"cell_type":"markdown","source":"** Some other logistic regressors **"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nfrom sklearn.linear_model import LogisticRegression\n# all parameters not specified are set to their defaults\nlogisticRegr = LogisticRegression(random_state=0,max_iter=1000)\nlogisticRegr.fit(x_train, y_train.ravel()) # train regressor \n# evaluating the algorithm\npredictions = logisticRegr.predict(x_test)\nfrom sklearn.metrics import classification_report, confusion_matrix\nprint(confusion_matrix(y_test,predictions))\nprint(classification_report(y_test,predictions))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nfrom sklearn.linear_model import LogisticRegression\n# all parameters not specified are set to their defaults\nlogisticRegr = LogisticRegression(random_state=0,solver='saga')\nlogisticRegr.fit(x_train, y_train.ravel()) # train regressor \n# evaluating the algorithm\npredictions = logisticRegr.predict(x_test)\nfrom sklearn.metrics import classification_report, confusion_matrix\nprint(confusion_matrix(y_test,predictions))\nprint(classification_report(y_test,predictions))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"** Trying out Neural Nets **"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nfrom sklearn.neural_network import MLPClassifier\nmlp=MLPClassifier(hidden_layer_sizes=(10,10,10),max_iter=20,random_state=0) # creates neural network\n# hidden_layer_sizes creates 3 layers of 10 nodes each; just try different combos and see what is best\n#max_iter = number of iterations of epochs (cycles of feed-forward and back propagation)\nmlp.fit(x_train, y_train.ravel())\n# make predictions to our test data\npredictions=mlp.predict(x_test)\n# evaluating the neural net algorithm\n\nfrom sklearn.metrics import classification_report, confusion_matrix\nprint(confusion_matrix(y_test,predictions))\nprint(classification_report(y_test,predictions))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nfrom sklearn.neural_network import MLPClassifier\nmlp=MLPClassifier(hidden_layer_sizes=(10,10,10),random_state=0) # creates neural network\n# hidden_layer_sizes creates 3 layers of 10 nodes each; just try different combos and see what is best\n#max_iter = number of iterations of epochs (cycles of feed-forward and back propagation)\nmlp.fit(x_train, y_train.ravel())\n# make predictions to our test data\npredictions=mlp.predict(x_test)\n# evaluating the neural net algorithm\n\nfrom sklearn.metrics import classification_report, confusion_matrix\nprint(confusion_matrix(y_test,predictions))\nprint(classification_report(y_test,predictions))","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}