--- title: "Diabetic retinopathy diagnosis from retinal texture descriptors" author: "Jesús Martín de la Sierra" output: html_document --- # Introduction This notebook is a tidy version of the draft solution submitted to the [APTOS 2019 Blindness Detection](https://www.kaggle.com/c/aptos2019-blindness-detection/) competition. Such solution was scored much lower than the expected from the validation set, but it's still interesting to share my approach since I didn't see any work based on texture patterns in this competition. At the same time, the solution shown here is only a part of a wider research that studied the automatic detection and removal of the optic disk and the retinal vessels to clean the images from interferer elements. The improvement with these operations wasn't really proven since the processing time per image was unacceptable for the competition scoring requirements. So, only a lightweight version of all this work was finally considered. The metric to score in this competition was the quadratic weighted kappa. While the validation set scored around 0.70 (in contrast to >0.80 in my original study), the private test score droped to ~0.55. I think that my approach didn't work properly because of the images in the test sets. I suspect the images in there show important differences in framing, focus and/or general quality, being the method sensitive to these changes if they are present in a significative fraction of the data set. My approach explores the idea of measuring the general texture for each image, so we can then use it for classification. The operation used was the Local Binary Pattern which, for each position in an image, calculates a binary sequence given by the neighbor pixels. This sequence is converted to decimal and a histogram of these values for the whole image is computed. Once we have the histogram for each image, then the solution basically reduces to a histogram classifier, task done ultimately by a neural network in my case. LBPs are widely used for different purposes like image segmentation, face recognition or surface diagnostic as an alternative, not so accurate but close in some cases, to other methods like those based on neural networks for image detection only. Specifically, this work is inspired by the research _Identify Diabetic Retinopathy by Color Textures (Vo & Verma, 2016)_ where the authors obtain an accuracy up to 71% using LBPs as a key part of the entire process. Here in my study, the accuracy is even greater for the validation test, but unkown for the test sets. # Packages ```{r message=FALSE, warning=FALSE} # Packages needed for this notebook library(imager) library(wvtool) library(keras) library(Metrics) ``` # Texture descriptors for training set The first step is to extract the texture descriptors, i.e. the LBPs histograms, for each image in the training data set. Before that, a reduction of the images is necessary. The green channel is extracted since it keeps a good contrast. A mask is defined to separate background from foreground, and a cropping and rescale is performed such that the greatest part of the retina in the original image is enclosed now in an area of 128x128 pixels. Then, two patterns are computed by image: one for radius 1 and another for radius 2. The radius indicates the neighbor distance from the pixel of interest. Both patterns are also computed in uniform mode, which uses a lower amount of levels because considers the number of transitions 0-1 and 1-0. This makes the texture descriptor more rotational invariant and the number of dimensions also decreases: 59 by descriptor, that is, 118 by image. Much lower than the 16384 dimensions if we consider each pixel to be processed by image. This part is not processed in this notebook because it's quite time consuming. Instead, this has been previously ran and the result saved to CSV files. Then, they are loaded in the next step to train the model. ```{r message=FALSE, warning=FALSE, eval=FALSE} # List all training files file_list <- dir(path="../input/aptos2019-blindness-detection/train_images/", pattern="*.png") # Patterns initialization for the complete 128x128 image complete_r1 <- matrix(nrow=0, ncol=59) complete_r2 <- matrix(nrow=0, ncol=59) for (i in 1:length(file_list)) { # Load image image <- load.image(paste0("../input/aptos2019-blindness-detection/train_images/", file_list[i])) # Resizing and removing original image image_res <- image %>% resize(128, 128) %>% channel(2) %>% renorm() rm(image) # Mask to separate foreground from background fundus_mask <- matrix(image_res, nrow=128, ncol=128) fundus_mask <- ifelse(fundus_mask<20, 0, 1) fundus_mask <- fundus_mask %>% as.cimg(dim(image_res)) # Image with mask, autocrop and new resize to fit in a 128x128 square image_mask <- (image_res*fundus_mask) %>% autocrop() %>% resize(128, 128) # LBPs image_mat <- image_mask %>% as.matrix(., nrow=128, ncol=128) image_lbp <- lbp(image_mat, r=1) complete_r1 <- rbind(complete_r1, t(data.matrix(prop.table(table(factor(image_lbp$lbp.u2, levels=0:58)))))) image_lbp <- lbp(image_mat, r=2) complete_r2 <- rbind(complete_r2, t(data.matrix(prop.table(table(factor(image_lbp$lbp.u2, levels=0:58)))))) } # Save patterns write.csv(complete_r1, file="complete_r1.csv", quote=FALSE, row.names=FALSE) write.csv(complete_r2, file="complete_r2.csv", quote=FALSE, row.names=FALSE) ``` # Training At this point, every training image has been converted to a 118 dimensional pattern vector. Here is configured a neural network with two 512- node dense layers to be trained with the input patterns. I know that the architecture isn't the most efficient for this problem, but it assures the training in a short number of epochs (60) with a batch size of 64 patterns each. A 20% of the training images are used for validation for this model. ```{r message=FALSE, warning=FALSE, fig.align="center"} # Training set file train <- read.csv2(file="../input/aptos2019-blindness-detection/train.csv", sep=",") # Training patterns files complete_r1 <- as.matrix(read.table(file="../input/retinal-patterns/complete_r1.csv", sep=",")) complete_r1 <- complete_r1[-1,] complete_r2 <- as.matrix(read.table(file="../input/retinal-patterns/complete_r2.csv", sep=",")) complete_r2 <- complete_r2[-1,] # Combined training patterns train_patterns <- cbind(complete_r1, complete_r2) # Classes definition for output train_output <- to_categorical(train$diagnosis, 5) # Model model <- keras_model_sequential() model %>% layer_dense(units=512, activation="relu", input_shape=c(118)) %>% layer_dropout(rate=0.1) %>% layer_dense(units=512, activation="relu") %>% layer_dropout(rate=0.1) %>% layer_dense(units=5, activation="softmax") # Compilation model %>% compile( loss="categorical_crossentropy", optimizer=optimizer_rmsprop(), metrics=c("accuracy") ) # Training history <- model %>% fit(train_patterns, train_output, epochs=60, batch_size=64, validation_split=0.2) plot(history) ``` The plot above shows how the validation loss stops descending, so we can consider this as the point to complete the training. Otherwise, the model starts overfitting. The accuracy for the validation set prediction turns out to be: ```{r message=FALSE, warning=FALSE} # Predictions predictions <- model %>% predict_classes(train_patterns) train <- cbind(train, predictions) train$predictions <- as.integer(train$predictions) # Predictions for the bottom 20% (validation set) validation <- tail(train, 730) # Accuracy accuracy(validation$diagnosis, validation$predictions) ``` And the quadratic weighted kappa: ```{r message=FALSE, warning=FALSE} # Quadratic weighted kappa ScoreQuadraticWeightedKappa(validation$diagnosis, validation$predictions) ``` It's also interesting to show the confussion matrix and the amount of predictions by diagnosis category: ```{r message=FALSE, warning=FALSE} # Confussion table table(diagnosis=validation$diagnosis, prediction=validation$predictions) # Amount of predictions by category table(prediction=validation$predictions) ``` As we can see, the classification isn't so bad taking into account the low resolution of the images transformed to a 118 dimensional vector. # Texture descriptors for test set The LBPs for the test set are computed exactly the same way as the training set were: ```{r message=FALSE, warning=FALSE} # List all training files file_list <- dir(path="../input/aptos2019-blindness-detection/test_images/", pattern="*.png") # Patterns initialization for the complete 128x128 image complete_r1 <- matrix(nrow=0, ncol=59) complete_r2 <- matrix(nrow=0, ncol=59) for (i in 1:length(file_list)) { # Load image image <- load.image(paste0("../input/aptos2019-blindness-detection/test_images/", file_list[i])) # Resizing and removing original image image_res <- image %>% resize(128, 128) %>% channel(2) %>% renorm() rm(image) # Mask to separate foreground from background fundus_mask <- matrix(image_res, nrow=128, ncol=128) fundus_mask <- ifelse(fundus_mask<20, 0, 1) fundus_mask <- fundus_mask %>% as.cimg(dim(image_res)) # Image with mask, autocrop and new resize to fit in a 128x128 square image_mask <- (image_res*fundus_mask) %>% autocrop() %>% resize(128, 128) # LBPs image_mat <- image_mask %>% as.matrix(., nrow=128, ncol=128) image_lbp <- lbp(image_mat, r=1) complete_r1 <- rbind(complete_r1, t(data.matrix(prop.table(table(factor(image_lbp$lbp.u2, levels=0:58)))))) image_lbp <- lbp(image_mat, r=2) complete_r2 <- rbind(complete_r2, t(data.matrix(prop.table(table(factor(image_lbp$lbp.u2, levels=0:58)))))) } # Merge the patterns test_patterns <- cbind(complete_r1, complete_r2) ``` # Test set prediction The trained neural network is used here to make predictions on the test set: ```{r message=FALSE, warning=FALSE} # Test set file test <- read.csv2(file="../input/aptos2019-blindness-detection/test.csv", sep=",") # Predictions diagnosis <- model %>% predict_classes(test_patterns) test <- cbind(test, diagnosis) test$diagnosis <- as.integer(test$diagnosis) ``` The model predicts the following amounts by category: ```{r message=FALSE, warning=FALSE} # Amount of predictions by category table(prediction=test$diagnosis) ``` And finally, a chi-squared test is performed to prove if the categories predicted in the training and test sets follow the same distribution, allowing us to infer the accuracy for the test set: ```{r message=FALSE, warning=FALSE} # Chi-square test between validation and test predictions pred_counts <- as.table(rbind(table(factor(validation$predictions, levels=0:4)), table(factor(test$diagnosis, levels=0:4)))) chisq.test(pred_counts) ``` The obtained p-value suggests to reject the null hypothesis. For us, this means that it doesn't give information about the accuracy of the test predictions. However, a p-value much greater than the significance level (typically 0.05) would have told us that the accuracies are likely to be similar, and ultimately the images in both data sets also resemble. Definitely, this isn't the case. ```{r message=FALSE, warning=FALSE} # Save to submission file write.csv(test, file="submission.csv", quote=FALSE, row.names=FALSE) ``` # Conclusions My solution was far from aspiring to a respectable position in the leaderboard. The chosen approach wasn't a good choice for that goal, but I was confident to obtain a 0.75 score before reducing complexity and image pre-processing. The fact is that with similar image formatting between the training and test sets and a higher resolution for the LBPs, maybe the final score could have been closer to my expectations. However, there are good news too. The accuracy for the validation set is acceptable for a such simple pre-processing and low resolution images. I got an accuracy even better than the achieved in _Identify Diabetic Retinopathy by Color Textures (Vo & Verma, 2016)_ using less patterns, and there is still room for research and improvement.